Collapse together variant calls from multiple approaches into single data item with `variants`.
969349
def combine_multiple_callers(samples):
    
    by_bam = collections.OrderedDict()
    for data in (x[0] for x in samples):
        work_bam = tz.get_in(("combine", "work_bam", "out"), data, data.get("align_bam"))
        # For pre-computed VCF inputs, we don't have BAM files
        if not work_bam:
            work_bam = dd.get_sample_name(data)
        jointcaller = tz.get_in(("config", "algorithm", "jointcaller"), data)
        variantcaller = get_variantcaller(data)
        key = (multi.get_batch_for_key(data), work_bam)
        if key not in by_bam:
            by_bam[key] = []
        by_bam[key].append((variantcaller, jointcaller, data))
    out = []
    for callgroup in by_bam.values():
        ready_calls = []
        for variantcaller, jointcaller, data in callgroup:
            if variantcaller:
                cur = data.get("vrn_file_plus", {})
                cur.update({"variantcaller": variantcaller,
                            "vrn_file": data.get("vrn_file_orig") if jointcaller else data.get("vrn_file"),
                            "vrn_file_batch": data.get("vrn_file_batch") if not jointcaller else None,
                            "vrn_stats": data.get("vrn_stats"),
                            "validate": data.get("validate") if not jointcaller else None})
                if jointcaller:
                    cur["population"] = False
                ready_calls.append(cur)
            if jointcaller:
                cur = {"variantcaller": jointcaller,
                       "vrn_file": data.get("vrn_file"),
                       "vrn_file_batch": data.get("vrn_file_batch"),
                       "validate": data.get("validate"),
                       "do_upload": False}
                if not variantcaller:
                    cur["population"] = {"vcf": data.get("vrn_file")}
                ready_calls.append(cur)
            if not jointcaller and not variantcaller:
                ready_calls.append({"variantcaller": "precalled",
                                    "vrn_file": data.get("vrn_file"),
                                    "validate": data.get("validate"),
                                    "do_upload": False})
        final = callgroup[0][-1]
        def orig_variantcaller_order(x):
            try:
                return final["config"]["algorithm"]["orig_variantcaller"].index(x["variantcaller"])
            except ValueError:
                return final["config"]["algorithm"]["orig_jointcaller"].index(x["variantcaller"])
        if len(ready_calls) > 1 and "orig_variantcaller" in final["config"]["algorithm"]:
            final["variants"] = sorted(ready_calls, key=orig_variantcaller_order)
            final["config"]["algorithm"]["variantcaller"] = final["config"]["algorithm"].pop("orig_variantcaller")
            if "orig_jointcaller" in final["config"]["algorithm"]:
                final["config"]["algorithm"]["jointcaller"] = final["config"]["algorithm"].pop("orig_jointcaller")
        else:
            final["variants"] = ready_calls
        final.pop("vrn_file_batch", None)
        final.pop("vrn_file_orig", None)
        final.pop("vrn_file_plus", None)
        final.pop("vrn_stats", None)
        out.append([final])
    return out
1032462
def cli(variant_file, vep, split, outfile, verbose, silent, check_info,
        allele_symbol, logfile, loglevel):
    
    from vcf_parser import logger, init_log

    if not loglevel:
        if verbose:
            loglevel = 'INFO'

    init_log(logger, logfile, loglevel)
    nr_of_variants = 0
    start = datetime.now()
    
    # with open(variant_file, 'r', encoding="utf-8") as f:
    #     for line in f:
    #         if not line.startswith('#'):
    #             nr_of_variants += 1
    
    if variant_file == '-':
        logger.info("Start parsing variants from stdin")
        my_parser = VCFParser(
            fsock=sys.stdin, 
            split_variants=split,
            check_info=check_info, 
            allele_symbol=allele_symbol
        )
    else:
        logger.info("Start parsing variants from file {0}".format(variant_file))
        my_parser = VCFParser(
            infile = variant_file,
            split_variants=split, 
            check_info=check_info, 
            allele_symbol=allele_symbol
        )

    if outfile:
        f = open(outfile, 'w', encoding='utf-8')
        logger.info("Printing vcf to file {0}".format(outfile))
    
    if not silent:
        logger.info("Printing vcf to stdout")
    else:
        logger.info("Skip printing since silent is active")
    
    for line in my_parser.metadata.print_header():
        if outfile:
            f.write(line+'\n')
        else:
            if not silent:
                print(line)
    try:
        for variant in my_parser:
            variant_line = '\t'.join([variant[head] for head in my_parser.header])
            if outfile:
                f.write(variant_line + '\n')
            else:
                if not silent:
                    print(variant_line)
            nr_of_variants += 1
    except SyntaxError as e:
        print(e)

    logger.info('Number of variants: {0}'.format(nr_of_variants))
    logger.info('Time to parse file: {0}'.format(str(datetime.now() - start)))
