Commit b5c09390 authored by Aiden Lab's avatar Aiden Lab
Browse files

added in differing collisions treatment, taking out collisions script

parent b3aaa058

CPU/common/check.sh

+1 −1
Changes for CPU/common/check.sh: 1 added line, 1 removed line.
Original line number Diff line number Diff line
@@ -29,7 +29,7 @@ 


# Start by checking the statistics to see if they add up 

res1=`cat ${splitdir}/*norm*res*`

check1=`cat ${splitdir}/*norm*res* | awk '{s2+=$2; s3+=$3; s4+=$4; s5+=$5; s6+=$6; s7+=$7}END{if (s2 != s3+s4+s5+s6+s7){print 0}else{print 1}}'`

check1=`cat ${splitdir}/*norm*res* | awk '{s2+=$2; s3+=$3; s4+=$4; s5+=$5; s6+=$6; s7+=$7; s8+=$8}END{if (s2 != s3+s4+s5+s6+s7+s8){print 0}else{print 1}}'`

if [ $check1 -eq 0 ] || [ -z "$res1" ]

then

    echo "***! Error! The statistics do not add up. Alignment likely failed to complete on one or more files. Run relaunch_prep.sh"

CPU/common/chimeric_blacklist.awk

+208 −52
Changes for CPU/common/chimeric_blacklist.awk: 208 added lines, 52 removed lines.
Original line number Diff line number Diff line
@@ -76,7 +76,8 @@ BEGIN{ 
  count_unmapped = -1; # will count first non-group

  norm=fname"_norm.txt"

  alignable=fname"_alignable.sam"

  abnorm=fname"_abnorm.sam"

  collisions=fname"_collisions.sam"

  lowqcollisions=fname"_collisions_low_mapq.sam"

  unmapped=fname"_unmapped.sam"

  mapq0=fname"_mapq0.sam"

  dup=fname"_dup.sam" # for after duplicates have been marked
@@ -114,7 +115,6 @@ $0 !~ /^@/{ 
      mapped[j] = and(tmp[2],4) == 0;

      sum_mapped = sum_mapped+mapped[j];

    }



    if (sum_mapped < 2) {

      if (count_unmapped == -1) {

	print header > unmapped;
@@ -127,15 +127,13 @@ $0 !~ /^@/{ 
      count_unmapped++;

    }

    else {

      # find minimum mapq; if it's 0, send to mapq0 file

      # find if there are at least two mapq > 0 reads

      sum_mapq=0;

      minmapq=1;

      for (j=1; j <= count; j++) {

        split(c[j],tmp);

	if (tmp[5] > 0) {

          sum_mapq=sum_mapq+1;

	}

	else minmapq=0;

      }

      if (sum_mapq < 2) {

	if (count_mapq0 == 0) {
@@ -146,16 +144,42 @@ $0 !~ /^@/{ 
          print c[j] > mapq0;

	}

      }

      else if ((count > 4) || (count == 1)) {

      else if (count == 1) {

        # count==1 shouldn't happen; occurs with alternate aligners at times

	  if (count_abnorm == 0) {

	      print header > abnorm;

	# otherwise can signal a problem matching the fastq readnames

	if (count_lowqcollisions == 0) {

	  print header > lowqcollisions;

	}	

	print c[1] > lowqcollisions;

	count_lowqcollisions++;

      }

      else if (count > 4) {	  

	# check that MAPQ >= 10 for all reads

	mapq10=1;

	for (j=1; j <= count; j++) {

          print c[j] > abnorm;

	  split(c[j],tmp);

	  if (tmp[5] < 10) {

	    mapq10=0;

	  }

	}

	if (mapq10){

	  if (count_collisions == 0) {

	    print header > collisions;

	  }

	  for (j=1; j <= count; j++) {

	    print c[j] > collisions;

	  }

	  count_collisions++;

	}

	else {

	  if (count_lowqcollisions == 0) {

	    print header > lowqcollisions;

	  }

	  for (j=1; j <= count; j++) {

	    print c[j] > lowqcollisions;

	  }

	  count_lowqcollisions++;

	}

	count_abnorm++;

	

      }

      else {

        # count is 2,3,4.
@@ -251,26 +275,62 @@ $0 !~ /^@/{ 
	    }

	    else {

	      # abnormal

	      if (count_abnorm == 0) {

		print header > abnorm;

	      # check that MAPQ >= 10 for all reads

	      mapq10=1;

	      for (j=1; j <= count; j++) {

		split(c[j],tmp);

		if (tmp[5] < 10) {

		  mapq10=0;

		}

	      }

	      if (mapq10) {

		if (count_collisions == 0) {

		  print header > collisions;

		}

		for (j=1; j <= count; j++) {

		print c[j] > abnorm;

		  print c[j] > collisions;

		}

		count_collisions++;

	      }

	      else {

		if (count_lowqcollisions == 0) {

		  print header > lowqcollisions;

		}

		for (j=1; j <= count; j++) {

		  print c[j] > lowqcollisions;

		}

		count_lowqcollisions++;

	      }

	      count_abnorm++;

	

	    }

	  }

	  else { # count 4, not close together

	    # abnormal

	    if (count_abnorm == 0) {

	      print header > abnorm;

	    # check that MAPQ >= 10 for all reads

	    mapq10=1;

	    for (j=1; j <= count; j++) {

	      split(c[j],tmp);

	      if (tmp[5] < 10) {

		mapq10=0;

	      }

	    }

	    if (mapq10) {

	      if (count_collisions == 0) {

	        print header > collisions;

	      }

	      for (j=1; j <= count; j++) {

	      print c[j] > abnorm;

		print c[j] > collisions;

	      }

	      count_collisions++;

	    }

	    else {

	      if (count_lowqcollisions == 0) {

		print header > lowqcollisions;

	      }

	      for (j=1; j <= count; j++) {

		print c[j] > lowqcollisions;

	      }

	      count_lowqcollisions++;

	    }

	    count_abnorm++;

	

	  }

	} # if sum_mapq == 4

	else { # sum_mapq == 3
@@ -308,12 +368,31 @@ $0 !~ /^@/{ 
	  }

	  else {

	    # chimeric read with the 3 ends > 1KB apart

	    if (count_abnorm == 0) {

	      print header > abnorm;

	    # check that MAPQ >= 10 for all reads

	    mapq10=1;

	    for (j=1; j <= count; j++) {

	      split(c[j],tmp);

	      if (tmp[5] < 10) {

		mapq10=0;

	      }

	    }

	    if (mapq10) {

	      if (count_collisions == 0) {

	        print header > collisions;

	      }

	      for (j=1; j <= count; j++) {

		print c[j] > collisions;

	      }

	      count_collisions++;

	    }

	    else {

	      if (count_lowqcollisions == 0) {

		print header > lowqcollisions;

	      }

	    count_abnorm++;

	      for (j=1; j <= count; j++) {

	      print c[j] > abnorm;

		print c[j] > lowqcollisions;

	      }

	      count_lowqcollisions++;

	    }

	  }

	} # sum_mapq == 3
@@ -374,10 +453,10 @@ END{ 
    }



    if (sum_mapped < 2) {

      # unmapped

      if (count_unmapped == -1) {

	print header > unmapped;

      }

      # unmapped

      for (j=1; j <= count; j++) {

	print c[j] > unmapped;

      }
@@ -404,16 +483,42 @@ END{ 
          print c[j] > mapq0;

	}

      }

      else if ((count > 4) || (count == 1)) {

      else if (count == 1) {

        # count==1 shouldn't happen; occurs with alternate aligners at times

	if (count_abnorm == 0) {

	  print header > abnorm;

	# otherwise can signal a problem matching the fastq readnames

	if (count_lowqcollisions == 0) {

	  print header > lowqcollisions;

	}	

	print c[1] > lowqcollisions;

	count_lowqcollisions++;

      }

      else if (count > 4) {	  

	# check that MAPQ >= 10 for all reads

	mapq10=1;

	for (j=1; j <= count; j++) {

          print c[j] > abnorm;

	  split(c[j],tmp);

	  if (tmp[5] < 10) {

	    mapq10=0;

	  }

	}

	if (mapq10){

	  if (count_collisions == 0) {

	    print header > collisions;

	  }

	  for (j=1; j <= count; j++) {

	    print c[j] > collisions;

	  }

	  count_collisions++;

	}

	else {

	  if (count_lowqcollisions == 0) {

	    print header > lowqcollisions;

	  }

	  for (j=1; j <= count; j++) {

	    print c[j] > lowqcollisions;

	  }

	  count_lowqcollisions++;

	}

	count_abnorm++;

	

      }

      else {

        # count is 2,3,4.
@@ -480,10 +585,6 @@ END{ 
	      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;
@@ -513,26 +614,62 @@ END{ 
	    }

	    else {

	      # abnormal

	      if (count_abnorm == 0) {

	        print header > abnorm;

	      # check that MAPQ >= 10 for all reads

	      mapq10=1;

	      for (j=1; j <= count; j++) {

		split(c[j],tmp);

		if (tmp[5] < 10) {

		  mapq10=0;

		}

	      }

	      if (mapq10) {

		if (count_collisions == 0) {

		  print header > collisions;

		}

		for (j=1; j <= count; j++) {

		print c[j] > abnorm;

		  print c[j] > collisions;

		}

		count_collisions++;

	      }

	      else {

		if (count_lowqcollisions == 0) {

		  print header > lowqcollisions;

		}

		for (j=1; j <= count; j++) {

		  print c[j] > lowqcollisions;

		}

		count_lowqcollisions++;

	      }

	      count_abnorm++;

	

	    }

	  }

	  else { # count 4, not close together

	    # abnormal

	    if (count_abnorm == 0) {

	      print header > abnorm;

	    # check that MAPQ >= 10 for all reads

	    mapq10=1;

	    for (j=1; j <= count; j++) {

	      split(c[j],tmp);

	      if (tmp[5] < 10) {

		mapq10=0;

	      }

	    }

	    if (mapq10) {

	      if (count_collisions == 0) {

	        print header > collisions;

	      }

	      for (j=1; j <= count; j++) {

	      print c[j] > abnorm;

		print c[j] > collisions;

	      }

	      count_collisions++;

	    }

	    else {

	      if (count_lowqcollisions == 0) {

		print header > lowqcollisions;

	      }

	      for (j=1; j <= count; j++) {

		print c[j] > lowqcollisions;

	      }

	      count_lowqcollisions++;

	    }

	    count_abnorm++;

	

	  }

	} # if sum_mapq == 4

	else { # sum_mapq == 3
@@ -570,19 +707,38 @@ END{ 
	  }

	  else {

	    # chimeric read with the 3 ends > 1KB apart

	    if (count_abnorm == 0) {

	      print header > abnorm;

	    # check that MAPQ >= 10 for all reads

	    mapq10=1;

	    for (j=1; j <= count; j++) {

	      split(c[j],tmp);

	      if (tmp[5] < 10) {

		mapq10=0;

	      }

	    }

	    if (mapq10) {

	      if (count_collisions == 0) {

	        print header > collisions;

	      }

	      for (j=1; j <= count; j++) {

		print c[j] > collisions;

	      }

	      count_collisions++;

	    }

	    else {

	      if (count_lowqcollisions == 0) {

		print header > lowqcollisions;

	      }

	    count_abnorm++;

	

	      for (j=1; j <= count; j++) {

	      print c[j] > abnorm;

		print c[j] > lowqcollisions;

	      }

	      count_lowqcollisions++;

	    }

	  }

	} # sum_mapq == 3

	if (normal) {

	  # still need to check mapq0

	  if (m[read1]==0 || m[read2]==0) {

	    print "This shouldn't happen"

            # mapq0

	    for (j=1; j <= count; j++) {

	      print c[j] > mapq0;
@@ -614,5 +770,5 @@ END{ 
      } # else (mapq above 0) 

    } # else (sum_mapped >= 2)

    resfile = norm".res.txt"

    printf("%d %d %d %d %d %d\n", tottot, count_unmapped, count_reg, count_norm, count_abnorm, count_mapq0) >> resfile;

    printf("%d %d %d %d %d %d %d\n", tottot, count_unmapped, count_reg, count_norm, count_collisions, count_lowqcollisions, count_mapq0) >> resfile;

}

CPU/common/collisions.awk

deleted100755 → 0
+0 −121
Changes for CPU/common/collisions.awk: 0 added lines, 121 removed lines.
Original line number Diff line number Diff line 
#!/usr/bin/awk -f

##########

#The MIT License (MIT)

#

# Copyright (c) 2015 Aiden Lab

#

# Permission is hereby granted, free of charge, to any person obtaining a copy

# of this software and associated documentation files (the "Software"), to deal

# in the Software without restriction, including without limitation the rights

# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

# copies of the Software, and to permit persons to whom the Software is

# furnished to do so, subject to the following conditions:

#

# The above copyright notice and this permission notice shall be included in

# all copies or substantial portions of the Software.

#

#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN

#  THE SOFTWARE.

##########

# Takes a SAM file, looks for interesting abnormal chimeric reads.

# only those reads mapped to chromosomes 1-24 with MAPQ >= 10, and

# with at least a triple where the minimum distance >= 20Kb

# Juicer version 1.5

function abs(v) {

    return v<0?-v:v;

}

BEGIN{

  OFS="\t";

  mapq_thresh=-1;

  tottot = -1; # will count first non-group

}

{

  # input file is sorted by read name.  Look at read name to group 

  # appropriately

  split($1,a,"/");

  if(a[1]==prev){

    # move on to next record.  look below this block for actions that occur

    # regardless.

    count++;

  }

  else {

    # deal with read pair group

    tottot++;

    printme = 1;

    j = 0;

    for (i in c) {

      split(c[i], tmp);

      mapped[j] = tmp[3] ~ /^[1-9,X,Y,M][0-9,T]?$/ || tmp[3] ~ /^chr[1-9,X,Y,M][0-9,T]?$/;

      name[j] = tmp[1];

      str[j] = tmp[2];

      chr[j] = tmp[3];

      pos[j] = tmp[4];

      m[j] = tmp[5];

      cig[j] = tmp[6];

      seq[j] = tmp[10];

      j++;

    }

    len = j;

    bigdist = 0;

    for (j=0; j<len; j++) {

        printme = printme && mapped[j] && (m[j] >= 10);

        for (k=j+1; k<len; k++) {

            if (abs(pos[j]-pos[k]) >= 20000) {

                bigdist++;

            }

        }

    }

#    printme = printme && (bigdist >= 3);

    if (printme) {

      for (j = 0; j < len; j++) {

          printf("%s %d %s %d %d %s %s ", name[j], str[j], chr[j], pos[j], m[j], cig[j], seq[j]);

      }

      printf("\n");

    }



    # reset variables

    delete c;

    count=1;

    prev=a[1];

  }

  # these happen no matter what, after the above processing

  c[count] = $0;

}

END {

    printme = 1;

    j = 0;

    for (i in c) {

      split(c[i], tmp);

      mapped[j] = tmp[3] ~ /^[1-9,X,Y,M][0-9,T]?$/ || tmp[3] ~ /^chr[1-9,X,Y,M][0-9,T]?$/;

      name[j] = tmp[1];

      str[j] = tmp[2];

      chr[j] = tmp[3];

      pos[j] = tmp[4];

      m[j] = tmp[5];

      cig[j] = tmp[6];

      seq[j] = tmp[10];

      j++;

    }

    len = j;

    bigdist = 0;

    for (j=0; j<len; j++) {

        printme = printme && mapped[j] && (m[j] >= 10);

        for (k=j+1; k<len; k++) {

            if (abs(pos[j]-pos[k]) >= 20000) {

                bigdist++;

            }

        }

    }

#    printme = printme && (bigdist >= 3);

    if (printme) {

      for (j = 0; j < len; j++) {

          printf("%s %d %s %d %d %s %s ", name[j], str[j], chr[j], pos[j], m[j], cig[j], seq[j]);

      }

      printf("\n");

    }

}

CPU/common/stats_sub.awk

+10 −8
Changes for CPU/common/stats_sub.awk: 10 added lines, 8 removed lines.
Original line number Diff line number Diff line
@@ -25,14 +25,16 @@ 
# Script to read in individual split outputs and put stats together

# Juicer version 1.5

{

a1+=$2; # total reads

a2+=$3; # normal

a3+=$4; # chimeric paired

a4+=$5; # chimeric ambiguous

a5+=$6; # unmapped

a6+=$1; # ligations 

a7+=$7; # mapq0

tot+=$2; # total reads

unm+=$3; # unmapped

norm+=$4; # normal paired

chim+=$5; # chimeric paired

coll+=$6; # collisions

lowc+=$7; # low mapq collisions 

mapq0+=$8; # mapq0

lig+=$1; #ligations

}

END{

    printf("Sequenced Read Pairs:  %'d\n Normal Paired: %'d (%0.2f%)\n Chimeric Paired: %'d (%0.2f%)\n Chimeric Ambiguous: %'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", a1, a3, a3*100/a1, a4, a4*100/a1, a5, a5*100/a1, a2, a2*100/a1, a7, a7*100/a1, a6, a6*100/a1, a3+a4, (a3+a4)*100/a1);

    if (tot==0) tot=1;

    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", 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);

}

CPU/juicer.sh

+17 −12
Changes for CPU/juicer.sh: 17 added lines, 12 removed lines.
Original line number Diff line number Diff line
@@ -60,8 +60,7 @@ 
shopt -s extglob

juicer_version="1.5.6" 

### LOAD BWA AND SAMTOOLS

bwa_cmd="/home/neva/smalltest/bwa/bwa"

#bwa_cmd="bwa"

bwa_cmd="/home/neva/bwa-0.7.17/bwa"

# fastq files should look like filename_R1.fastq and filename_R2.fastq 

# if your fastq files look different, change this value

read1str="_R1"
@@ -71,7 +70,7 @@ read2str="_R2" 


# Juicer directory, contains scripts/, references/, and restriction_sites/

# can also be set in options via -D

juiceDir="/opt/juicer"

juiceDir="/aidenlab"

# top level directory, can also be set in options

topDir=$(pwd)

# restriction enzyme, can also be set in options
@@ -360,7 +359,7 @@ then 
	export LC_ALL=C

        # call chimeric_blacklist.awk to deal with chimeric reads; 

        # sorted file is sorted by read name at this point

	touch $name${ext}_abnorm.sam $name${ext}_unmapped.sam $name${ext}_mapq0.sam

	touch $name${ext}_collisions.sam $name${ext}_collisions_low_mapq.sam $name${ext}_collisions $name${ext}_unmapped.sam $name${ext}_mapq0.sam

	# chimeric takes in $name$ext

	awk -v "fname"=$name${ext} -f ${juiceDir}/scripts/common/chimeric_blacklist.awk $name$ext.sam

	if [ $? -ne 0 ]
@@ -368,7 +367,6 @@ then 
            echo "***! Failure during chimera handling of $name${ext}"

            exit 1

	fi	

	#samtools fixmate tmp_se.bam tmp_se_fixmate.bam

        # if any normal reads were written, find what fragment they correspond

	# to and store that

	if [ -e "$name${ext}_norm.txt" ] && [ "$site" != "none" ]
@@ -388,10 +386,16 @@ then 
	fi                              



	# convert sams to bams and delete the sams

	samtools view -hb $name${ext}_abnorm.sam > $name${ext}_abnorm.bam

	samtools view -hb $name${ext}_collisions.sam > $name${ext}_collisions.bam

	if [ $? -ne 0 ]

	then

            echo "***! Failure during bam write of $name${ext}_abnorm.sam"

            echo "***! Failure during bam write of $name${ext}_collisions.sam"

            exit 1

	fi	

	samtools view -hb $name${ext}_collisions_low_mapq.sam > $name${ext}_collisions_low_mapq.bam

	if [ $? -ne 0 ]

	then

            echo "***! Failure during bam write of $name${ext}_collisions_low_mapq.sam"

            exit 1

	fi	

	samtools view -hb $name${ext}_unmapped.sam > $name${ext}_unmapped.bam
@@ -413,7 +417,7 @@ then 
            exit 1

	fi	

	# remove all sams EXCEPT alignable, which we need for deduping

	rm $name${ext}_abnorm.sam $name${ext}_unmapped.sam $name${ext}_mapq0.sam

	rm $name${ext}_collisions.sam $name${ext}_collisions_low_mapq.sam $name${ext}_unmapped.sam $name${ext}_mapq0.sam



        # sort by chromosome, fragment, strand, and position

	sort -T $tmpdir -k2,2d -k6,6d -k4,4n -k8,8n -k1,1n -k5,5n -k3,3n $name${ext}.frag.txt > $name${ext}.sort.txt
@@ -477,7 +481,7 @@ then 
	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

#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}/alignable.bam ${splitdir}/*_alignable.bam
@@ -487,15 +491,16 @@ then 
    # combine bams

    if [ `ls -1 $splitdir/*_abnorm.bam 2>/dev/null | wc -l ` -gt 1 ]

    then

	samtools merge -n $outputdir/abnormal.bam $splitdir/*_abnorm.bam 

	samtools merge -n $outputdir/collisions.bam $splitdir/*_collisions.bam 

	samtools merge -n $outputdir/collisions_low_mapq.bam $splitdir/*_collisions_low_mapq.bam 

	samtools merge -n $outputdir/unmapped.bam $splitdir/*_unmapped.bam 

	samtools merge -n $outputdir/mapq0.bam $splitdir/*_mapq0.bam

    else

	cat $splitdir/*_abnorm.bam > $outputdir/abnormal.bam

	cat $splitdir/*_collisions.bam > $outputdir/collisions.bam

	cat $splitdir/*_collisions_low_mapq.bam > $outputdir/collisions_low_mapq.bam

	cat $splitdir/*_unmapped.bam > $outputdir/unmapped.bam

	cat $splitdir/*_mapq0.bam > $outputdir/mapq0.bam

    fi

    samtools view $outputdir/abnormal.bam | awk -f ${juiceDir}/scripts/common/collisions.awk  > $outputdir/collisions.txt

fi



if [ -n "$deduponly" ]

CRA Git | Maintained and supported by SUSTech CRA and CCSE