bwStats.c

#include "bigWig.h"
#include "bwCommon.h"
#include <errno.h>
#include <stdlib.h>
#include <zlib.h>
#include <math.h>
#include <string.h>

//Returns -1 if there are no applicable levels, otherwise an integer indicating the most appropriate level.
//Like Kent's library, this divides the desired bin size by 2 to minimize the effect of blocks overlapping multiple bins
static int32_t determineZoomLevel(const bigWigFile_t *fp, int basesPerBin) {
    int32_t out = -1;
    int64_t diff;
    uint32_t bestDiff = -1;
    uint16_t i;

    basesPerBin/=2;
    for(i=0; i<fp->hdr->nLevels; i++) {
        diff = basesPerBin - (int64_t) fp->hdr->zoomHdrs->level[i];
        if(diff >= 0 && diff < bestDiff) {
            bestDiff = diff;
            out = i;
        }
    }
    return out;
}

/// @cond SKIP
struct val_t {
    uint32_t nBases;
    float min, max, sum, sumsq;
    double scalar;
};

struct vals_t {
    uint32_t n;
    struct val_t **vals;
};
/// @endcond

void destroyVals_t(struct vals_t *v) {
    uint32_t i;
    if(!v) return;
    for(i=0; i<v->n; i++) free(v->vals[i]);
    if(v->vals) free(v->vals);
    free(v);
}

//Determine the base-pair overlap between an interval and a block
double getScalar(uint32_t i_start, uint32_t i_end, uint32_t b_start, uint32_t b_end) {
    double rv = 0.0;
    if(b_start <= i_start) {
        if(b_end > i_start) rv = ((double)(b_end - i_start))/(b_end-b_start);
    } else if(b_start < i_end) {
        if(b_end < i_end) rv = ((double)(b_end - b_start))/(b_end-b_start);
        else rv = ((double)(i_end - b_start))/(b_end-b_start);
    }

    return rv;
}

//Returns NULL on error
static struct vals_t *getVals(bigWigFile_t *fp, bwOverlapBlock_t *o, int i, uint32_t tid, uint32_t start, uint32_t end) {
    void *buf = NULL, *compBuf = NULL;
    uLongf sz = fp->hdr->bufSize;
    int compressed = 0, rv;
    uint32_t *p, vtid, vstart, vend;
    struct vals_t *vals = NULL;
    struct val_t *v = NULL;

    if(sz) {
        compressed = 1;
        buf = malloc(sz); 
    }
    sz = 0; //This is now the size of the compressed buffer

    if(bwSetPos(fp, o->offset[i])) goto error;

    vals = calloc(1,sizeof(struct vals_t));
    if(!vals) goto error;

    v = malloc(sizeof(struct val_t));
    if(!v) goto error;

    if(sz < o->size[i]) compBuf = malloc(o->size[i]);
    if(!compBuf) goto error;

    if(bwRead(compBuf, o->size[i], 1, fp) != 1) goto error;
    if(compressed) {
        sz = fp->hdr->bufSize;
        rv = uncompress(buf, &sz, compBuf, o->size[i]);
        if(rv != Z_OK) goto error;
    } else {
        buf = compBuf;
        sz = o->size[i];
    }

    p = buf;
    while(((uLongf) ((void*)p-buf)) < sz) {
        vtid = p[0];
        vstart = p[1];
        vend = p[2];
        v->nBases = p[3];
        v->min = ((float*) p)[4];
        v->max = ((float*) p)[5];
        v->sum = ((float*) p)[6];
        v->sumsq = ((float*) p)[7];
        v->scalar = getScalar(start, end, vstart, vend);

        if(tid == vtid) {
            if((start <= vstart && end > vstart) || (start < vend && start >= vstart)) {
                vals->vals = realloc(vals->vals, sizeof(struct val_t*)*(vals->n+1));
                if(!vals->vals) goto error;
                vals->vals[vals->n++] = v;
                v = malloc(sizeof(struct val_t));
                if(!v) goto error;
            }
            if(vstart > end) break;
        } else if(vtid > tid) {
            break;
        }
        p+=8;
    }

    free(v);
    free(buf);
    if(compressed) free(compBuf);
    return vals;

error:
    if(buf) free(buf);
    if(compBuf && compressed) free(compBuf);
    if(v) free(v);
    destroyVals_t(vals);
    return NULL;
}