1129037
def _create_variant(self,
                        start,
                        end,
                        ref,
                        alt,
                        fsext_len=None,
                        is_dup=False,
                        acc=None,
                        is_ambiguous=False,
                        is_sub=False,
                        is_ext=False,
                        is_no_protein=False,
                        is_init_met=False):
        """"""

        if is_init_met:
            posedit = AARefAlt(ref=ref, alt=alt, init_met=True)
        elif is_ambiguous:
            posedit = None
        else:
            interval = Interval(start=start, end=end)
            # Note - order matters
            if is_no_protein:
                edit = '0'
            elif is_sub:
                edit = AASub(ref=ref, alt=alt)
            elif is_ext:
                edit = AAExt(ref=ref, alt=alt, aaterm='*', length=fsext_len)
            elif self._is_frameshift:
                edit = AAFs(ref=ref, alt=alt, length=fsext_len)
            elif is_dup:
                edit = Dup()
            elif ref == alt == '':
                edit = AARefAlt(ref='', alt='')
            else:
                edit = AARefAlt(ref=ref, alt=alt)
            posedit = PosEdit(
                pos=interval,
                edit=edit,
                uncertain=hgvs.global_config.mapping.inferred_p_is_uncertain)
        var_p = hgvs.sequencevariant.SequenceVariant(acc, 'p', posedit)
        return var_p
842500
def _get_variants(data):
    
    active_vs = []
    if "variants" in data:
        variants = data["variants"]
        # CWL based list of variants
        if isinstance(variants, dict) and "samples" in variants:
            variants = variants["samples"]
        for v in variants:
            # CWL -- a single variant file
            if isinstance(v, six.string_types) and os.path.exists(v):
                active_vs.append(_add_filename_details(v))
            elif (isinstance(v, (list, tuple)) and len(v) > 0 and
                  isinstance(v[0], six.string_types) and os.path.exists(v[0])):
                for subv in v:
                    active_vs.append(_add_filename_details(subv))
            elif isinstance(v, dict):
                if v.get("vrn_file"):
                    active_vs.append(v)
                elif v.get("population"):
                    vrnfile = v.get("population").get("vcf")
                    active_vs.append(_add_filename_details(vrnfile))
                elif v.get("vcf"):
                    active_vs.append(_add_filename_details(v.get("vcf")))
    return active_vs
236239
def _finalize_sv(solution_file, data):
    
    out = {"variantcaller": "titancna"}
    with open(solution_file) as in_handle:
        solution = dict(zip(in_handle.readline().strip("\r\n").split("\t"),
                            in_handle.readline().strip("\r\n").split("\t")))
    if solution.get("path"):
        out["purity"] = solution["purity"]
        out["ploidy"] = solution["ploidy"]
        out["cellular_prevalence"] = [x.strip() for x in solution["cellPrev"].split(",")]
        base = os.path.basename(solution["path"])
        out["plot"] = dict([(n, solution["path"] + ext) for (n, ext) in [("rplots", ".Rplots.pdf"),
                                                                         ("cf", "/%s_CF.pdf" % base),
                                                                         ("cna", "/%s_CNA.pdf" % base),
                                                                         ("loh", "/%s_LOH.pdf" % base)]
                            if os.path.exists(solution["path"] + ext)])
        out["subclones"] = "%s.segs.txt" % solution["path"]
        out["hetsummary"] = solution_file
        out["vrn_file"] = to_vcf(out["subclones"], "TitanCNA", _get_header, _seg_to_vcf, data)
        out["lohsummary"] = loh.summary_status(out, data)
    return out
