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; ichrom; 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; inext) - { - double val = iv->val; - int end = iv->end; - for (i=iv->start; istart, 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; iblockCount; ++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");