43bbf45d9f7b380ea7394b7a8628be0ad3ae1b8e angie Tue Jul 5 16:10:06 2011 -0700 MLQ #4463 (problem w/intersection):enabling intersection for bigWig uncovered a couple new bugs: 1. The secondary-table-fetching code didn't handle the bigWig case. Fix: added bigWigIntervalsToBed and hooked it up in bitsForIntersectingTable. 2. The code for intersecting w/bigWig primary table had been written but never tested, and the code for incompletely covered items from the primary table was incorrect. diff --git src/hg/hgTables/bigWig.c src/hg/hgTables/bigWig.c index e13564e..ad80fdc 100644 --- src/hg/hgTables/bigWig.c +++ src/hg/hgTables/bigWig.c @@ -110,38 +110,37 @@ next = iv->next; int start = iv->start; int size = iv->end - start; int overlap = bitCountRange(bits2, start, size); if (isBpWise) { if (overlap == size) { slAddHead(&newList, iv); } else if (overlap > 0) { /* Here we have to break things up. */ double val = iv->val; struct bbiInterval *partIv = iv; // Reuse memory for first interval - int s = start, end = start+size; + int s = iv->start, end = iv->end; for (;;) { - int bitsClear = bitFindClear(bits2, s, end-s); - s += bitsClear; + s = bitFindSet(bits2, s, end); if (s >= end) break; - int bitsSet = bitFindSet(bits2, s, end-s); // Always > 0 + int bitsSet = bitFindClear(bits2, s, end) - s; if (partIv == NULL) lmAllocVar(lm, partIv); partIv->start = s; partIv->end = s + bitsSet; partIv->val = val; slAddHead(&newList, partIv); partIv = NULL; s += bitsSet; if (s >= end) break; } } } else { @@ -289,30 +288,52 @@ numberStatRow("bases with data", validCount); long long regionSize = basesInRegion(regionList,0); long long gapTotal = gapsInRegion(conn, regionList,0); numberStatRow("bases with sequence", regionSize - gapTotal); numberStatRow("bases in region", regionSize); wigFilterStatRow(conn); stringStatRow("intersection", cartUsualString(cart, hgtaIntersectTable, "off")); long wigFetchTime = clock1000() - startTime; floatStatRow("load and calc time", 0.001*wigFetchTime); hTableEnd(); bbiFileClose(&bwf); htmlClose(); } +struct bed *bigWigIntervalsToBed(struct sqlConnection *conn, char *table, struct region *region, + struct lm *lm) +/* Return a list of unfiltered, unintersected intervals in region as bed (for + * secondary table in intersection). */ +{ +struct bed *bed, *bedList = NULL; +char *fileName = bigWigFileName(table, conn); +struct bbiFile *bwf = bigWigFileOpen(fileName); +struct bbiInterval *iv, *ivList = bigWigIntervalQuery(bwf, region->chrom, region->start, + region->end, lm); +for (iv = ivList; iv != NULL; iv = iv->next) + { + lmAllocVar(lm, bed); + bed->chrom = region->chrom; + bed->chromStart = iv->start; + bed->chromEnd = iv->end; + slAddHead(&bedList, bed); + } +slReverse(&bedList); +return bedList; +} + void bigWigFillDataVector(char *table, struct region *region, struct sqlConnection *conn, struct dataVector *vector) /* Fill in data vector with bigWig info on region. Handles filters and intersections. */ { /* Figure out filter values if any. */ double ll, ul; enum wigCompare cmp; getWigFilter(database, curTable, &cmp, &ll, &ul); /* Get intervals that pass filter and intersection. */ struct lm *lm = lmInit(0); char *fileName = bigWigFileName(table, conn); struct bbiFile *bwf = bigWigFileOpen(fileName); struct bbiInterval *iv, *ivList; ivList = intersectedFilteredBbiIntervalsOnRegion(conn, bwf, region, cmp, ll, ul, lm);