//On error, errno is set to ENOMEM and NaN is returned (though NaN can be returned normally)
static double blockMean(bigWigFile_t *fp, bwOverlapBlock_t *blocks, uint32_t tid, uint32_t start, uint32_t end) {
    uint32_t i, j;
    double output = 0.0, coverage = 0.0;
    struct vals_t *v = NULL;

    if(!blocks->n) return strtod("NaN", NULL);

    //Iterate over the blocks
    for(i=0; i<blocks->n; i++) {
        v = getVals(fp, blocks, i, tid, start, end);
        if(!v) goto error;
        for(j=0; j<v->n; j++) {
            output += v->vals[j]->sum * v->vals[j]->scalar;
            coverage += v->vals[j]->nBases * v->vals[j]->scalar;
        }
        destroyVals_t(v);
    }


    if(!coverage) return strtod("NaN", NULL);

    return output/coverage;

error:
    if(v) free(v);
    errno = ENOMEM;
    return strtod("NaN", NULL);
}

static double intMean(bwOverlappingIntervals_t* ints, uint32_t start, uint32_t end) {
    double sum = 0.0;
    uint32_t nBases = 0, i, start_use, end_use;

    if(!ints->l) return strtod("NaN", NULL);

    for(i=0; i<ints->l; i++) {
        start_use = ints->start[i];
        end_use = ints->end[i];
        if(ints->start[i] < start) start_use = start;
        if(ints->end[i] > end) end_use = end;
        nBases += end_use-start_use;
        sum += (end_use-start_use)*((double) ints->value[i]);
    }

    return sum/nBases;
}

//Does UCSC compensate for partial block/range overlap?
static double blockDev(bigWigFile_t *fp, bwOverlapBlock_t *blocks, uint32_t tid, uint32_t start, uint32_t end) {
    uint32_t i, j;
    double mean = 0.0, ssq = 0.0, coverage = 0.0, diff;
    struct vals_t *v = NULL;

    if(!blocks->n) return strtod("NaN", NULL);

    //Iterate over the blocks
    for(i=0; i<blocks->n; i++) {
        v = getVals(fp, blocks, i, tid, start, end);
        if(!v) goto error;
        for(j=0; j<v->n; j++) {
            coverage += v->vals[j]->nBases * v->vals[j]->scalar;
            mean += v->vals[j]->sum * v->vals[j]->scalar;
            ssq += v->vals[j]->sumsq * v->vals[j]->scalar;
        }
        destroyVals_t(v);
        v = NULL;
    }

    if(coverage<=1.0) return strtod("NaN", NULL);
    diff = ssq-mean*mean/coverage;
    if(coverage > 1.0) diff /= coverage-1;
    if(fabs(diff) > 1e-8) { //Ignore floating point differences
        return sqrt(diff);
    } else {
        return 0.0;
    }

error:
    if(v) destroyVals_t(v);
    errno = ENOMEM;
    return strtod("NaN", NULL);
}

//This uses compensated summation to account for finite precision math
static double intDev(bwOverlappingIntervals_t* ints, uint32_t start, uint32_t end) {
    double v1 = 0.0, mean, rv;
    uint32_t nBases = 0, i, start_use, end_use;

    if(!ints->l) return strtod("NaN", NULL);
    mean = intMean(ints, start, end);

    for(i=0; i<ints->l; i++) {
        start_use = ints->start[i];
        end_use = ints->end[i];
        if(ints->start[i] < start) start_use = start;
        if(ints->end[i] > end) end_use = end;
        nBases += end_use-start_use;
        v1 += (end_use-start_use) * pow(ints->value[i]-mean, 2.0); //running sum of squared difference
    }

    if(nBases>=2) rv = sqrt(v1/(nBases-1));
    else if(nBases==1) rv = sqrt(v1);
    else rv = strtod("NaN", NULL);

    return rv;
}

static double blockMax(bigWigFile_t *fp, bwOverlapBlock_t *blocks, uint32_t tid, uint32_t start, uint32_t end) {
    uint32_t i, j, isNA = 1;
    double o = strtod("NaN", NULL);
    struct vals_t *v = NULL;

    if(!blocks->n) return o;

    //Iterate the blocks
    for(i=0; i<blocks->n; i++) {
        v = getVals(fp, blocks, i, tid, start, end);
        if(!v) goto error;
        for(j=0; j<v->n; j++) {
            if(isNA) {
                o = v->vals[j]->max;
                isNA = 0;
            } else if(v->vals[j]->max > o) {
                o = v->vals[j]->max;
            }
        }
        destroyVals_t(v);
    }

    return o;

error:
    destroyVals_t(v);
    errno = ENOMEM;
    return strtod("NaN", NULL);
}

