Index: src/hg/tcga/bamBam/bamBam.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/tcga/bamBam/bamBam.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -b -B -U 1000000 -r1.3 -r1.4
--- src/hg/tcga/bamBam/bamBam.c	25 Sep 2009 18:44:27 -0000	1.3
+++ src/hg/tcga/bamBam/bamBam.c	25 Sep 2009 21:16:13 -0000	1.4
@@ -1,263 +1,108 @@
 /* bamBam -- discovering variants from BAM file
  *   code adapted from Heng Li's Samtools library functions
  */
 
 #ifdef USE_BAM
 
 #include "common.h"
 #include "hCommon.h"
 #include "linefile.h"
 #include "hash.h"
 #include "hdb.h"
 #include "options.h"
 #include "jksql.h"
 
 #define _IOLIB 2
 #include "bam.h"
 #include "sam.h"
 #include "bam_dual_pileup.h"
 
 static char const rcsid[] = "$Id";
 
 void usage()
 /* Explain usage and exit. */
 {
 errAbort(
   "bamBam - Perform analysis on two BAM files.\n"
   "usage:\n"
-  "   bamBam [options] <in.bam1> <in.bam2>\n"
+  "   bamBam [options] <left.bam> <right.bam>\n"
   "options:\n"
   "\n"
   );
 }
 
 static struct optionSpec options[] = {
     {NULL, 0}
 };
 
 /* Option variables */
 char *position = NULL;
 
 typedef struct {
-    int beg, end;
-    samfile_t *in;
-} tmpstruct_t;
+    int minQual;
+    int minMapQual;
+} analysis_params;
 
-typedef struct {
-    uint32_t format;
-    int tid, len, last_pos;
-    int mask;
-    char *ref;
-
-    /* Left pileup */
-    bam_header_t *hL;
-    int32_t tidL;
-    int32_t posL;
-    int nL;
-    bam_pileup1_t *puL;
-
-    /* Right pileup */
-    bam_header_t *hR;
-    int32_t tidR;
-    int32_t posR;
-    int nR;
-    bam_pileup1_t *puR;
-
-} pileup_data_t;
-
-static void reset_pileups(pileup_data_t *d)
+static void perform_analysis(void *data)
 {
-d->posL = -1;
-d->tidL = -1;
-d->nL   = -1;
-d->puL  = NULL;
-
-d->posR = -1;
-d->tidR = -1;
-d->nR   = -1;
-d->puR  = NULL;
-}
+pileup_data_t *d = (pileup_data_t*)data;
+analysis_params *ap = (analysis_params*)d->func_data;
 
-static int perform_analysis(pileup_data_t *d)
-{
-printf("%d\t%d\t%d\t", d->tidL, d->posL, d->nL);
+printf("%d\t%d\t%d\t%d/%d\t", d->tidL, d->posL, d->nL, ap->minQual, ap->minMapQual);
 
 int i;
 const bam_pileup1_t *p;
 for (i = 0; i < d->nL; ++i)
     {
     p = d->puL + i;
     int c = toupper(bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)]);
     putchar(c);
     }
 putchar('\t');
 printf("%d\t", d->nR);
 
 for (i = 0; i < d->nR; ++i)
     {
     p = d->puR + i;
     int c = toupper(bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)]);
     putchar(c);
     }
 putchar('\n');
-    
-reset_pileups(d);
-
-return 0;
 }
 
