src/hg/tcga/pBamBam/pBamBamMuts.c 1.2

Index: src/hg/tcga/pBamBam/pBamBamMuts.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/tcga/pBamBam/pBamBamMuts.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -b -B -U 4 -r1.1 -r1.2
--- src/hg/tcga/pBamBam/pBamBamMuts.c	23 Jan 2010 04:24:49 -0000	1.1
+++ src/hg/tcga/pBamBam/pBamBamMuts.c	4 Mar 2010 04:03:51 -0000	1.2
@@ -328,8 +328,10 @@
 static int model_select(dual_pileup_t *d, analysis_params_t *ap, 
 			char *callL1, char *callL2, char *callR1, char *callR2, 
 			float *hetL, float *hetR, int *stateL, int *stateR, int *sites)
 {
+*sites += 1;
+
 if (majorR.vote == 0 || majorL.vote == 0)
     return 0;  // most likely a bad quality region.
 
 if (sig_diff(d, ap) < ap->minChiSq)
@@ -342,48 +344,76 @@
 
 if (rCount < 10)
     return 0;
 
+int retVal1 = 0, retVal2 = 0;
 for (i = 0; i < 4; i++)
     {
     if (NUCS[i] != majorL.nuc && NUCS[i] != minorL.nuc)
 	continue;
     
     if (NUCS[i] == majorL.nuc && majorL.nuc != majorR.nuc && majorL.nuc != minorR.nuc)
 	{
-	if (countsL1[i] < 2 || countsL2[i] < 1)
-	    return 0;
-	if (countsL1[i] == 0 || countsL2[i] == 0 || countsL3[i] == 0)
-	    return 0;
-	if (rCount >= 30 && countsR[i] > 1) 
-	    return 0;
-	if (rCount < 30 && countsR[i] > 0)
-	    return 0;
+//	if (countsL1[i] >= 2 && countsL2[i] >= 1)
+//	    retVal1 = 1;
+//	else if (countsL1[i] > 0 && countsL2[i] > 0 && countsL3[i] > 0)
+//	    retVal1 = 1;
+
+	if (countsL1[i] > 0 && countsL2[i] > 0 && countsL1[i] + countsL2[i] >= 3)
+	    retVal1 = 1;
+	else if (countsL1[i] > 0 && countsL3[i] > 0 && countsL1[i] + countsL3[i] >= 3)
+	    retVal1 = 1;
+	else if (countsL2[i] > 0 && countsL3[i] > 0 && countsL2[i] + countsL3[i] >= 3)
+	    retVal1 = 1;
+
+	if (retVal1)
+	    {
+	    if (rCount >= 30 && countsR[i] <= 1) 
+		retVal1 = 1;
+	    else if (rCount < 30 && countsR[i] == 0)
+		retVal1 = 1;
+	    else
+		retVal1 = 0;
+	    }
 	}
 
 
     if (NUCS[i] == minorL.nuc && (minorL.nuc != majorR.nuc && minorL.nuc != minorR.nuc))
 	{
-	if (countsL1[i] < 2 || countsL2[i] < 1)
-	    return 0;
-	if (countsL1[i] == 0 || countsL2[i] == 0 || countsL3[i] == 0)
-	    return 0;
-	if (rCount >= 30 && countsR[i] > 1) 
-	    return 0;
-	if (rCount < 30 && countsR[i] > 0)
-	    return 0;
+//	if (countsL1[i] >= 2 && countsL2[i] >= 1)
+//	    retVal2 = 1;
+//	else if (countsL1[i] > 0 && countsL2[i] > 0 && countsL3[i] > 0)
+//	    retVal2 = 1;
+
+	if (countsL1[i] > 0 && countsL2[i] > 0 && countsL1[i] + countsL2[i] >= 3)
+	    retVal2 = 1;
+	else if (countsL1[i] > 0 && countsL3[i] > 0 && countsL1[i] + countsL3[i] >= 3)
+	    retVal2 = 1;
+	else if (countsL2[i] > 0 && countsL3[i] > 0 && countsL2[i] + countsL3[i] >= 3)
+	    retVal2 = 1;
+
+	if (retVal2) 
+	    {
+	    if (rCount >= 30 && countsR[i] <= 1) 
+		retVal2 = 1;
+	    else if (rCount < 30 && countsR[i] == 0)
+		retVal2 = 1;
+	    else
+		retVal2 = 0;
+	    }
 	}
     }    
     
+if (!retVal1 && !retVal2)
+    return 0;
+
 	    
 float freqL = ((float) minorL.count) / ((float) majorL.count + (float) minorL.count);
 float freqR = ((float) minorR.count) / ((float) majorR.count + (float) minorR.count);
 
 *hetL = freqL;
 *hetR = freqR;
 
-*sites += 1;
-
 *callL1 = majorL.nuc;
 *callL2 = minorL.nuc;
 *callR1 = majorR.nuc;
 *callR2 = minorR.nuc;