Commit 66e9807c authored by nchernia's avatar nchernia
Browse files

added new handling of four way reads to END block

parent a55bfe32

CPU/common/chimeric_blacklist.awk

+70 −10
Changes for CPU/common/chimeric_blacklist.awk: 70 added lines, 10 removed lines.
Original line number Diff line number Diff line
@@ -355,15 +355,15 @@ print str[1],chr[1],pos[1],str[0],chr[0],pos[0],m[1],cigarstr[1],seq[1],m[0],cig 
END{

    # deal with read pair group

    tottot++;

  if (count==3) {

    if (count==3 || count==4) {

      # chimeric read

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

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

	split(c[j], tmp);

	split(tmp[1],readname,"/");

	read[j] = readname[2];

	name[j] = tmp[1];

	

      # strand

	# strand; Bit 16 set means reverse strand

	str[j] = and(tmp[2],16);

	# chromosome

	chr[j] = tmp[3];
@@ -376,6 +376,7 @@ END{ 
	# sequence

	seq[j] = tmp[10];

        #qual[j] = tmp[11];

        

	# get rid of soft clipping to know correct position

	if (str[j] == 0 && tmp[6] ~/^[0-9]+S/) {

	  split(tmp[6], cigar, "S");
@@ -384,6 +385,7 @@ END{ 
	    pos[j] = 1;

	  }

	}

	# get rid of hard clipping to know correct position

	else if (str[j] == 0 && tmp[6] ~/^[0-9]+H/) {

	  split(tmp[6], cigar, "H");

	  pos[j] = pos[j] - cigar[1];
@@ -409,19 +411,64 @@ END{ 
	    cigloc = substr(tmp[6],where,RLENGTH-1) + 0;

	    pos[j] = pos[j] + cigloc;

	  }

	if (tmp[6] ~ /[0-9]+H$/) {

          else if (tmp[6] ~ /[0-9]+H$/) {

	    where = match(tmp[6],/[0-9]+H$/);

	    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;

	}

	}

        # blacklist - if 3rd bit set (=4) it means unmapped

        mapped[j] = and(tmp[2],4) == 0; 

      }

      if (count == 4) {

        # looking for A/B...A/B

	dist[12] = abs(chr[1]-chr[2])*10000000 + abs(pos[1]-pos[2]);

	dist[13] = abs(chr[1]-chr[3])*10000000 + abs(pos[1]-pos[3]);

	dist[14] = abs(chr[1]-chr[4])*10000000 + abs(pos[1]-pos[4]);

	dist[23] = abs(chr[2]-chr[3])*10000000 + abs(pos[2]-pos[3]);

	dist[24] = abs(chr[2]-chr[4])*10000000 + abs(pos[2]-pos[4]);

	dist[34] = abs(chr[3]-chr[4])*10000000 + abs(pos[3]-pos[4]);

	

	# A/B, A/B

	if ((dist[13] < 1000 && dist[24] < 1000) || (dist[14] < 1000 && dist[23] < 1000)) {

	  read1 = 1;

	  read2 = 2;

	}

	else if ((dist[12] < 1000 && dist[34] < 1000)) {

	  read1 = 1;

	  read2 = 3;

	}

	else {

	  read1 = 0;

	}

	if (read1 != 0) {

	  if (mapped[read1] && mapped[read2]) {

	    count_norm++;

	    if (less_than(str[read1],chr[read1],pos[read1],str[read2],chr[read2],pos[read2])) {

	      #  print str[read1],chr[read1],pos[read1],str[read2],chr[read2],pos[read2],m[read1],cigarstr[read1],seq[read1],m[read2],cigarstr[read2],seq[read2],name[read1],name[read2],qual[read1],qual[read2] > fname1;

	      print str[read1],chr[read1],pos[read1],str[read2],chr[read2],pos[read2],m[read1],cigarstr[read1],seq[read1],m[read2],cigarstr[read2],seq[read2],name[read1],name[read2] > fname1;

	    }

	    else {

	      #  print str[read2],chr[read2],pos[read2],str[read1],chr[read1],pos[read1],m[read2],cigarstr[read2],seq[read2],m[read1],cigarstr[read1],seq[read1],name[read2],name[read1],qual[read2],qual[read1] > fname1;

	      print str[read2],chr[read2],pos[read2],str[read1],chr[read1],pos[read1],m[read2],cigarstr[read2],seq[read2],m[read1],cigarstr[read1],seq[read1],name[read2],name[read1] > fname1;

	    }

	  }

	  else {

	    for (i in c) {

	      print c[i] > fname3;

	    }	

	    count_unmapped++;

	  }

	}	

	else { 

	  # chimeric read with the 4 ends > 1KB apart

	  count_abnorm++;

	  for (i in c) {

	    print c[i] > fname2;

	  }

	}

      }

      else {

	dist[12] = abs(chr[1]-chr[2])*10000000 + abs(pos[1]-pos[2]);

	dist[23] = abs(chr[2]-chr[3])*10000000 + abs(pos[2]-pos[3]);

	dist[13] = abs(chr[1]-chr[3])*10000000 + abs(pos[1]-pos[3]);
@@ -450,9 +497,11 @@ END{ 
	  if (mapped[read1] && mapped[read2]) {

	    count_norm++;

	    if (less_than(str[read1],chr[read1],pos[read1],str[read2],chr[read2],pos[read2])) {

	      #  print str[read1],chr[read1],pos[read1],str[read2],chr[read2],pos[read2],m[read1],cigarstr[read1],seq[read1],m[read2],cigarstr[read2],seq[read2],name[read1],name[read2],qual[read1],qual[read2] > fname1;

	      print str[read1],chr[read1],pos[read1],str[read2],chr[read2],pos[read2],m[read1],cigarstr[read1],seq[read1],m[read2],cigarstr[read2],seq[read2],name[read1],name[read2] > fname1;

	    }

	    else {

	      #  print str[read2],chr[read2],pos[read2],str[read1],chr[read1],pos[read1],m[read2],cigarstr[read2],seq[read2],m[read1],cigarstr[read1],seq[read1],name[read2],name[read1],qual[read2],qual[read1] > fname1;

	      print str[read2],chr[read2],pos[read2],str[read1],chr[read1],pos[read1],m[read2],cigarstr[read2],seq[read2],m[read1],cigarstr[read1],seq[read1],name[read2],name[read1] > fname1;

	    }

	  }
@@ -471,6 +520,7 @@ END{ 
	  }

	}

      }

    }

    else if (count > 3) {

      # chimeric read > 3, too many to deal with

      count_abnorm++;
@@ -493,6 +543,9 @@ END{ 
        #qual[j] = tmp[11];

	name[j] = tmp[1];

	

        # blacklist - if 3rd bit set (=4) it means unmapped

        mapped[j] = and(tmp[2],4) == 0; 

	

	if (str[j] == 0 && tmp[6] ~/^[0-9]+S/) {

	  split(tmp[6], cigar, "S");

	  pos[j] = pos[j] - cigar[1];
@@ -532,17 +585,17 @@ END{ 
	    pos[j] = pos[j] - 16569;

	  }

	}

      mapped[j] = and(tmp[2],4) == 0;

	j++;

      }

      if (mapped[0] && mapped[1]) {

	count_reg++;

	if (less_than(str[0],chr[0],pos[0],str[1],chr[1],pos[1])) {

	  # ideally we'll get rid of printing out cigar string at some point

	# qual

	  #print str[0],chr[0],pos[0],str[1],chr[1],pos[1],m[0],cigarstr[0],seq[0],m[1],cigarstr[1],seq[1],name[0],name[1],qual[0],qual[1] > fname1;

	  print str[0],chr[0],pos[0],str[1],chr[1],pos[1],m[0],cigarstr[0],seq[0],m[1],cigarstr[1],seq[1],name[0],name[1] > fname1;

	}

	else {

	  #print str[1],chr[1],pos[1],str[0],chr[0],pos[0],m[1],cigarstr[1],seq[1],m[0],cigarstr[0],seq[0],name[1],name[0],qual[1],qual[0] > fname1;

print str[1],chr[1],pos[1],str[0],chr[0],pos[0],m[1],cigarstr[1],seq[1],m[0],cigarstr[0],seq[0],name[1],name[0] > fname1;

	}

      }
@@ -553,6 +606,13 @@ END{ 
	count_unmapped++;

      }

    }

    else if (count == 1) {

      # this actually shouldn't happen, but it happens with alternate aligners on occasion

      count_abnorm++;

      for (i in c) {

	print c[i] > fname2;

      }

    }

  resfile=fname1".res.txt";

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

}

CRA Git | Maintained and supported by SUSTech CRA and CCSE