Merge branch 'encode' of https://github.com/aidenlab/juicer into encode

# Juicer version 1.5

$1=="Sequenced"{

gsub(/,/,"",$4);

 total=total+$4;

 tot=tot+$4;

}

$1=="Unmapped:"{

gsub(/,/,"",$2);

 unmapped=unmapped+$2

 unm=unm+$2

}

$1=="Normal"{

gsub(/,/,"",$3);

 regular=regular+$3}

 norm=norm+$3}

$1=="Chimeric" && $2=="Paired:"{

gsub(/,/,"",$3);

 normal=normal+$3}

$2=="Ambiguous:"{

 chim=chim+$3}

$1=="MAPQ" && $2=="0:"{

gsub(/,/,"",$3);

 abnorm=abnorm+$3}

$1=="Alignable"{

 mapq0=mapq0+$3}

$1=="Collisions:"{

gsub(/,/,"",$2);

 coll=coll+$2}

$3=="Collisions:"{

gsub(/,/,"",$4);

 lowc=lowc+$4}

$1=="Ligation" {

gsub(/,/,"",$4);

 alignable=alignable+$4}

 lig=lig+$4}

$1=="Unique" {

gsub(/,/,"",$3);

 dedup=dedup+$3}

gsub(/,/,"",$3);

 optdups=optdups+$3}

END{

 printf("%s %'d\n", "Sequenced Read Pairs:", total);

 printf(" %s %'d (%0.2f%)\n", "Normal Paired:", regular, regular*100/total);

 printf(" %s %'d (%0.2f%)\n", "Chimeric Paired:", normal, normal*100/total);

 printf(" %s %'d (%0.2f%)\n", "Chimeric Ambiguous:", abnorm, abnorm*100/total);

 printf(" %s %'d (%0.2f%)\n", "Unmapped:", unmapped, unmapped*100/total);

 printf(" %s %'d (%0.2f%)\n", "Alignable (Normal+Chimeric Paired):", alignable, alignable*100/total);

 printf("%s %'d\n", "Unique Reads:", dedup);

 printf("%s %'d\n", "PCR Duplicates:", dups);

  printf("%s %'d\n", "Optical Duplicates:", optdups);

  printf("Sequenced Read Pairs:  %'d\n Normal Paired: %'d (%0.2f%)\n Chimeric Paired: %'d (%0.2f%)\n Collisions: %'d (%0.2f%)\n Low MAPQ Collisions: %'d (%0.2f%)\n Unmapped: %'d (%0.2f%)\n MAPQ 0: %'d (%0.2f%)\n Ligation Motif Present: %'d (%0.2f%)\nAlignable (Normal+Chimeric Paired): %'d (%0.2f%)\n Unique Reads: %'d\n PCR Duplicates: %'d\n Optical Duplicates: %'d\n", tot, norm, norm*100/tot, chim, chim*100/tot, coll, coll*100/tot, lowc, lowc*100/tot, unm, unm*100/tot, mapq0, mapq0*100/tot, lig, lig*100/tot, norm+chim, (norm+chim)*100/tot, dedup, dups, optdups );

}

juiceDir="/aidenlab"

usageHelp="Usage: ${0##*/} -g genomeID [-d topDir] [-s site] [-h]"

usageHelp="Usage: ${0##*/} -g genomeID [-d topDir] [-s site] [-b ligation] [-h]"

genomeHelp="   genomeID must be defined in the script, e.g. \"hg19\" or \"mm10\" (default \"$genomeID\")"

dirHelp="   [topDir] is the top level directory (default \"$topDir\") and must contain links to all merged_nodups files underneath it"

siteHelp="   [site] must be defined in the script, e.g.  \"HindIII\" or \"MboI\" (default \"$site\"); alternatively, this can be the restriction site file"

ligationHelp="* [ligation junction]: use this string when counting ligation junctions"

stageHelp="* [stage]: must be one of \"final\", \"postproc\", or \"early\".\n    -Use \"final\" when the reads have been combined into merged_nodups but the\n     final stats and hic files have not yet been created.\n    -Use \"postproc\" when the hic files have been created and only\n     postprocessing feature annotation remains to be completed.\n    -Use \"early\" for an early exit, before the final creation of the stats and\n     hic files"

excludeHelp="   -x: exclude fragment-delimited maps from Hi-C mega map (will run much faster)"

helpHelp="   -h: print this help and exit"

    echo "$genomeHelp"

    echo "$dirHelp"

    echo "$siteHelp"

    echo "$ligationHelp"

    echo "$stageHelp"

    echo "$excludeHelp"

    echo "$helpHelp"

    exit "$1"

}

while getopts "d:g:hxs:S:" opt; do

while getopts "d:g:hxs:S:b:" opt; do

    case $opt in

	g) genomeID=$OPTARG ;;

	h) printHelpAndExit 0;;

	s) site=$OPTARG ;;

	x) exclude=1 ;;

	S) stage=$OPTARG ;;

	b) ligation=$OPTARG ;;

	[?]) printHelpAndExit 1;;

    esac

done

## Set ligation junction based on restriction enzyme

if [ -z "$ligation" ]

then

    case $site in

	HindIII) ligation="AAGCTAGCTT";;

	DpnII) ligation="GATCGATC";;

	    echo "$site not listed as recognized enzyme, so trying it as site file."

	    echo "Ligation junction is undefined";;

    esac

fi

if [ -z "$site_file" ]

then

#!/usr/bin/perl

##########

#The MIT License (MIT)

#

# Check arguments

getopts('s:l:o:q:h');

my $site_file;

my $site_file = "/projects/ea14/juicer/restriction_sites/hg19_DpnII.txt";

my $ligation_junction = "GATCGATC";

my $stats_file = "stats.txt";

my $mapq_threshold = 1;

  print STDOUT "No input file specified, reading from input stream\n";

}

# Remove parenthesis

$ligation_junction =~ s/\(//;

$ligation_junction =~ s/\)//;

# with OR | symbol, this won't work, need to explicitly fix

my $dangling_junction = substr $ligation_junction, length($ligation_junction)/2;

# Global variables for calculating statistics

        }

      }

      # read pair contains ligation junction

      if ($record[10] =~ m/$ligation_junction/ || $record[13] =~ m/$ligation_junction/) {

      if ($record[10] =~ m/($ligation_junction)/ || $record[13] =~ m/($ligation_junction)/) {

        $ligation++;

      }

    }

}

printf FILE "%0.2f\%)\n", $intra_fragment*100/$unique; 

print FILE "Below MAPQ Threshold: " . commify($under_mapq);

if ($seq == 1) {

  printf FILE " (%0.2f\% / ", $under_mapq*100/$reads; 

}

else {

  print FILE "(";

}

printf FILE "%0.2f\%)\n", $under_mapq*100/$unique; 

print FILE "Hi-C Contacts: " . commify($total_current);

if ($seq == 1) {

  printf FILE " (%0.2f\% / ", $total_current*100/$reads; 

else {

  print FILE "(";

}

printf FILE "%0.2f\%)\n", $ligation*100/$unique; 

if ($five_prime_end + $three_prime_end > 0) {