src/hg/tcga/bamBam/bamBam.c 1.1
1.1 2009/09/25 05:46:04 jsanborn
initial commit
Index: src/hg/tcga/bamBam/bamBam.c
===================================================================
RCS file: src/hg/tcga/bamBam/bamBam.c
diff -N src/hg/tcga/bamBam/bamBam.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ src/hg/tcga/bamBam/bamBam.c 25 Sep 2009 05:46:04 -0000 1.1
@@ -0,0 +1,253 @@
+/* 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"
+ "options:\n"
+ "\n"
+ );
+}
+
+static struct optionSpec options[] = {
+ {NULL, 0}
+};
+
+/* Option variables */
+char *position = NULL;
+
+typedef struct {
+ int beg, end;
+ samfile_t *in;
+} tmpstruct_t;
+
+typedef struct {
+ uint32_t format;
+ int tid, len, last_pos;
+ int mask;
+ char *ref;
+
+ /* Left pileup */
+ bam_header_t *leftH;
+ uint32_t leftTid;
+ uint32_t leftPos;
+ int leftN;
+ bam_pileup1_t *leftPu;
+
+ /* Right pileup */
+ bam_header_t *rightH;
+ uint32_t rightTid;
+ uint32_t rightPos;
+ int rightN;
+ bam_pileup1_t *rightPu;
+
+} pileup_data_t;
+
+static void reset_pileups(pileup_data_t *d)
+{
+d->leftPos = -1;
+d->leftTid = -1;
+d->leftN = -1;
+d->leftPu = NULL;
+
+d->rightPos = -1;
+d->rightTid = -1;
+d->rightN = -1;
+d->rightPu = NULL;
+}
+
+static int perform_analysis(pileup_data_t *d)
+{
+printf("%d\t%d\t%d\t", d->leftTid, d->leftPos, d->leftN);
+
+int i;
+const bam_pileup1_t *p;
+for (i = 0; i < d->leftN; ++i)
+ {
+ p = d->leftPu + i;
+ int c = toupper(bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)]);
+ putchar(c);
+ }
+putchar('\t');
+printf("%d\t", d->rightN);
+
+for (i = 0; i < d->rightN; ++i)
+ {
+ p = d->rightPu + 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)
+{ // will put everything to stdout.
+pileup_data_t *d = (pileup_data_t*)data;
+
+d->leftTid = tid;
+d->leftPos = pos;
+d->leftN = n;
+d->leftPu = (bam_pileup1_t*)pu;
+
+if (!d->rightPu)
+ return 1; // return '1' to stop reading data to let right gather some data
+else if (d->rightPos < d->leftPos)
+ return 1; // return '1' to stop reading data to let right catch up
+else if (d->rightPos > d->leftPos)
+ return 0; // return '0' to continue reading on left side data
+
+return perform_analysis(d);
+}
+
+static int pileup_func_right(uint32_t tid, uint32_t pos, int n,
+ const bam_pileup1_t *pu, void *data)
+{ // will put everything to stdout.
+pileup_data_t *d = (pileup_data_t*)data;
+
+d->rightTid = tid;
+d->rightPos = pos;
+d->rightN = n;
+d->rightPu = (bam_pileup1_t*)pu;
+
+if (!d->leftPu)
+ return 1; // return '1' to stop reading data to let left gather some data
+else if (d->leftPos < d->rightPos)
+ return 1; // return '1' to stop reading data to let left catch up
+else if (d->leftPos > d->rightPos)
+ return 0; // return '0' to continue reading on right side data
+
+return perform_analysis(d);
+}
+
+void bamBam(char *inbam1, char *inbam2)
+{
+tmpstruct_t left, right;
+left.beg = 0; left.end = 0x7fffffff;
+left.in = samopen(inbam1, "rb", NULL);
+if (left.in == NULL)
+ errAbort("samopen(%s, \"rb\") returned NULL", inbam1);
+
+right.beg = 0; right.end = 0x7fffffff;
+right.in = samopen(inbam2, "rb", NULL);
+if (right.in == NULL)
+ errAbort("samopen(%s, \"rb\") returned NULL", inbam2);
+
+pileup_data_t *d = (pileup_data_t*)calloc(1, sizeof(pileup_data_t));
+
+/* Init pileup_data object */
+d->tid = -1;
+d->mask = BAM_DEF_MASK;
+d->leftH = left.in->header;
+d->rightH = right.in->header;
+reset_pileups(d);
+
+bam1_t *bl = bam_init1();
+bam_plbuf_t *bufLeft = bam_plbuf_init(pileup_func_left, d);
+bam_plbuf_set_mask(bufLeft, -1); // mask = -1, no masking?
+
+bam1_t *br = bam_init1();
+bam_plbuf_t *bufRight = bam_plbuf_init(pileup_func_right, d);
+bam_plbuf_set_mask(bufRight, -1); // mask = -1, no masking?
+
+int retLeft = 1;
+int retRight = 1;
+int stopLeft = 0;
+int stopRight = 0;
+
+do
+ {
+ while (!stopLeft && ((retLeft = samread(left.in, bl)) >= 0))
+ stopLeft = bam_dual_plbuf_push(bl, bufLeft);
+
+ while (!stopRight && ((retRight = samread(right.in, br)) >= 0))
+ stopRight = bam_dual_plbuf_push(br, bufRight);
+
+ while (stopLeft && stopRight)
+ {
+ if (!d->leftPu && !d->rightPu)
+ {
+ stopLeft = bam_dual_plbuf_resume(bufLeft);
+ stopRight = bam_dual_plbuf_resume(bufRight);
+ }
+ else if (!d->leftPu)
+ stopLeft = bam_dual_plbuf_resume(bufLeft);
+ else if (!d->rightPu)
+ stopRight = bam_dual_plbuf_resume(bufRight);
+ else if (d->leftPos > 0 && d->leftPos < d->rightPos)
+ stopLeft = bam_dual_plbuf_resume(bufLeft);
+ else if (d->rightPos > 0 && d->rightPos < d->leftPos)
+ stopRight = bam_dual_plbuf_resume(bufRight);
+ }
+
+ } while (retLeft >= 0 && retRight >= 0);
+
+if (retLeft < 0)
+ {
+ bam_dual_plbuf_push(0, bufLeft); // finalize pileup
+ bam_dual_plbuf_destroy(bufLeft);
+ }
+
+if (retRight < 0)
+ {
+ bam_dual_plbuf_push(0, bufRight); // finalize pileup
+ bam_dual_plbuf_destroy(bufRight);
+ }
+
+bam_destroy1(bl);
+bam_destroy1(br);
+
+samclose(left.in);
+samclose(right.in);
+}
+
+int main(int argc, char *argv[])
+/* Process command line. */
+{
+optionInit(&argc, argv, options);
+if (argc != 3)
+ usage();
+
+char *inbam1 = argv[1];
+char *inbam2 = argv[2];
+bamBam(inbam1, inbam2);
+
+return 0;
+}
+
+#else // USE_BAM not defined
+
+int main(int argc, char *argv[])
+/* Process command line. */
+{
+printf("BAM not installed");
+
+return 0;
+}
+#endif