-static int pileup_func_left(uint32_t tid, uint32_t pos, int n, 
-			    const bam_pileup1_t *pu, void *data)
-{
-pileup_data_t *d = (pileup_data_t*)data;
-
-d->tidL = tid;
-d->posL = pos;
-d->nL   = n;
-d->puL  = (bam_pileup1_t*)pu;
-
-if (!d->puR) 
-    return 1;  // return '1' to stop reading to let right get some data
-else if (d->tidR < d->tidL)
-    return 1;  // return '1' to stop reading to let right catch up to chrom
-else if (d->tidR > d->tidL)
-    return 0;  // return '0' to catch up to right's chrom 
-else if (d->posR < d->posL)
-    return 1;  // return '1' to stop reading data to let right catch up
-else if (d->posR > d->posL)
-    return 0;  // return '0' to catch up to left's pos
-
-// only perform analysis when L & R positions line up.
-return perform_analysis(d);
-}
-
-static int pileup_func_right(uint32_t tid, uint32_t pos, int n, 
-			     const bam_pileup1_t *pu, void *data)
-{
-pileup_data_t *d = (pileup_data_t*)data;
-
-d->tidR = tid;
-d->posR = pos;
-d->nR   = n;
-d->puR  = (bam_pileup1_t*)pu;
-
-if (!d->puL) 
-    return 1;  // return '1' to stop reading data to let left get data
-else if (d->tidL < d->tidR)
-    return 1;  // return '1' to stop reading data to let left catch up to chrom
-else if (d->tidR > d->tidR)
-    return 0; // return '0' to catch up to left's chrom.
-else if (d->posL < d->posR)
-    return 1;  // return '1' to stop reading data to let left catch up
-else if (d->posL > d->posR)
-    return 0;  // return '0' to catch up to left's pos
-
-// only perform analysis when L & R positions line up.
-return perform_analysis(d);
-}
 
 void bamBam(char *inbamL, char *inbamR)
 {
-tmpstruct_t left, right;  
-left.beg = 0; left.end = 0x7fffffff;
-left.in = samopen(inbamL, "rb", NULL);
-if (left.in == NULL)
-    errAbort("samopen(%s, \"rb\") returned NULL", inbamL);
-
-right.beg = 0; right.end = 0x7fffffff;
-right.in = samopen(inbamR, "rb", NULL);
-if (right.in == NULL)
-    errAbort("samopen(%s, \"rb\") returned NULL", inbamR);
-
-// Init pileup_data object
-pileup_data_t *d = (pileup_data_t*)calloc(1, sizeof(pileup_data_t));
-d->tid = -1; 
-d->mask = BAM_DEF_MASK;
-d->hL = left.in->header;
-d->hR = right.in->header;
-reset_pileups(d);
-
-// Initialize read & pileup buffers for left side
-bam1_t *bL = bam_init1();
-bam_plbuf_t *bufL = bam_plbuf_init(pileup_func_left, d);
-bam_plbuf_set_mask(bufL, -1);  // mask = -1, no masking?
-
-// Initialize read & pileup buffers for right side
-bam1_t *bR = bam_init1();
-bam_plbuf_t *bufR = bam_plbuf_init(pileup_func_right, d);
-bam_plbuf_set_mask(bufR, -1);  // mask = -1, no masking?
-
-int retL = 1; 
-int retR = 1;
-int pauseL = 0;
-int pauseR = 0;
-
-do 
-    {
-    while (!pauseL && ((retL = samread(left.in, bL)) >= 0))
-        pauseL = bam_dual_plbuf_push(bL, bufL);
-    
-    while (!pauseR && ((retR = samread(right.in, bR)) >= 0))
-	pauseR = bam_dual_plbuf_push(bR, bufR);
-
-    while (pauseL && pauseR)
-	{  // only allow one change per round.
-	if (!d->puL)                // L needs data
-	    pauseL = bam_dual_plbuf_resume(bufL);
-	else if (!d->puR)           // R needs data
-	    pauseR = bam_dual_plbuf_resume(bufR);
-	else if (d->tidL < d->tidR) // L needs to catch up to R's chrom
-	    pauseL = bam_dual_plbuf_resume(bufL);
-	else if (d->tidL > d->tidR) // R needs to catch up to L's chrom
-	    pauseR = bam_dual_plbuf_resume(bufR);
-	else if (d->posL > 0 && d->posL < d->posR) // L needs to catch up to R's pos
-	    pauseL = bam_dual_plbuf_resume(bufL);
-	else if (d->posR > 0 && d->posR < d->posL) // R needs to catch up to L's pos
-	    pauseR = bam_dual_plbuf_resume(bufR);
-	}
 
-    } while (retL >= 0 && retR >= 0);
+analysis_params aps;
+aps.minQual = 50;
+aps.minMapQual = 25;
 
-if (retL < 0)
-    {
-    bam_dual_plbuf_push(0, bufL);  // finalize pileup
-    bam_dual_plbuf_destroy(bufL);
-    }
-
-if (retR < 0)
-    {
-    bam_dual_plbuf_push(0, bufR);  // finalize pileup
-    bam_dual_plbuf_destroy(bufR);
-    }
-
-bam_destroy1(bL);
-bam_destroy1(bR);
-
-samclose(left.in);
-samclose(right.in);
+bam_bam_file(inbamL, inbamR, perform_analysis, &aps);
 }
 
+
 int main(int argc, char *argv[])
 /* Process command line. */
 {
 optionInit(&argc, argv, options);
 if (argc != 3)
     usage();
 
 char *inbamL = argv[1];
 char *inbamR = argv[2];
 bamBam(inbamL, inbamR);
 
 return 0;
 }
 
 #else // USE_BAM not defined
 
 int main(int argc, char *argv[])
 /* Process command line. */
 {
 printf("BAM not installed");
 
 return 0;
 }
 #endif