542809
def parallel_combine_variants(orig_files, out_file, ref_file, config, run_parallel):
    
    file_key = "vcf_files"
    def split_by_region(data):
        base, ext = utils.splitext_plus(os.path.basename(out_file))
        args = []
        for region in [x.name for x in ref.file_contigs(ref_file, config)]:
            region_out = os.path.join(os.path.dirname(out_file), "%s-regions" % base,
                                      "%s-%s%s" % (base, region, ext))
            utils.safe_makedir(os.path.dirname(region_out))
            args.append((region_out, ref_file, config, region))
        return out_file, args
    config = copy.deepcopy(config)
    config["file_key"] = file_key
    prep_files = run_multicore(p_bgzip_and_index, [[x, config] for x in orig_files], config)
    items = [[{file_key: prep_files}]]
    parallel_split_combine(items, split_by_region, run_parallel,
                           "merge_variant_files", "concat_variant_files",
                           file_key, ["region", "sam_ref", "config"], split_outfile_i=0)
    return out_file
95450
def postprocess_variants(items):
    
    vrn_key = "vrn_file"
    if not isinstance(items, dict):
        items = [utils.to_single_data(x) for x in items]
        if "vrn_file_joint" in items[0]:
            vrn_key = "vrn_file_joint"
    data, items = _get_batch_representative(items, vrn_key)
    items = cwlutils.unpack_tarballs(items, data)
    data = cwlutils.unpack_tarballs(data, data)
    cur_name = "%s, %s" % (dd.get_sample_name(data), get_variantcaller(data, require_bam=False))
    logger.info("Finalizing variant calls: %s" % cur_name)
    orig_vrn_file = data.get(vrn_key)
    data = _symlink_to_workdir(data, [vrn_key])
    data = _symlink_to_workdir(data, ["config", "algorithm", "variant_regions"])
    if data.get(vrn_key):
        logger.info("Calculating variation effects for %s" % cur_name)
        ann_vrn_file, vrn_stats = effects.add_to_vcf(data[vrn_key], data)
        if ann_vrn_file:
            data[vrn_key] = ann_vrn_file
        if vrn_stats:
            data["vrn_stats"] = vrn_stats
        orig_items = _get_orig_items(items)
        logger.info("Annotate VCF file: %s" % cur_name)
        data[vrn_key] = annotation.finalize_vcf(data[vrn_key], get_variantcaller(data, require_bam=False),
                                                orig_items)
        if cwlutils.is_cwl_run(data):
            logger.info("Annotate with population level variation data")
            ann_file = population.run_vcfanno(data[vrn_key], data)
            if ann_file:
                data[vrn_key] = ann_file
        logger.info("Filtering for %s" % cur_name)
        data[vrn_key] = variant_filtration(data[vrn_key], dd.get_ref_file(data),
                                           tz.get_in(("genome_resources", "variation"), data, {}),
                                           data, orig_items)
        logger.info("Prioritization for %s" % cur_name)
        prio_vrn_file = prioritize.handle_vcf_calls(data[vrn_key], data, orig_items)
        if prio_vrn_file != data[vrn_key]:
            data[vrn_key] = prio_vrn_file
            logger.info("Germline extraction for %s" % cur_name)
            data = germline.extract(data, orig_items)

        if dd.get_align_bam(data):
            data = damage.run_filter(data[vrn_key], dd.get_align_bam(data), dd.get_ref_file(data),
                                     data, orig_items)
    if orig_vrn_file and os.path.samefile(data[vrn_key], orig_vrn_file):
        data[vrn_key] = orig_vrn_file
    return [[data]]
