CNVkit segment confidence values & weighted segment mean

cnv_segment_stats.py

from __future__ import print_function

import sys

import numpy as np
import cnvlib
from cnvlib.metrics import biweight_location

cnr, cns = sys.argv[1:]

probes = cnvlib.read(cnr)
segments = cnvlib.read(cns)

# print("Segment", "nprobes", "CI-lo", "CI-hi", "orig",
#       "mean", "wt.mean", "biloc", sep='\t')
out_cns_mean = []
out_cns_ci = []
for seg, segcna in probes.by_segment(segments):
    label = "{}:{}-{}".format(seg["chromosome"], seg["start"], seg["end"])
    k = len(segcna)
    if k == 0:
        print("Skipping 0-probe segment", label, file=sys.stderr)
        continue

    # Bootstrap for CI
    bootstrap_dist = [np.random.choice(segcna["coverage"], k).mean()
                      for i in xrange(100)]
    ci = np.percentile(bootstrap_dist, [2.5, 97.5])
    wt_mean = np.average(segcna.coverage, weights=segcna["weight"])
    # print(label, k,
    #       ci[0], ci[1],
    #       seg["coverage"],
    #       # Recalculate segment value
    #       segcna.coverage.mean(),
    #       wt_mean,
    #       biweight_location(segcna.coverage),
    #       sep='\t')
    out_cns_mean.append(wt_mean)
    out_cns_ci.append(ci)

segments["coverage"] = np.asarray(out_cns_mean)
segments["gene"] = np.asarray(["CI=%.3g_%.3g" % tuple(ci) for ci in out_cns_ci])
segments.write("out.cns")