dkfzbiasfilter_summarize.py 文件源码

def main(filter_vcf, opts):
    changes_pcr = collections.defaultdict(int)
    changes_seq = collections.defaultdict(int)
    depths_pcr = collections.defaultdict(int)
    depths_seq = collections.defaultdict(int)
    depths_pass = collections.defaultdict(int)
    with pysam.VariantFile(filter_vcf) as bcf_in:
        damage_filters = set([f.name for f in bcf_in.header.filters.values()
                              if f.description.startswith("Variant allele shows a bias")])
        for rec in bcf_in.fetch():
            if len(set(rec.filter) & damage_filters) == len(rec.filter):
                if len(rec.ref) == 1 and len(rec.alts[0]) == 1:
                    change = "%s->%s" % (rec.ref, rec.alts[0])
                    if "bPcr" in set(rec.filter):
                        changes_pcr[change] += 1
                        depths_pcr = damage_support(rec, depths_pcr)
                    if "bSeq" in set(rec.filter):
                        changes_seq[change] += 1
                        seqerror_support(rec, depths_seq)
            elif list(rec.filter) == ["PASS"]:
                depths_pass = pass_support(rec, depths_pass)
    changes = {}
    depths = {}
    for val, change in _organize_changes(changes_pcr, opts):
        changes["%s damage" % change] = val
    depths["damage"] = {}
    for depth, count in sorted(depths_pcr.items()):
        if depth < opts.maxdepth:
            depths["damage"][depth] = count
    for val, change in _organize_changes(changes_seq, opts):
        changes["%s bias" % change] = val
    depths["bias"] = {}
    for depth, count in sorted(depths_seq.items()):
        if depth < opts.maxdepth:
            depths["bias"][depth] = count
    depths["pass"] = {}
    for depth, count in sorted(depths_pass.items()):
        if depth < opts.maxdepth:
            depths["pass"][depth] = count
    out = {"changes": changes, "depths": depths}
    if opts.sample:
        out["sample"] = opts.sample
    out_handle = open(opts.outfile, "w") if opts.outfile else sys.stdout
    yaml.safe_dump(out, out_handle, default_flow_style=False, allow_unicode=False)