src/hg/tcga/bamBam/bamBam.c 1.2
1.2 2009/09/25 18:08:57 jsanborn
it is alive
Index: src/hg/tcga/bamBam/bamBam.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/tcga/bamBam/bamBam.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -b -B -U 1000000 -r1.1 -r1.2
--- src/hg/tcga/bamBam/bamBam.c 25 Sep 2009 05:46:04 -0000 1.1
+++ src/hg/tcga/bamBam/bamBam.c 25 Sep 2009 18:08:57 -0000 1.2
@@ -1,253 +1,295 @@
/* 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;
+ int32_t leftTid;
+ int32_t leftPos;
int leftN;
bam_pileup1_t *leftPu;
/* Right pileup */
bam_header_t *rightH;
- uint32_t rightTid;
- uint32_t rightPos;
+ int32_t rightTid;
+ int32_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)
{
+//return 0;
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->rightTid < d->leftTid)
+ return 1;
+else if (d->rightTid > d->leftTid)
+ return 0;
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->leftTid < d->rightTid)
+ return 1;
+else if (d->leftTid > d->rightTid)
+ return 0;
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;
-
+int count = 0;
+int debug = 0;
do
{
+
+ if (d->leftPos >= 247194573 || d->rightPos >= 247194573)
+ debug = 0;
+
+ if (debug)
+ {
+ count++;
+ printf("tid: %d, %d\n", d->leftTid, d->rightTid);
+ printf("pos: %d, %d\n", d->leftPos, d->rightPos);
+ printf("stop: %d, %d\n", stopLeft, stopRight);
+ }
+
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 (debug)
{
- if (!d->leftPu && !d->rightPu)
+ count++;
+ printf("tid: %d, %d\n", d->leftTid, d->rightTid);
+ printf("pos: %d, %d\n", d->leftPos, d->rightPos);
+ printf("stop: %d, %d\n", stopLeft, stopRight);
+ }
+
+ while (stopLeft && stopRight && count <= 50)
{
- stopLeft = bam_dual_plbuf_resume(bufLeft);
- stopRight = bam_dual_plbuf_resume(bufRight);
+ if (debug)
+ {
+ count++;
+ printf("loop tid: %d, %d\n", d->leftTid, d->rightTid);
+ printf("loop pos: %d, %d\n", d->leftPos, d->rightPos);
+ printf("loop stop: %d, %d\n", stopLeft, stopRight);
}
- else if (!d->leftPu)
+// if (!d->leftPu && !d->rightPu)
+// {
+// stopLeft = bam_dual_plbuf_resume(bufLeft);
+// stopRight = bam_dual_plbuf_resume(bufRight);
+// }
+ if (!d->leftPu)
stopLeft = bam_dual_plbuf_resume(bufLeft);
else if (!d->rightPu)
stopRight = bam_dual_plbuf_resume(bufRight);
+ else if (d->leftTid < d->rightTid)
+ stopLeft = bam_dual_plbuf_resume(bufLeft);
+ else if (d->leftTid > d->rightTid)
+ 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);
+// count++;
+ } while (retLeft >= 0 && retRight >= 0 && count <= 500);
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