static double intMax(bwOverlappingIntervals_t* ints) {
    uint32_t i;
    double o;

    if(ints->l < 1) return strtod("NaN", NULL);

    o = ints->value[0];
    for(i=1; i<ints->l; i++) {
        if(ints->value[i] > o) o = ints->value[i];
    }

    return o;
}

static double blockMin(bigWigFile_t *fp, bwOverlapBlock_t *blocks, uint32_t tid, uint32_t start, uint32_t end) {
    uint32_t i, j, isNA = 1;
    double o = strtod("NaN", NULL);
    struct vals_t *v = NULL;

    if(!blocks->n) return o;

    //Iterate the blocks
    for(i=0; i<blocks->n; i++) {
        v = getVals(fp, blocks, i, tid, start, end);
        if(!v) goto error;
        for(j=0; j<v->n; j++) {
            if(isNA) {
                o = v->vals[j]->min;
                isNA = 0;
            } else if(v->vals[j]->min < o) o = v->vals[j]->min;
        }
        destroyVals_t(v);
    }

    return o;

error:
    destroyVals_t(v);
    errno = ENOMEM;
    return strtod("NaN", NULL);
}

static double intMin(bwOverlappingIntervals_t* ints) {
    uint32_t i;
    double o;

    if(ints->l < 1) return strtod("NaN", NULL);

    o = ints->value[0];
    for(i=1; i<ints->l; i++) {
        if(ints->value[i] < o) o = ints->value[i];
    }

    return o;
}

//Does UCSC compensate for only partial block/interval overlap?
static double blockCoverage(bigWigFile_t *fp, bwOverlapBlock_t *blocks, uint32_t tid, uint32_t start, uint32_t end) {
    uint32_t i, j;
    double o = 0.0;
    struct vals_t *v = NULL;

    if(!blocks->n) return strtod("NaN", NULL);

    //Iterate over the blocks
    for(i=0; i<blocks->n; i++) {
        v = getVals(fp, blocks, i, tid, start, end);
        if(!v) goto error;
        for(j=0; j<v->n; j++) {
            o+= v->vals[j]->nBases * v->vals[j]->scalar;
        }
        destroyVals_t(v);
    }

    if(o == 0.0) return strtod("NaN", NULL);
    return o;

error:
    destroyVals_t(v);
    errno = ENOMEM;
    return strtod("NaN", NULL);
}

static double intCoverage(bwOverlappingIntervals_t* ints, uint32_t start, uint32_t end) {
    uint32_t i, start_use, end_use;
    double o = 0.0;

    if(!ints->l) return strtod("NaN", NULL);

    for(i=0; i<ints->l; i++) {
        start_use = ints->start[i];
        end_use = ints->end[i];
        if(start_use < start) start_use = start;
        if(end_use > end) end_use = end;
        o += end_use - start_use;
    }

    return o/(end-start);
}

static double blockSum(bigWigFile_t *fp, bwOverlapBlock_t *blocks, uint32_t tid, uint32_t start, uint32_t end) {
    uint32_t i, j, sizeUse;
    double o = 0.0;
    struct vals_t *v = NULL;

    if(!blocks->n) return strtod("NaN", NULL);

    //Iterate over the blocks
    for(i=0; i<blocks->n; i++) {
        v = getVals(fp, blocks, i, tid, start, end);
        if(!v) goto error;
        for(j=0; j<v->n; j++) {
            //Multiply the block average by min(bases covered, block overlap with interval)
            sizeUse = v->vals[j]->scalar;
            if(sizeUse > v->vals[j]->nBases) sizeUse = v->vals[j]->nBases;
            o+= (v->vals[j]->sum * sizeUse) / v->vals[j]->nBases;
        }
        destroyVals_t(v);
    }

    if(o == 0.0) return strtod("NaN", NULL);
    return o;

error:
    destroyVals_t(v);
    errno = ENOMEM;
    return strtod("NaN", NULL);
}

