@@ -0,0 +1,1675 @@
+%22%22%22%0Acompare what we get from tabix with what we get from the encoded CADD stuff%0Aand make sure it is %3C 0.05 difference%0A%22%22%22%0Aimport sys%0Aimport toolshed as ts%0Aimport os.path as op%0A%0Acadd_prefix = sys.argv%5B1%5D%0Atsv = sys.argv%5B2%5D%0A%0Aassert op.exists(%22%25s.idx%22 %25 cadd_prefix)%0Aassert op.exists(%22%25s.bin%22 %25 cadd_prefix)%0A%0Aslen = %7Bx%5B0%5D: int(x%5B1%5D) for x in (l.rstrip().split() for l in open(cadd_prefix + %22.idx%22))%7D%0A%0A%0Adef read_cadd(path, chrom, pos):%0A    cmd = %22cadd %7Bpath%7D %7Bchrom%7D %7Bpos%7D %7Bletter%7D%22%0A    cmds = %5B%5D%0A    d = locals()%0A    for letter in %22ACGT%22:%0A        d%5B'letter'%5D = letter%0A        cmds.append(cmd.format(**d))%0A    cmd = %22 && %22.join(cmds)%0A    vals = %5Bx.strip().split() for x in ts.nopen('%7C' + cmd)%5D%0A    for i in range(len(vals)):%0A        vals%5Bi%5D%5B0%5D = float(vals%5Bi%5D%5B0%5D)%0A    return dict(zip(%22ACGT%22, vals))%0A%0Adef read_tsv(path, chrom, pos):%0A    cmd = %22tabix %7Bpath%7D %7Bchrom%7D:%7Bpos%7D-%7Bpos%7D%22%0A    res = %5Bx.strip().split() for x in ts.nopen(%22%7C%25s%22 %25 cmd.format(**locals()))%5D%0A    v = %7Br%5B3%5D: float(r%5B5%5D) for r in res%7D%0A    for k in %22ACGT%22:%0A        if not k in v:%0A            v%5Bk%5D = 0.0%0A    return v%0A%0Adef runner(chrom, length):%0A    step = max(2, length/10000)%0A    for i, pos in enumerate(range(10000, length, step)):%0A%0A        cadd = read_cadd(cadd_prefix + %22.idx%22, chrom, pos)%0A        tabx = read_tsv(tsv, chrom, pos)%0A        for letter in %22ACGT%22:%0A%0A            if abs(cadd%5Bletter%5D%5B0%5D - tabx%5Bletter%5D) %3E 0.05:%0A                print %22BAD:%22, chrom, pos, cadd, tabx%0A                1/0%0A%0A        else:%0A            if i %25 100 == 0:%0A                print %22%25s %25d/%25d ok%22 %25 (chrom, pos, length)%0A%0Alist(ts.pmap(runner, ((chrom, length) for chrom, length in slen.items() if chrom.isdigit() or chrom in %22XY%22), n=12))%0A

add cadd test script