allele balance invenstigation in CEPH

README.md

nim c -d:danger -d:release -r abs.nim > abs.txt
mkdir -p pngs
python abs.py
rm -f pngs/NA*
rm -f pngs/HG*
convert -delay 25 pngs/*.png -loop 0 abs.gif

abs.nim

import hts/vcf

let path = "/scratch/ucgd/lustre-work/u1006375/ceph/vcf/FQF-1.2.1_16-08-06_WashU-Yandell-CEPH_1000Genomes_1512787121.bcf"

var ivcf:VCF
if not ivcf.open(path, threads=3):
  quit "couldn't open vcf"

const N = 50
var ab_hist = newSeq[array[N + 1, uint32]](ivcf.n_samples)

proc get_ab(v:Variant, ints: var seq[int32]): seq[float32] =
  doAssert v.format.get("AD", ints) == Status.OK
  result = newSeq[float32](v.n_samples)

  for i, ab in result.mpairs:
    ab = max(0, ints[2*i + 1]) / (max(0, ints[2*i+1]) + max(0, ints[2*i]))


var vqslod: seq[float32]
var ints: seq[int32]
var x: seq[int32]

for v in ivcf.query("19"):
  if v.ALT.len != 1: continue
  if v.REF.len != 1: continue
  if v.FILTER != "PASS": continue

  if v.info.get("VQSLOD", vqslod) != Status.OK:
    continue
  if vqslod[0] < 1'f32: continue

  var alts = v.format.genotypes(x).alts

  var ab_vals = get_ab(v, ints)

  for i, ab in ab_vals:
    if alts[i] == 1:
      ab_hist[i][(N.float32*ab).int].inc

for i, h in ab_hist:
  echo ivcf.samples[i], " ", h

abs.py

from matplotlib import pyplot as plt
import numpy as np

import matplotlib.backends.backend_pdf
pdf = matplotlib.backends.backend_pdf.PdfPages("abs.pdf")

for line in open("abs.txt"):

    sample, vals = line.split(" ", 1)
    # drop 1kg samples
    if sample[0] in "HN": continue
    vals = [int(x.strip()) for x in vals.strip("[] \n").split(",")]
    xs = np.arange(0, 1, 1.0/len(vals))
    plt.step(xs, vals)
    plt.title("sample:" + sample)
    plt.xlabel("het allele balance")
    plt.ylabel("variant count")
    f = plt.gcf()
    pdf.savefig(f)
    f.savefig("pngs/" + sample + ".png")
    plt.close()

pdf.close()