Changed chimeric blacklist to handle quadruple reads and eliminate MT exception (a55bfe32) · Commits · Chaos / chaos.juicer

CPU/common/chimeric_blacklist.awk

+51 −6

Changes for CPU/common/chimeric_blacklist.awk: 51 added lines, 6 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];
		@@ -147,15 +147,59 @@ BEGIN{
		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]);
		@@ -207,6 +251,7 @@ BEGIN{
		}
		}
		}
		}
		else if (count > 3) {
		# chimeric read > 3, too many to deal with
		count_abnorm++;

Admin message