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

1.7 2010/03/08 19:56:11 jsanborn
updated
Index: src/hg/tcga/pBamBam/asBamBam.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/tcga/pBamBam/asBamBam.c,v
retrieving revision 1.6
retrieving revision 1.7
diff -b -B -U 4 -r1.6 -r1.7
--- src/hg/tcga/pBamBam/asBamBam.c	7 Mar 2010 05:49:45 -0000	1.6
+++ src/hg/tcga/pBamBam/asBamBam.c	8 Mar 2010 19:56:11 -0000	1.7
@@ -27,9 +27,10 @@
   "usage:\n"
   "   asBamBam [options] <left.bam> <right.bam>\n"
   "options:\n"
   "   -position=chr1:1-100 - Position to do analysis. [NULL=process all reads]\n" 
-  "   -minSupport=INT   - Minimum Support in Germline. [4]\n"
+  "   -minSuppSNP=INT   - Minimum Support for SNP/Het Calling. [4]\n"
+  "   -minSuppSV=INT   - Minimum Support for struct vars. [2]\n"
   "   -minQ=INT         - Minimum acceptable base quality. [10]\n"
   "   -minMapQ=INT      - Minimum acceptable mapping quality. [20]\n"
   "   -avgMapQ=INT      - Minimum acceptable average mapping quality. [30]\n"
   "   -minChiSq=INT     - Minimum acceptable chi square of nucleotide freqs. [250]\n"
@@ -40,9 +41,10 @@
 }
 
 static struct optionSpec options[] = {
     {"position", OPTION_STRING},
-    {"minSupport", OPTION_INT},
+    {"minSuppSNP", OPTION_INT},
+    {"minSuppSV", OPTION_INT},
     {"minQ", OPTION_INT},
     {"minMapQ", OPTION_INT},
     {"avgMapQ", OPTION_INT},
     {"maxISize", OPTION_INT},
@@ -54,9 +56,10 @@
 #define BAD_TID -999
 
 /* Option variables */
 char *position  = NULL;
-int minSupport  = 4;
+int minSuppSNP  = 4;
+int minSuppSV   = 2;
 int minQ        = 10;
 int minMapQ     = 20;
 int avgMapQ     = 30;
 int maxISize    = 2000;
@@ -68,9 +71,10 @@
     int minMapQ;
     int avgMapQ;
     int maxISize;
     int minChiSq;
-    int minSupp;
+    int minSuppSNP;
+    int minSuppSV;
     int sites;
 
     int readsPerBin;
 
@@ -371,9 +375,9 @@
 
 majorL.nuc   = 'N';
 majorL.vote  = 0;
 majorL.count = 0;
-if (iL1 != -1 && countL1 >= ap->minSupp)
+if (iL1 != -1 && countL1 >= ap->minSuppSNP)
     {
     majorL.nuc   = NUCS[iL1];
     majorL.vote  = maxL1;
     majorL.count = countL1;
@@ -381,9 +385,9 @@
 
 minorL.nuc   = 'N';
 minorL.vote  = 0;
 minorL.count = 0;
-if (iL2 != -1 && countL2 >= ap->minSupp)
+if (iL2 != -1 && countL2 >= ap->minSuppSNP)
     {
     minorL.nuc   = NUCS[iL2];
     minorL.vote  = maxL2;
     minorL.count = countL2;
@@ -391,9 +395,9 @@
 
 majorR.nuc   = 'N';
 majorR.vote  = 0;
 majorR.count = 0;
-if (iR1 != -1 && countR1 >= ap->minSupp)
+if (iR1 != -1 && countR1 >= ap->minSuppSNP)
     {
     majorR.nuc   = NUCS[iR1];
     majorR.vote  = maxR1;
     majorR.count = countR1;
@@ -401,9 +405,9 @@
 
 minorR.nuc   = 'N';
 minorR.vote  = 0;
 minorR.count = 0;
-if (iR2 != -1 && countR2 >= ap->minSupp)
+if (iR2 != -1 && countR2 >= ap->minSuppSNP)
     {
     minorR.nuc   = NUCS[iR2];
     minorR.vote  = maxR2;
     minorR.count = countR2;
@@ -681,9 +685,9 @@
 int i = 0;
 int rCount = 0, maxCountR = 0;
 for (i = 0; i < 4; i++)
     {
-    if (countsR[i] < ap->minSupp)
+    if (countsR[i] < ap->minSuppSNP)
 	continue;
     rCount += countsR[i];
     if (countsR[i] > maxCountR)
 	maxCountR = countsR[i];
@@ -920,9 +924,9 @@
 	return 0;  // there still might be reads to overlap
 
     /* check average mapping quality of discordants */
     int qual = ceil((double) interL.qual / (double) interL.supp); 
-    if (interL.supp >= ap->minSupp && qual >= ap->avgMapQ)
+    if (interL.supp >= ap->minSuppSV && qual >= ap->avgMapQ)
 	return 1;
     reset_inter();
     }
 
@@ -1088,9 +1092,9 @@
 	return 0;  // there still might be reads to overlap
 
     /* check average mapping quality of discordants */
     int qual = ceil((double) intraL.qual / (double) intraL.supp); 
-    if (intraL.supp >= ap->minSupp && qual >= ap->avgMapQ)
+    if (intraL.supp >= ap->minSuppSV && qual >= ap->avgMapQ)
 	return 1;
     reset_intra();
     }
 
@@ -1099,9 +1103,9 @@
 
 static int indel_variant(dual_pileup_t *d, analysis_params_t *ap,
 			 int *isIns, int *isDel, int *insGerm, int *delGerm)
 {
-if (insL >= ap->minSupp)
+if (insL >= ap->minSuppSNP)
     {
     *isIns = 1;
     if (insR > 0)  // if you see even once, believe it.
 	*insGerm = 1;
@@ -1109,9 +1113,9 @@
 	*insGerm = 0;
     return 1;
     }
 
-if (delL >= ap->minSupp)
+if (delL >= ap->minSuppSNP)
     {
     *isDel = 1;
     if (delR > 0)  // if you see even once, believe it.
 	*delGerm = 1;
@@ -1264,9 +1268,9 @@
 	minSa = NUCS[i];
 	}
     }
 
-if (minSc < ap->minSupp)
+if (minSc < ap->minSuppSNP)
     {
     minSc = 0;
     minSa = 'N';
     }
@@ -1379,9 +1383,10 @@
 ap.minMapQ  = minMapQ;
 ap.avgMapQ  = avgMapQ;
 ap.maxISize = maxISize;
 ap.minChiSq = minChiSq;
-ap.minSupp  = minSupport;
+ap.minSuppSNP  = minSuppSNP;
+ap.minSuppSV   = minSuppSV;
 ap.sites    = 0;
 
 ap.readsPerBin = readsPerBin;
 
@@ -1433,10 +1438,12 @@
     usage();
 
 if (optionExists("position"))
     position = optionVal("position", NULL);
-if (optionExists("minSupport"))
-    minSupport = optionInt("minSupport", 4);
+if (optionExists("minSuppSNP"))
+    minSuppSNP = optionInt("minSuppSNP", 4);
+if (optionExists("minSuppSV"))
+    minSuppSV = optionInt("minSuppSV", 2);
 if (optionExists("minQ"))
     minQ = optionInt("minQ", 10);
 if (optionExists("minMapQ"))
     minMapQ = optionInt("minMapQ", 20);