further bug fixes to encode pipeline (b3aaa058) · Commits · Chaos / chaos.juicer

CPU/common/chimeric_blacklist.awk

+45 −25

Changes for CPU/common/chimeric_blacklist.awk: 45 added lines, 25 removed lines.
Original line number	Diff line number	Diff line
		@@ -79,23 +79,18 @@ BEGIN{
		abnorm=fname"_abnorm.sam"
		unmapped=fname"_unmapped.sam"
		mapq0=fname"_mapq0.sam"
		dup=fname"_dup.sam" # for after duplicates have been marked
		linenum=0
		}
		$0 ~ /^@/{
		# print SAM header to SAM files
		print > alignable;
		linenum++;
		print > abnorm;
		print > unmapped;
		print > mapq0;
		header = header""$0"\n";
		}
		$0 !~ /^@/{
		if (tottot == -1) {
		print "@PG", "ID:Juicer", "VN:1.6" > alignable;
		linenum++;
		print "@PG", "ID:Juicer", "VN:1.6" > abnorm;
		print "@PG", "ID:Juicer", "VN:1.6" > unmapped;
		print "@PG", "ID:Juicer", "VN:1.6" > mapq0;
		header = header"@PG\tID:Juicer\tVN:1.6";
		}
		# input file is sorted by read name. Look at read name to group
		# appropriately
		@@ -121,6 +116,9 @@ $0 !~ /^@/{
		}

		if (sum_mapped < 2) {
		if (count_unmapped == -1) {
		print header > unmapped;
		}
		# unmapped
		for (j=1; j <= count; j++) {
		print c[j] > unmapped;
		@@ -140,7 +138,9 @@ $0 !~ /^@/{
		else minmapq=0;
		}
		if (sum_mapq < 2) {
		# if (sum_mapq < 2 \|\| minmapq==0) {
		if (count_mapq0 == 0) {
		print header > mapq0;
		}
		count_mapq0++;
		for (j=1; j <= count; j++) {
		print c[j] > mapq0;
		@@ -148,6 +148,9 @@ $0 !~ /^@/{
		}
		else if ((count > 4) \|\| (count == 1)) {
		# count==1 shouldn't happen; occurs with alternate aligners at times
		if (count_abnorm == 0) {
		print header > abnorm;
		}
		for (j=1; j <= count; j++) {
		print c[j] > abnorm;
		}
		@@ -219,10 +222,6 @@ $0 !~ /^@/{
		cigloc = substr(tmp[6],where,RLENGTH-1) + 0;
		pos[j] = pos[j] + cigloc;
		}
		# Mitochrondria loops around
		if (chr[j] ~ /MT/ && pos[j] >= 16569) {
		pos[j] = pos[j] - 16569;
		}
		} # else if str[j]==16
		} # end for loop, now have all info about each read end
		read1=0;
		@@ -252,6 +251,9 @@ $0 !~ /^@/{
		}
		else {
		# abnormal
		if (count_abnorm == 0) {
		print header > abnorm;
		}
		for (j=1; j <= count; j++) {
		print c[j] > abnorm;
		}
		@@ -261,6 +263,9 @@ $0 !~ /^@/{
		}
		else { # count 4, not close together
		# abnormal
		if (count_abnorm == 0) {
		print header > abnorm;
		}
		for (j=1; j <= count; j++) {
		print c[j] > abnorm;
		}
		@@ -303,8 +308,10 @@ $0 !~ /^@/{
		}
		else {
		# chimeric read with the 3 ends > 1KB apart
		if (count_abnorm == 0) {
		print header > abnorm;
		}
		count_abnorm++;

		for (j=1; j <= count; j++) {
		print c[j] > abnorm;
		}
		@@ -321,15 +328,13 @@ $0 !~ /^@/{
		count_mapq0++;
		}
		else {
		if (count_reg == 0 && count_norm == 0) {
		print header > alignable;
		}
		if (sum_mapped == 2 && count==2) {
		count_reg++;
		}
		else count_norm++;
		# or do we want to print the whole read group including chimeric?
		#print c[read1] > alignable;
		#linenum++;
		#print c[read2] > alignable;
		#linenum++;
		prevlinenum=linenum;
		for (j=1; j <= count; j++) {
		print c[j] > alignable;
		@@ -370,6 +375,9 @@ END{

		if (sum_mapped < 2) {
		# unmapped
		if (count_unmapped == -1) {
		print header > unmapped;
		}
		for (j=1; j <= count; j++) {
		print c[j] > unmapped;
		}
		@@ -388,7 +396,9 @@ END{
		else minmapq=0;
		}
		if (sum_mapq < 2) {
		# if (sum_mapq < 2 \|\| minmapq == 0) {
		if (count_mapq0 == 0) {
		print header > mapq0;
		}
		count_mapq0++;
		for (j=1; j <= count; j++) {
		print c[j] > mapq0;
		@@ -396,6 +406,9 @@ END{
		}
		else if ((count > 4) \|\| (count == 1)) {
		# count==1 shouldn't happen; occurs with alternate aligners at times
		if (count_abnorm == 0) {
		print header > abnorm;
		}
		for (j=1; j <= count; j++) {
		print c[j] > abnorm;
		}
		@@ -500,6 +513,9 @@ END{
		}
		else {
		# abnormal
		if (count_abnorm == 0) {
		print header > abnorm;
		}
		for (j=1; j <= count; j++) {
		print c[j] > abnorm;
		}
		@@ -509,6 +525,9 @@ END{
		}
		else { # count 4, not close together
		# abnormal
		if (count_abnorm == 0) {
		print header > abnorm;
		}
		for (j=1; j <= count; j++) {
		print c[j] > abnorm;
		}
		@@ -551,6 +570,9 @@ END{
		}
		else {
		# chimeric read with the 3 ends > 1KB apart
		if (count_abnorm == 0) {
		print header > abnorm;
		}
		count_abnorm++;

		for (j=1; j <= count; j++) {
		@@ -568,15 +590,13 @@ END{
		count_mapq0++;
		}
		else {
		if (count_reg == 0 && count_norm == 0) {
		print header > alignable;
		}
		if (sum_mapped == 2 && count==2) {
		count_reg++;
		}
		else count_norm++;
		# or do we want to print the whole read group including chimeric?
		#print c[read1] > alignable;
		#linenum++;
		#print c[read2] > alignable;
		#linenum++;
		prevlinenum=linenum;
		for (j=1; j <= count; j++) {
		print c[j] > alignable;

CPU/common/cleanup.sh

+5 −5

Changes for CPU/common/cleanup.sh: 5 added lines, 5 removed lines.
Original line number	Diff line number	Diff line
		@@ -25,22 +25,22 @@
		# Script to clean up big repetitive files and zip fastqs. Run after you are
		# sure the pipeline ran successfully. Run from top directory (HIC001 e.g.).
		# Juicer version 1.5
		total=`ls -l aligned/merged_sort.txt \| awk '{print $5}'`
		total2=`ls -l aligned/merged_nodups.txt aligned/dups.txt aligned/opt_dups.txt \| awk '{sum = sum + $5}END{print sum}'`
		if [ $total -eq $total2 ]
		then
		rm aligned/merged_sort.txt
		rm -r splits
		rm -rf splits
		testname=$(ls -l fastq \| awk 'NR==1{print $9}')
		if [ "${testname: -5}" == ".fastq" ]
		then
		for i in fastq/*.fastq
		do
		gzip $i
		done
		fi
		gzip aligned/merged_nodups.txt
		gzip aligned/dups.txt
		gzip aligned/opt_dups.txt
		gzip aligned/abnormal.sam
		gzip aligned/collisions.txt
		gzip aligned/unmapped.sam
		else
		echo "Problem: The sum of merged_nodups and the dups files is not the same size as merged_sort.txt"
		echo "Did NOT clean up";

CPU/common/juicer_postprocessing.sh

+3 −3

Changes for CPU/common/juicer_postprocessing.sh: 3 added lines, 3 removed lines.
Original line number	Diff line number	Diff line
		@@ -67,7 +67,7 @@ fi

		echo -e "${juicer_tools_path} is post-processing Hi-C for ${genomeID}\nData read from ${hic_file_path}.\nMotifs read from ${bed_file_dir}\n"
		echo -e "ARROWHEAD:\n"
		${juicer_tools_path} arrowhead ${hic_file_path} ${hic_file_path%.*}"_contact_domains.txt"
		${juicer_tools_path} arrowhead ${hic_file_path} ${hic_file_path%.*}"_contact_domains"
		if [ $? -ne 0 ]; then
		echo "***! Problem while running Arrowhead";
		exit 1
		@@ -75,7 +75,7 @@ fi
		echo -e "\nHiCCUPS:\n"
		if hash nvcc 2>/dev/null
		then
		${juicer_tools_path} hiccups ${hic_file_path} ${hic_file_path%.*}"_loops.txt"
		${juicer_tools_path} hiccups ${hic_file_path} ${hic_file_path%.*}"_loops"
		if [ $? -ne 0 ]; then
		echo "***! Problem while running HiCCUPS";
		exit 1
		@@ -87,7 +87,7 @@ fi
		if [ -f ${hic_file_path%.*}"_loops.txt" ]
		then
		echo -e "\nAPA:\n"
		${juicer_tools_path} apa ${hic_file_path} ${hic_file_path%.*}"_loops.txt" "apa_results"
		${juicer_tools_path} apa ${hic_file_path} ${hic_file_path%.*}"_loops" "apa_results"
		## Check that bed folder exists
		if [ ! -e "${bed_file_dir}" ]; then
		echo "***! Can't find folder ${bed_file_dir}";

CPU/juicer.sh

+31 −19

Changes for CPU/juicer.sh: 31 added lines, 19 removed lines.
Original line number	Diff line number	Diff line
		@@ -83,15 +83,16 @@ genomeID="hg19"
		shortreadend=0
		# description, default empty
		about=""
		nofrag=0
		nofrag=1

		## TODO Change usage to be correct and print nicely
		## Read arguments
		usageHelp="Usage: ${0##*/} [-g genomeID] [-d topDir] [-s site] [-a about] [-R end]\n [-S stage] [-p chrom.sizes path] [-y restriction site file]\n [-z reference genome file] [-D Juicer scripts directory]\n [-b ligation] [-t threads] [-r] [-h] [-x]"
		genomeHelp="* [genomeID] must be defined in the script, e.g. \"hg19\" or \"mm10\" (default \n \"$genomeID\"); alternatively, it can be defined using the -z command"
		dirHelp="* [topDir] is the top level directory (default\n \"$topDir\")\n [topDir]/fastq must contain the fastq files\n [topDir]/splits will be created to contain the temporary split files\n [topDir]/aligned will be created for the final alignment"
		siteHelp="* [site] must be defined in the script, e.g. \"HindIII\" or \"MboI\" \n (default \"$site\")"
		aboutHelp="* [about]: enter description of experiment, enclosed in single quotes"
		stageHelp="* [stage]: must be one of \"merge\", \"dedup\", \"final\", \"postproc\", or \"early\".\n -Use \"merge\" when alignment has finished but the merged_sort file has not\n yet been created.\n -Use \"dedup\" when the files have been merged into merged_sort but\n merged_nodups has not yet been created.\n -Use \"final\" when the reads have been deduped 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"
		stageHelp="* [stage]: must be one of \"merge\", \"dedup\", \"final\", \"postproc\", \"early\", \"alignonly\", .\n -Use \"merge\" when alignment has finished but the merged_sort file has not\n yet been created.\n -Use \"dedup\" when the files have been merged into merged_sort but\n merged_nodups has not yet been created.\n -Use \"final\" when the reads have been deduped 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"
		pathHelp="* [chrom.sizes path]: enter path for chrom.sizes file"
		siteFileHelp="* [restriction site file]: enter path for restriction site file (locations of\n restriction sites in genome; can be generated with the script\n misc/generate_site_positions.py)"
		scriptDirHelp="* [Juicer scripts directory]: set the Juicer directory,\n which should have scripts/ references/ and restriction_sites/ underneath it\n (default ${juiceDir})"
		@@ -119,7 +120,7 @@ printHelpAndExit() {
		exit "$1"
		}

		while getopts "d:g:a:hs:p:y:z:S:D:xt:b:" opt; do
		while getopts "d:g:a:hs:p:y:z:S:D:ft:b:" opt; do
		case $opt in
		g) genomeID=$OPTARG ;;
		h) printHelpAndExit 0;;
		@@ -131,7 +132,7 @@ while getopts "d:g:a:hs:p:y:z:S:D:xt:b:" opt; do
		z) refSeq=$OPTARG ;;
		S) stage=$OPTARG ;;
		D) juiceDir=$OPTARG ;;
		x) nofrag=1 ;; #no fragment maps
		f) nofrag=0 ;; #include fragment maps
		b) ligation=$OPTARG ;;
		t) threads=$OPTARG ;;
		[?]) printHelpAndExit 1;;
		@@ -147,6 +148,8 @@ then
		final) final=1 ;;
		postproc) postproc=1 ;;
		alignonly) alignonly=1 ;;
		mergeonly) mergeonly=1 ;;
		deduponly) deduponly=1 ;;
		*) echo "$usageHelp"
		echo "$stageHelp"
		exit 1
		@@ -161,7 +164,6 @@ then
		mm10) refSeq="${juiceDir}/references/Mus_musculus_assembly10.fasta";;
		hg38) refSeq="${juiceDir}/references/Homo_sapiens_assembly38.fasta";;
		hg19) refSeq="${juiceDir}/references/Homo_sapiens_assembly19.fasta";;

		*) echo "$usageHelp"
		echo "$genomeHelp"
		exit 1
		@@ -257,13 +259,13 @@ else
		fi

		## Create output directory, only if not in dedup, final, or postproc stages
		if [[ -d "$outputdir" && -z "$final" && -z "$dedup" && -z "$postproc" ]]
		if [[ -d "$outputdir" && -z "$final" && -z "$dedup" && -z "$postproc" && -z "$deduponly" ]]
		then
		echo "***! Move or remove directory \"$outputdir\" before proceeding."
		echo "***! Type \"juicer.sh -h \" for help"
		exit 1
		else
		if [[ -z "$final" && -z "$dedup" && -z "$postproc" ]]; then
		if [[ -z "$final" && -z "$dedup" && -z "$postproc" && -z "$deduponly" ]]; then
		mkdir "$outputdir" \|\| { echo "***! Unable to create ${outputdir}, check permissions." ; exit 1; }
		fi
		fi
		@@ -276,7 +278,7 @@ else
		fi

		## Create temporary directory, used for sort later
		if [ ! -d "$tmpdir" ] && [ -z "$final" ] && [ -z "$dedup" ] && [ -z "$postproc" ]; then
		if [ ! -d "$tmpdir" ] && [ -z "$final" ] && [ -z "$dedup" ] && [ -z "deduponly" ] && [ -z "$postproc" ]; then
		mkdir "$tmpdir"
		chmod 777 "$tmpdir"
		fi
		@@ -314,7 +316,7 @@ echo "$0 $@" >> $headfile
		## ALIGN FASTQ AS SINGLE END, SORT BY READNAME, HANDLE CHIMERIC READS

		## Not in merge, dedup, final, or postproc stage, i.e. need to align files.
		if [ -z $merge ] && [ -z $final ] && [ -z $dedup ] && [ -z $postproc ]
		if [ -z $merge ] && [ -z $mergeonly ] && [ -z $final ] && [ -z $dedup ] && [ -z $deduponly ] && [ -z $postproc ]
		then
		echo -e "(-: Aligning files matching $fastqdir\n to genome $genomeID with site file $site_file"
		if [ ! $splitdirexists ]
		@@ -345,7 +347,7 @@ then
		source ${juiceDir}/scripts/common/countligations.sh

		# Align reads
		echo "Running command $bwa_cmd mem $threadstring $refSeq $name1$ext $name2$ext > $name$ext.sam"
		echo "Running command $bwa_cmd mem -SP5M $threadstring $refSeq $name1$ext $name2$ext > $name$ext.sam"
		$bwa_cmd mem -SP5M $threadstring $refSeq $name1$ext $name2$ext > $name$ext.sam
		if [ $? -ne 0 ]
		then
		@@ -431,7 +433,7 @@ then
		fi

		#MERGE SORTED AND ALIGNED FILES
		if [ -z $final ] && [ -z $dedup ] && [ -z $postproc ]
		if [ -z $final ] && [ -z $dedup ] && [ -z $deduponly ] && [ -z $postproc ]
		then
		if [ -d $donesplitdir ]
		then
		@@ -443,9 +445,15 @@ then
		exit 1
		else
		echo "(-: Finished sorting all sorted files into a single merge."
		rm -r ${tmpdir}
		rm -rf ${tmpdir}
		fi
		fi

		if [ -n "$mergeonly" ]
		then
		exit 0
		fi

		#REMOVE DUPLICATES
		if [ -z $final ] && [ -z $postproc ]
		then
		@@ -466,14 +474,15 @@ then
		# this could end up being computationally expensive memory-wise
		# could also consider sed solution, no hashtable required
		# could rewrite this to not need readname mod, but need filename
		awk 'FNR==NR{a[$1];next}$0 ~ /^@/{print}(FNR in a){print}' "${i}" "${j}_alignable.sam" > "${j}_alignable_dedup.sam"
		samtools view -hb "${j}_alignable_dedup.sam" > "${j}_alignable_dedup.bam"
		awk 'BEGIN{OFS="\t"}FNR==NR{for (i=$1; i<=$2; i++){a[i];} next}(!(FNR in a) && $1 !~ /^@/){$2=or($2,1024)}{print}' "${i}" "${j}_alignable.sam" > "${j}_alignable_dedup.sam"
		samtools view -hb "${j}_alignable_dedup.sam" > "${j}_alignable.bam"
		done
		if [ `ls -1 $splitdir/*_alignable_dedup.bam 2>/dev/null \| wc -l ` -gt 1 ]
		#hi neva start here. alignable.bam should be marked with duplicates after this point
		if [ `ls -1 $splitdir/*_alignable.bam 2>/dev/null \| wc -l ` -gt 1 ]
		then
		samtools merge -n ${outputdir}/merged_nodups.bam ${splitdir}/*_alignable_dedup.bam
		samtools merge -n ${outputdir}/alignable.bam ${splitdir}/*_alignable.bam
		else
		cat ${splitdir}/*_alignable_dedup.bam > ${outputdir}/merged_nodups.bam
		cat ${splitdir}/*_alignable.bam > ${outputdir}/alignable.bam
		fi
		# combine bams
		if [ `ls -1 $splitdir/*_abnorm.bam 2>/dev/null \| wc -l ` -gt 1 ]
		@@ -481,16 +490,19 @@ then
		samtools merge -n $outputdir/abnormal.bam $splitdir/*_abnorm.bam
		samtools merge -n $outputdir/unmapped.bam $splitdir/*_unmapped.bam
		samtools merge -n $outputdir/mapq0.bam $splitdir/*_mapq0.bam
		samtools merge -n $outputdir/alignable.bam $splitdir/*_alignable.bam
		else
		cat $splitdir/*_abnorm.bam > $outputdir/abnormal.bam
		cat $splitdir/*_unmapped.bam > $outputdir/unmapped.bam
		cat $splitdir/*_mapq0.bam > $outputdir/mapq0.bam
		cat $splitdir/*_alignable.bam > $outputdir/alignable.bam
		fi
		samtools view $outputdir/abnormal.bam \| awk -f ${juiceDir}/scripts/common/collisions.awk > $outputdir/collisions.txt
		fi

		if [ -n "$deduponly" ]
		then
		exit 0
		fi

		#CREATE HIC FILES
		# if early exit, we stop here, once the merged_nodups.txt file is created.
		if [ -z "$earlyexit" ]

Admin message