770607
def _handle_multiple_svcallers(data, stage):
    
    svs = get_svcallers(data)
    # special cases -- prioritization
    if stage == "ensemble" and dd.get_svprioritize(data):
        svs.append("prioritize")
    out = []
    for svcaller in svs:
        if svcaller in _get_callers([data], stage):
            base = copy.deepcopy(data)
            # clean SV callers present in multiple rounds and not this caller
            final_svs = []
            for sv in data.get("sv", []):
                if (stage == "ensemble" or sv["variantcaller"] == svcaller or sv["variantcaller"] not in svs
                      or svcaller not in _get_callers([data], stage, special_cases=True)):
                    final_svs.append(sv)
            base["sv"] = final_svs
            base["config"]["algorithm"]["svcaller"] = svcaller
            base["config"]["algorithm"]["svcaller_orig"] = svs
            out.append(base)
    return out
155376
def split_somatic(items):
    
    items = [_clean_flat_variantcaller(x) for x in items]
    somatic_groups, somatic, non_somatic = vcfutils.somatic_batches(items)
    # extract germline samples to run from normals in tumor/normal pairs
    germline_added = set([])
    germline = []
    for somatic_group in somatic_groups:
        paired = vcfutils.get_paired(somatic_group)
        if paired and paired.normal_data:
            cur = utils.deepish_copy(paired.normal_data)
            vc = dd.get_variantcaller(cur)
            if isinstance(vc, dict) and "germline" in vc:
                if cur["description"] not in germline_added:
                    germline_added.add(cur["description"])
                    cur["rgnames"]["sample"] = cur["description"]
                    cur["metadata"]["batch"] = "%s-germline" % cur["description"]
                    cur["metadata"]["phenotype"] = "germline"
                    cur = remove_align_qc_tools(cur)
                    cur["config"]["algorithm"]["variantcaller"] = vc["germline"]
                    germline.append(cur)
    # Fix variantcalling specification for only somatic targets
    somatic_out = []
    for data in somatic:
        vc = dd.get_variantcaller(data)
        if isinstance(vc, dict) and "somatic" in vc:
            data["config"]["algorithm"]["variantcaller"] = vc["somatic"]
        somatic_out.append(data)
    return non_somatic + somatic_out + germline
1095406
def trim_variant_sequences(variant_sequences, min_variant_sequence_coverage):
    
    n_total = len(variant_sequences)
    trimmed_variant_sequences = [
        variant_sequence.trim_by_coverage(min_variant_sequence_coverage)
        for variant_sequence in variant_sequences
    ]
    collapsed_variant_sequences = collapse_substrings(trimmed_variant_sequences)
    n_after_trimming = len(collapsed_variant_sequences)
    logger.info(
        "Kept %d/%d variant sequences after read coverage trimming to >=%dx",
        n_after_trimming,
        n_total,
        min_variant_sequence_coverage)
    return collapsed_variant_sequences
950799
def square_off(samples, run_parallel):
    
    to_process = []
    extras = []
    for data in [utils.to_single_data(x) for x in samples]:
        added = False
        if tz.get_in(("metadata", "batch"), data):
            for add in genotype.handle_multiple_callers(data, "jointcaller", require_bam=False):
                if _is_jointcaller_compatible(add):
                    added = True
                    to_process.append([add])
        if not added:
            extras.append([data])
    processed = grouped_parallel_split_combine(to_process, _split_by_callable_region,
                                               multi.group_batches_joint, run_parallel,
                                               "square_batch_region", "concat_variant_files",
                                               "vrn_file", ["region", "sam_ref", "config"])
    return _combine_to_jointcaller(processed) + extras
1232088
def visit_variants(self, func, variants=None, **kwargs):
        """"""
        if variants:
            present_variants = range(self.package.num_variants)
            invalid_variants = set(variants) - set(present_variants)
            if invalid_variants:
                raise BuildError(
                    "The package does not contain the variants: %s"
                    % ", ".join(str(x) for x in sorted(invalid_variants)))

        # iterate over variants
        results = []
        num_visited = 0

        for variant in self.package.iter_variants():
            if variants and variant.index not in variants:
                self._print_header(
                    "Skipping variant %s (%s)..."
                    % (variant.index, self._n_of_m(variant)))
                continue

            # visit the variant
            result = func(variant, **kwargs)
            results.append(result)
            num_visited += 1

        return num_visited, results