static double intSum(bwOverlappingIntervals_t* ints, uint32_t start, uint32_t end) {
    uint32_t i, start_use, end_use;
    double o = 0.0;

    if(!ints->l) return strtod("NaN", NULL);

    for(i=0; i<ints->l; i++) {
        start_use = ints->start[i];
        end_use = ints->end[i];
        if(start_use < start) start_use = start;
        if(end_use > end) end_use = end;
        o += (end_use - start_use) * ints->value[i];
    }

    return o;
}

//Returns NULL on error, otherwise a double* that needs to be free()d
static double *bwStatsFromZoom(bigWigFile_t *fp, int32_t level, uint32_t tid, uint32_t start, uint32_t end, uint32_t nBins, enum bwStatsType type) {
    bwOverlapBlock_t *blocks = NULL;
    double *output = NULL;
    uint32_t pos = start, i, end2;

    if(!fp->hdr->zoomHdrs->idx[level]) {
        fp->hdr->zoomHdrs->idx[level] = bwReadIndex(fp, fp->hdr->zoomHdrs->indexOffset[level]);
        if(!fp->hdr->zoomHdrs->idx[level]) return NULL;
    }
    errno = 0; //Sometimes libCurls sets and then doesn't unset errno on errors

    output = malloc(sizeof(double)*nBins);
    if(!output) return NULL;

    for(i=0, pos=start; i<nBins; i++) {
        end2 = start + ((double)(end-start)*(i+1))/((int) nBins);
        blocks = walkRTreeNodes(fp, fp->hdr->zoomHdrs->idx[level]->root, tid, pos, end2);
        if(!blocks) goto error;

        switch(type) {
        case 0:
            //mean
            output[i] = blockMean(fp, blocks, tid, pos, end2);
            break;
        case 1:
            //stdev
            output[i] = blockDev(fp, blocks, tid, pos, end2);
            break;
        case 2:
            //max
            output[i] = blockMax(fp, blocks, tid, pos, end2);
            break;
        case 3:
            //min
            output[i] = blockMin(fp, blocks, tid, pos, end2);
            break;
        case 4:
            //cov
            output[i] = blockCoverage(fp, blocks, tid, pos, end2)/(end2-pos);
            break;
        case 5:
            //sum
            output[i] = blockSum(fp, blocks, tid, pos, end2);
            break;
        default:
            goto error;
            break;
        }
        if(errno) goto error;
        destroyBWOverlapBlock(blocks);
        pos = end2;
    }

    return output;

error:
    fprintf(stderr, "got an error in bwStatsFromZoom in the range %"PRIu32"-%"PRIu32": %s\n", pos, end2, strerror(errno));
    if(blocks) destroyBWOverlapBlock(blocks);
    if(output) free(output);
    return NULL;
}

double *bwStatsFromFull(bigWigFile_t *fp, const char *chrom, uint32_t start, uint32_t end, uint32_t nBins, enum bwStatsType type) {
    bwOverlappingIntervals_t *ints = NULL;
    double *output = malloc(sizeof(double)*nBins);
    uint32_t i, pos = start, end2;
    if(!output) return NULL;

    for(i=0; i<nBins; i++) {
        end2 = start + ((double)(end-start)*(i+1))/((int) nBins);
        ints = bwGetOverlappingIntervals(fp, chrom, pos, end2);

        if(!ints) {
            output[i] = strtod("NaN", NULL);
            continue;
        }

        switch(type) {
        default :
        case 0:
            output[i] = intMean(ints, pos, end2);
            break;
        case 1:
            output[i] = intDev(ints, pos, end2);
            break;
        case 2:
            output[i] = intMax(ints);
            break;
        case 3:
            output[i] = intMin(ints);
            break;
        case 4:
            output[i] = intCoverage(ints, pos, end2);
            break;
        case 5:
            output[i] = intSum(ints, pos, end2);
            break;
        }
        bwDestroyOverlappingIntervals(ints);
        pos = end2;
    }

    return output;
}

//Returns a list of floats of length nBins that must be free()d
//On error, NULL is returned
double *bwStats(bigWigFile_t *fp, const char *chrom, uint32_t start, uint32_t end, uint32_t nBins, enum bwStatsType type) {
    int32_t level = determineZoomLevel(fp, ((double)(end-start))/((int) nBins));
    uint32_t tid = bwGetTid(fp, chrom);
    if(tid == (uint32_t) -1) return NULL;

    if(level == -1) return bwStatsFromFull(fp, chrom, start, end, nBins, type);
    return bwStatsFromZoom(fp, level, tid, start, end, nBins, type);
}