Commit 74ae9e04 authored by nchernia's avatar nchernia
Browse files

single dedup

parent b19b7a77
Loading
Loading
Loading
Loading
+41 −6
Changes for CPU/common/chimeric_sam.awk: 41 added lines, 6 removed lines.
Original line number Diff line number Diff line
@@ -596,7 +596,7 @@ $0 !~ /^@/{
	      interiorpos1 = adjust(pos[read1],str[read1],cigarstr[read1],notprimary[read1]);
	      interiorpos2 = adjust(pos[read2],str[read2],cigarstr[read2],notprimary[read2]);
	      sortpos=sprintf("%0" chrlen "d_%0" chrlen "d",interiorpos1,interiorpos2);
	      cb_str = "cb:Z:"chr[read1]"_"chr[2]"_"fragstr"_"str[read1]"_"outputstr"_"sortpos;
	      cb_str = "cb:Z:"chr[read1]"_"chr[read2]"_"fragstr"_"str[read1]"_"outputstr"_"sortpos;

              # assign mate mapping quality, "read type", "interior position", "mate interior position"
	      # unclear for singleend what the read type should be; both 0/1 and 2/3 are possible
@@ -643,11 +643,29 @@ $0 !~ /^@/{
	split(c[i], tmp);
        # blacklist - if 3rd bit set (=4) it means unmapped
        mapped[j] = and(tmp[2],4) == 0;
        str[j] = and(tmp[2],16);
        # chromosome
        chr[j] = tmp[3];
        # get rid of "_" in chromosome name, for cb_str
        gsub(/_/, "", chr[j]);
        # position
        pos[j] = tmp[4];
        # cigar string
        cigarstr[j] = tmp[6];
      }
      if (mapped[0]) { 
      if (mapped[0] && length(singleend)>0) { 
	count_singleton++;
	# get 3' end for second end for dedupping
	interiorpos1 = adjust(pos[0],str[0],cigarstr[0],0);
	sortpos=sprintf("%0" chrlen "d_%0" chrlen "d",pos[0],interiorpos1);
	cb_str = "cb:Z:"chr[0]"_"chr[0]"_0_1_"str[0]"_"str[0]"_"sortpos;
	# read type 7 is singleton reads
	print c[1],"rt:A:7",cb_str;
      }
      else if (mapped[0]) {
	count_singleton++;
	for (i in c) {
	  print c[i],"rt:A:8";
	  print c[i],"rt:A:7";
	}
      }
      else {
@@ -1106,14 +1124,31 @@ END{
      j=0;
      for (i in c) {
	split(c[i], tmp);
	split(tmp[1],readname,"/");
        # blacklist - if 3rd bit set (=4) it means unmapped
        mapped[j] = and(tmp[2],4) == 0;
        str[j] = and(tmp[2],16);
        # chromosome
        chr[j] = tmp[3];
        # get rid of "_" in chromosome name, for cb_str
        gsub(/_/, "", chr[j]);
        # position
        pos[j] = tmp[4];
        # cigar string
        cigarstr[j] = tmp[6];
      }
      if (mapped[0]) { 
      if (mapped[0] && length(singleend)>0) { 
	count_singleton++;
	# get 3' end for second end for dedupping
	interiorpos1 = adjust(pos[0],str[0],cigarstr[0],0);
	sortpos=sprintf("%0" chrlen "d_%0" chrlen "d",pos[0],interiorpos1);
	cb_str = "cb:Z:"chr[0]"_"chr[0]"_0_1_"str[0]"_"str[0]"_"sortpos;
	# read type 7 is singleton reads
	print c[1],"rt:A:7",cb_str;
      }
      else if (mapped[0]) {
	count_singleton++;
	for (i in c) {
	  print c[i],"rt:A:8";
	  print c[i],"rt:A:7";
	}
      }
      else {