1321978
def variant_filtration(call_file, ref_file, vrn_files, data, items):
    
    caller = data["config"]["algorithm"].get("variantcaller")
    if "gvcf" not in dd.get_tools_on(data):
        call_file = ploidy.filter_vcf_by_sex(call_file, items)
    if caller in ["freebayes"]:
        return vfilter.freebayes(call_file, ref_file, vrn_files, data)
    elif caller in ["platypus"]:
        return vfilter.platypus(call_file, data)
    elif caller in ["samtools"]:
        return vfilter.samtools(call_file, data)
    elif caller in ["gatk", "gatk-haplotype", "haplotyper"]:
        if dd.get_analysis(data).lower().find("rna-seq") >= 0:
            from bcbio.rnaseq import variation as rnaseq_variation
            return rnaseq_variation.gatk_filter_rnaseq(call_file, data)
        else:
            return gatkfilter.run(call_file, ref_file, vrn_files, data)
    # no additional filtration for callers that filter as part of call process
    else:
        return call_file
1169798
def combine_calls_parallel(samples, run_parallel):
    
    batch_groups, extras = _group_by_batches(samples, _has_ensemble)
    out = []
    if batch_groups:
        processed = run_parallel("combine_calls", ((b, xs, xs[0]) for b, xs in batch_groups.items()))
        for batch_id, callinfo in processed:
            for data in batch_groups[batch_id]:
                data["variants"].insert(0, callinfo)
                out.append([data])
    return out + extras
237456
async def variant(self, elem=None, elem_type=None, params=None):
        
        elem_type = elem_type if elem_type else elem.__class__

        if hasattr(elem_type, 'blob_serialize'):
            elem = elem_type() if elem is None else elem
            return await elem.bob_serialize(self, elem=elem, elem_type=elem_type, params=params)

        if self.writing:
            return await self.dump_variant(elem=elem,
                                           elem_type=elem_type if elem_type else elem.__class__, params=params)
        else:
            return await self.load_variant(elem_type=elem_type if elem_type else elem.__class__,
                                           params=params, elem=elem)
593742
def variants_in_coords(nucmer_matches, vcf_file):
        ''''''
        found_variants = {}
        f = pyfastaq.utils.open_file_read(vcf_file)
        for line in f:
            if line.startswith('#'):
                continue

            data = line.rstrip().split('\t')
            scaff = data[0]

            if scaff in nucmer_matches:
                position = int(data[1]) - 1
                i = pyfastaq.intervals.Interval(position, position)
                intersects = len([x for x in nucmer_matches[scaff] if x.intersects(i)]) > 0
                if intersects:
                    if scaff not in found_variants:
                        found_variants[scaff] = set()
                    found_variants[scaff].add(position)

        pyfastaq.utils.close(f)
        return found_variants
151157
def batch(samples):
    
    samples = [utils.to_single_data(x) for x in samples]
    sample_order = [dd.get_sample_name(x) for x in samples]
    batch_groups = collections.defaultdict(list)
    for data in samples:
        batch_samples = tuple(data.get("batch_samples", [dd.get_sample_name(data)]))
        batch_groups[(batch_samples, dd.get_phenotype(data))].append(data)

    out = []
    for (batch_samples, phenotype), gsamples in batch_groups.items():
        if len(gsamples) > 1:
            batches = set([])
            for d in gsamples:
                batches |= set(dd.get_batches(d))
            gsamples.sort(key=dd.get_variantcaller_order)
            cur = copy.deepcopy(gsamples[0])
            cur.update({"batch_id": sorted(list(batches))[0] if batches else "_".join(batch_samples),
                        "batch_samples": batch_samples,
                        "variants": {"variantcallers": [dd.get_variantcaller(d) for d in gsamples],
                                     "calls": [d.get("vrn_file") for d in gsamples]}})
            out.append(cur)

    def by_original_order(d):
        return min([sample_order.index(s) for s in d["batch_samples"] if s in sample_order])
    return sorted(out, key=by_original_order)