src/hg/tcga/pBamBam/asBamBam.c 1.9

Index: src/hg/tcga/pBamBam/asBamBam.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/tcga/pBamBam/asBamBam.c,v
retrieving revision 1.8
retrieving revision 1.9
diff -b -B -U 4 -r1.8 -r1.9
--- src/hg/tcga/pBamBam/asBamBam.c	14 Mar 2010 23:58:25 -0000	1.8
+++ src/hg/tcga/pBamBam/asBamBam.c	15 Mar 2010 00:24:27 -0000	1.9
@@ -230,10 +230,18 @@
 reset_globals(ap);
 
 int i, index;
 int32_t start, stop;
-int32_t prevStart = -1, prevStop = -1;
+int32_t prevStart = -1;
 const bam_pileup1_t *p;
+
+char best_n1  = 'N';
+int best_q1 = -1;
+char best_n2  = 'N';
+int best_q2 = -1;
+
+int insertion = 0, deletion = 0;
+ 
 for (i = 0; i < d->nL; ++i)
     { 
     p = d->puL + i;
     if (!check_read(p, ap))
@@ -243,36 +251,76 @@
 
     /* Check if start & stop has already been visited, to dupes missed in prev step*/    
     start = (c)->pos;
     stop  = start + (c)->l_qseq;    
-    if (start <= prevStart && stop <= prevStop)
-	continue;
-    prevStart = start;
-    prevStop  = stop;
 
+    //if (start <= prevStart && stop <= prevStop)
+//	continue;
     if (ap->pStart == -1)
 	{
 	ap->pStart = start;
 	ap->pStop  = stop;
 	}
+    else 
+	{
     if (start < ap->pStart)
 	ap->pStart = start;
     if (stop  > ap->pStop)
 	ap->pStop = stop;
+	}
+
+    if (prevStart == -1)
+	prevStart = start;
+
+    if (start < prevStart)
+	continue;
 
+    if (start == prevStart)
+	{
+	if (p->is_del)
+	    {
     if (p->indel > 0)
-	insL += 1;
+		insertion = 1;
     else if (p->indel < 0)
-	delL += 1;
-    
-    if (!p->is_del && p->indel == 0)
+		deletion = 1;
+	    }
+	else
 	{
-        double lseq = (double) p->b->core.l_qseq;
 	int b = toupper(bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)]);    
 	if ((index = nuc_index(b)) < 0)
 	    continue;
+	    if (best_b1 == 'N')
+		{
+		best_b1 = b;
+		best_q1 = qual;
+		}
+	    else if (qual > best_q1)
+		{
+		best_b2 = best_b1;
+		best_q2 = best_q1;
 
-	votesL[index] += qual;
+		best_b1 = b;
+		best_q1 = qual;
+		}
+	    else if (qual > best_q2)
+		{
+		best_b2 = b;
+		best_q2 = qual;
+		}
+	    }
+	continue;
+	}
+    prevStart = start;
+    
+    if (insertion)
+	insL += 1;
+    else if (deletion)
+	delL += 1;
+
+    double lseq = (double) p->b->core.l_qseq;
+    if (!(index = nuc_index(best_b1)) < 0)
+	{
+	votesL[index] += best_q1;
 	countsL[index]++;
 
         if (p->qpos < round(0.33 * lseq))
 	    countsL1[index]++;
@@ -280,12 +328,30 @@
 	    countsL2[index]++;
 	else
 	    countsL3[index]++;
 	}
+    
+    if (!(index = nuc_index(best_b2)) < 0)
+	{
+	votesL[index] += best_q2;
+	countsL[index]++;
+	
+        if (p->qpos < round(0.33 * lseq))
+	    countsL1[index]++;
+	else if (p->qpos < round(0.67 * lseq))
+	    countsL2[index]++;
+	else
+	    countsL3[index]++;
+	}
+    best_b1 = 'N';
+    best_q1 = -1;
+    best_b2 = 'N';
+    best_q2 = -1;
+    
+    insertion = 0, deletion = 0;
     }
 
 prevStart = -1;
-prevStop  = -1;
 for (i = 0; i < d->nR; ++i)
     {
     p = d->puR + i;
     if (!check_read(p, ap))