a97681a357b6f9849619f7ff1794d180d69c525c kent Tue Mar 29 14:30:04 2011 -0700 Making bigWigAverageOverBed use the bigWigValsOnChrom interface. diff --git src/utils/bigWigAverageOverBed/bigWigAverageOverBed.c src/utils/bigWigAverageOverBed/bigWigAverageOverBed.c index b4d6628..2097e98 100644 --- src/utils/bigWigAverageOverBed/bigWigAverageOverBed.c +++ src/utils/bigWigAverageOverBed/bigWigAverageOverBed.c @@ -141,112 +141,81 @@ *pSumCoverage += cov1; int i; double sum1 = 0; for (i=start; i<end; ++i) sum1 += valBuf[i]; *pSumVal += sum1; } void averageFetchingEachChrom(struct bbiFile *bbi, struct bed **pBedList, int fieldCount, FILE *f) /* Do the averaging by sorting bedList by chromosome, and then processing each chromosome * at once. Faster for long bedLists. */ { /* Sort by chromosome. */ slSort(pBedList, bedCmpChrom); -double *valBuf = NULL; -Bits *covBuf = NULL; -int bufSize = 0; +struct bigWigValsOnChrom *chromVals = bigWigValsOnChromNew(); struct bed *bed, *bedList, *nextChrom; verbose(1, "processing chromosomes"); for (bedList = *pBedList; bedList != NULL; bedList = nextChrom) { /* Figure out which chromosome we're working on, and the last bed using it. */ char *chrom = bedList->chrom; nextChrom = nextChromInList(bedList); + verbose(2, "Processing %s\n", chrom); - /* Check that bigWig file actually has data on this chromosome. If not - * just output as if coverage is 0 */ - int chromSize = bbiChromSize(bbi, chrom); - if (chromSize <= 0) + if (bigWigValsOnChromFetchData(chromVals, chrom, bbi)) { - for (bed = bedList; bed != nextChrom; bed = bed->next) - fprintf(f, "%s\t%d\t0\t0\t0\t0\n", bed->name, bedTotalBlockSize(bed)); - } - else - { - verbose(2, "Processing %s (%d bases)\n", chrom, chromSize); - - /* Make sure merge buffers are big enough. */ - if (chromSize > bufSize) - { - bufSize = chromSize; - freeMem(covBuf); - freeMem(valBuf); - valBuf = needHugeMem(bufSize * sizeof(double)); - covBuf = bitAlloc(bufSize); - } - - /* Zero out buffers */ - bitClear(covBuf, chromSize); - int i; - for (i=0; i<chromSize; ++i) - valBuf[i] = 0.0; - - /* Fetch intervals for this chromosome and fold into buffers. */ - struct lm *lm = lmInit(0); - struct bbiInterval *iv, *ivList = bigWigIntervalQuery(bbi, chrom, 0, chromSize, lm); - for (iv = ivList; iv != NULL; iv = iv->next) - { - double val = iv->val; - int end = iv->end; - for (i=iv->start; i<end; ++i) - valBuf[i] = val; - bitSetRange(covBuf, iv->start, iv->end - iv->start); - } - lmCleanup(&lm); + double *valBuf = chromVals->valBuf; + Bits *covBuf = chromVals->covBuf; /* Loop through beds doing sums and outputting. */ for (bed = bedList; bed != nextChrom; bed = bed->next) { int size = 0, coverage = 0; double sum = 0.0; if (fieldCount < 12) { addBufIntervalInfo(valBuf, covBuf, bed->chromStart, bed->chromEnd, &size, &coverage, &sum); } else { int i; for (i=0; i<bed->blockCount; ++i) { int start = bed->chromStart + bed->chromStarts[i]; int end = start + bed->blockSizes[i]; addBufIntervalInfo(valBuf, covBuf, start, end, &size, &coverage, &sum); } } /* Print out result, fudging mean to 0 if no coverage at all. */ double mean = 0; if (coverage > 0) mean = sum/coverage; fprintf(f, "%s\t%d\t%d\t%g\t%g\t%g\n", bed->name, size, coverage, sum, sum/size, mean); } verboseDot(); } + else + { + /* If no bigWig data on this chromosome, just output as if coverage is 0 */ + for (bed = bedList; bed != nextChrom; bed = bed->next) + fprintf(f, "%s\t%d\t0\t0\t0\t0\n", bed->name, bedTotalBlockSize(bed)); + } } verbose(1, "\n"); } void bigWigAverageOverBed(char *inBw, char *inBed, char *outTab) /* bigWigAverageOverBed - Compute average score of big wig over each bed, which may have introns. */ { struct bed *bedList; int fieldCount; bedLoadAllReturnFieldCount(inBed, &bedList, &fieldCount); checkUniqueNames(bedList); struct bbiFile *bbi = bigWigFileOpen(inBw); FILE *f = mustOpen(outTab, "w");