Commit c9de8924 authored by nchernia's avatar nchernia
Browse files

Now properly handle quadruple reads. Also eliminated MT check, which was not doing the right thing

parent 66e9807c

AWS/scripts/chimeric_blacklist.awk

+130 −27
Changes for AWS/scripts/chimeric_blacklist.awk: 130 added lines, 27 removed lines.
Original line number Diff line number Diff line
@@ -86,9 +86,9 @@ BEGIN{ 
  else {

    # 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];
@@ -106,7 +106,7 @@ BEGIN{ 
	cigarstr[j] = tmp[6];

	# sequence

	seq[j] = tmp[10];

        qual[j] = tmp[11];

        #qual[j] = tmp[11];

        

	# get rid of soft clipping to know correct position

	if (str[j] == 0 && tmp[6] ~/^[0-9]+S/) {
@@ -146,16 +146,60 @@ BEGIN{ 
	    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]);
@@ -184,9 +228,11 @@ BEGIN{ 
	  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;

	    }

	  }
@@ -205,6 +251,7 @@ BEGIN{ 
	  }

	}

      }

    }

    else if (count > 3) {

      # chimeric read > 3, too many to deal with

      count_abnorm++;
@@ -224,7 +271,7 @@ BEGIN{ 
	m[j] = tmp[5];

	cigarstr[j] = tmp[6];

	seq[j] = tmp[10];

        qual[j] = tmp[11];

        #qual[j] = tmp[11];

	name[j] = tmp[1];

	

        # blacklist - if 3rd bit set (=4) it means unmapped
@@ -265,9 +312,6 @@ BEGIN{ 
	    cigloc = substr(tmp[6],where,RLENGTH-1) + 0;

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

	  }

					if (chr[j] ~ /MT/ && pos[j] >= 16569) {

						pos[j] = pos[j] - 16569;

					}

	}

	j++;

      }
@@ -275,9 +319,11 @@ BEGIN{ 
	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

	  #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;

	}

      }
@@ -306,15 +352,15 @@ BEGIN{ 
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];
@@ -326,7 +372,8 @@ END{ 
	cigarstr[j] = tmp[6];

	# sequence

	seq[j] = tmp[10];

        qual[j] = tmp[11];

        #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");
@@ -335,6 +382,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];
@@ -360,19 +408,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]);
@@ -401,9 +494,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;

	    }

	  }
@@ -422,6 +517,7 @@ END{ 
	  }

	}

      }

    }

    else if (count > 3) {

      # chimeric read > 3, too many to deal with

      count_abnorm++;
@@ -441,9 +537,11 @@ END{ 
	m[j] = tmp[5];

	cigarstr[j] = tmp[6];

	seq[j] = tmp[10];

        qual[j] = tmp[11];

        #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");
@@ -480,20 +578,18 @@ END{ 
	    cigloc = substr(tmp[6],where,RLENGTH-1) + 0;

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

	  }

					if (chr[j] ~ /MT/ && pos[j] >= 16569) {

						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

	  #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;

	}

      }
@@ -504,6 +600,13 @@ END{ 
	count_unmapped++;

      }

    }

    

  printf("%d %d %d %d %d\n", tottot, count_unmapped, count_reg, count_norm, count_abnorm) >> fname1".res.txt";

    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;

}

CPU/common/chimeric_blacklist.awk

+0 −6
Changes for CPU/common/chimeric_blacklist.awk: 0 added lines, 6 removed lines.
Original line number Diff line number Diff line
@@ -312,9 +312,6 @@ BEGIN{ 
	    cigloc = substr(tmp[6],where,RLENGTH-1) + 0;

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

	  }

	  if (chr[j] ~ /MT/ && pos[j] >= 16569) {

	    pos[j] = pos[j] - 16569;

	  }

	}

	j++;

      }
@@ -581,9 +578,6 @@ END{ 
	    cigloc = substr(tmp[6],where,RLENGTH-1) + 0;

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

	  }

	  if (chr[j] ~ /MT/ && pos[j] >= 16569) {

	    pos[j] = pos[j] - 16569;

	  }

	}

	j++;

      }

LSF/scripts/chimeric_blacklist.awk

+130 −27
Changes for LSF/scripts/chimeric_blacklist.awk: 130 added lines, 27 removed lines.
Original line number Diff line number Diff line
@@ -86,9 +86,9 @@ BEGIN{ 
  else {

    # 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];
@@ -106,7 +106,7 @@ BEGIN{ 
	cigarstr[j] = tmp[6];

	# sequence

	seq[j] = tmp[10];

        qual[j] = tmp[11];

        #qual[j] = tmp[11];

        

	# get rid of soft clipping to know correct position

	if (str[j] == 0 && tmp[6] ~/^[0-9]+S/) {
@@ -146,16 +146,60 @@ BEGIN{ 
	    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]);
@@ -184,9 +228,11 @@ BEGIN{ 
	  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;

	    }

	  }
@@ -205,6 +251,7 @@ BEGIN{ 
	  }

	}

      }

    }

    else if (count > 3) {

      # chimeric read > 3, too many to deal with

      count_abnorm++;
@@ -224,7 +271,7 @@ BEGIN{ 
	m[j] = tmp[5];

	cigarstr[j] = tmp[6];

	seq[j] = tmp[10];

        qual[j] = tmp[11];

        #qual[j] = tmp[11];

	name[j] = tmp[1];

	

        # blacklist - if 3rd bit set (=4) it means unmapped
@@ -265,9 +312,6 @@ BEGIN{ 
	    cigloc = substr(tmp[6],where,RLENGTH-1) + 0;

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

	  }

					if (chr[j] ~ /MT/ && pos[j] >= 16569) {

						pos[j] = pos[j] - 16569;

					}

	}

	j++;

      }
@@ -275,9 +319,11 @@ BEGIN{ 
	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

	  #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;

	}

      }
@@ -306,15 +352,15 @@ BEGIN{ 
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];
@@ -326,7 +372,8 @@ END{ 
	cigarstr[j] = tmp[6];

	# sequence

	seq[j] = tmp[10];

        qual[j] = tmp[11];

        #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");
@@ -335,6 +382,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];
@@ -360,19 +408,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]);
@@ -401,9 +494,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;

	    }

	  }
@@ -422,6 +517,7 @@ END{ 
	  }

	}

      }

    }

    else if (count > 3) {

      # chimeric read > 3, too many to deal with

      count_abnorm++;
@@ -441,9 +537,11 @@ END{ 
	m[j] = tmp[5];

	cigarstr[j] = tmp[6];

	seq[j] = tmp[10];

        qual[j] = tmp[11];

        #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");
@@ -480,20 +578,18 @@ END{ 
	    cigloc = substr(tmp[6],where,RLENGTH-1) + 0;

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

	  }

					if (chr[j] ~ /MT/ && pos[j] >= 16569) {

						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

	  #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;

	}

      }
@@ -504,6 +600,13 @@ END{ 
	count_unmapped++;

      }

    }

    

  printf("%d %d %d %d %d\n", tottot, count_unmapped, count_reg, count_norm, count_abnorm) >> fname1".res.txt";

    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;

}

PBS/scripts/chimeric_blacklist.awk

+130 −27

File changed.

Preview size limit exceeded, changes collapsed.

SLURM/scripts/chimeric_blacklist.awk

+131 −28

File changed.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE