Produce joint called variants from input gVCF files.
498737
def run_gvcftyper(vrn_files, out_file, region, data):
    
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            license = license_export(data)
            ref_file = dd.get_ref_file(data)
            input_files = " ".join(vrn_files)
            region = bamprep.region_to_gatk(region)
            cmd = ("{license}sentieon driver -r {ref_file} --interval {region} "
                   "--algo GVCFtyper {tx_out_file} {input_files}")
            do.run(cmd.format(**locals()), "Sentieon GVCFtyper")
    return out_file
207636
def summarize_vc(items):
    
    items = [utils.to_single_data(x) for x in utils.flatten(items)]
    items = [_normalize_vc_input(x) for x in items]
    items = validate.summarize_grading(items)
    items = [utils.to_single_data(x) for x in items]
    out = {"validate": validate.combine_validations(items),
           "variants": {"calls": [], "gvcf": [], "samples": []}}
    added = set([])
    variants_by_sample = collections.defaultdict(list)
    sample_order = []
    for data in items:
        batch_samples = data.get("batch_samples", [dd.get_sample_name(data)])
        for s in batch_samples:
            if s not in sample_order:
                sample_order.append(s)
        if data.get("vrn_file"):
            # Only get batches if we're actually doing variantcalling in bcbio
            # otherwise we'll be using the original files
            names = dd.get_batches(data) if dd.get_variantcaller(data) else None
            if not names:
                names = [dd.get_sample_name(data)]
            batch_name = names[0]
            if data.get("vrn_file_joint") is not None:
                to_add = [("vrn_file", "gvcf", dd.get_sample_name(data)),
                          ("vrn_file_joint", "calls", batch_name)]
            else:
                to_add = [("vrn_file", "calls", batch_name)]
            for vrn_key, out_key, name in to_add:
                cur_name = "%s-%s" % (name, dd.get_variantcaller(data))
                out_file = os.path.join(utils.safe_makedir(os.path.join(dd.get_work_dir(data),
                                                                        "variants", out_key)),
                                        "%s.vcf.gz" % cur_name)
                for s in batch_samples:
                    variants_by_sample[s].append(out_file)
                if cur_name not in added:
                    added.add(cur_name)
                    # Ideally could symlink here but doesn't appear to work with
                    # Docker container runs on Toil where PATHs don't get remapped
                    utils.copy_plus(os.path.realpath(data[vrn_key]), out_file)
                    vcfutils.bgzip_and_index(out_file, data["config"])
                    out["variants"][out_key].append(out_file)
    for sample in sample_order:
        out["variants"]["samples"].append(variants_by_sample[sample])
    return [out]
282948
def run_vep(in_file, data):
    
    if not vcfutils.vcf_has_variants(in_file):
        return None
    out_file = utils.append_stem(in_file, "-vepeffects")
    assert in_file.endswith(".gz") and out_file.endswith(".gz")
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            vep_dir, ensembl_name = prep_vep_cache(data["genome_build"],
                                                   tz.get_in(["reference", "fasta", "base"], data))
            if vep_dir:
                cores = tz.get_in(("config", "algorithm", "num_cores"), data, 1)
                fork_args = ["--fork", str(cores)] if cores > 1 else []
                vep = config_utils.get_program("vep", data["config"])
                # HGVS requires a bgzip compressed, faidx indexed input file or is unusable slow
                if dd.get_ref_file_compressed(data):
                    hgvs_compatible = True
                    config_args = ["--fasta", dd.get_ref_file_compressed(data)]
                else:
                    hgvs_compatible = False
                    config_args = ["--fasta", dd.get_ref_file(data)]
                if vcfanno.is_human(data):
                    plugin_fns = {"loftee": _get_loftee, "maxentscan": _get_maxentscan,"genesplicer": _get_genesplicer,
                                  "spliceregion": _get_spliceregion, "G2P": _get_G2P}
                    plugins = ["loftee", "G2P"]
                    if "vep_splicesite_annotations" in dd.get_tools_on(data):
                        # "genesplicer" too unstable so currently removed
                        plugins += ["maxentscan", "spliceregion"]
                    for plugin in plugins:
                        plugin_args = plugin_fns[plugin](data)
                        config_args += plugin_args
                    config_args += ["--sift", "b", "--polyphen", "b"]
                    if hgvs_compatible:
                        config_args += ["--hgvsg","--hgvs", "--shift_hgvs", "1"]
                if (dd.get_effects_transcripts(data).startswith("canonical")
                      or tz.get_in(("config", "algorithm", "clinical_reporting"), data)):
                    config_args += ["--most_severe"]
                else:
                    config_args += ["--flag_pick_allele_gene"]
                if ensembl_name.endswith("_merged"):
                    config_args += ["--merged"]
                    ensembl_name = ensembl_name.replace("_merged", "")
                resources = config_utils.get_resources("vep", data["config"])
                extra_args = [str(x) for x in resources.get("options", [])]
                cmd = [vep, "--vcf", "-o", "stdout", "-i", in_file] + fork_args + extra_args + \
                      ["--species", ensembl_name,
                       "--no_stats", "--cache",
                        "--offline", "--dir", vep_dir,
                       "--symbol", "--numbers", "--biotype", "--total_length", "--canonical",
                       "--gene_phenotype", "--ccds", "--uniprot", "--domains", "--regulatory",
                       "--protein", "--tsl", "--appris", "--af", "--max_af", "--af_1kg", "--af_esp", "--af_gnomad",
                       "--pubmed", "--variant_class", "--allele_number"] + config_args
                perl_exports = utils.get_perl_exports()
                # Remove empty fields (';;') which can cause parsing errors downstream
                cmd = "%s && %s | sed '/^#/! s/;;/;/g' | bgzip -c > %s" % (perl_exports, " ".join(cmd), tx_out_file)
                do.run(cmd, "Ensembl variant effect predictor", data)
    if utils.file_exists(out_file):
        return vcfutils.bgzip_and_index(out_file, data["config"])
988722
def parallel_variantcall_region(samples, run_parallel):
    
    to_process, extras = _dup_samples_by_variantcaller(samples)
    split_fn = _split_by_ready_regions(".vcf.gz", "work_bam", get_variantcaller)
    samples = _collapse_by_bam_variantcaller(
        grouped_parallel_split_combine(to_process, split_fn,
                                       multi.group_batches, run_parallel,
                                       "variantcall_sample", "concat_variant_files",
                                       "vrn_file", ["region", "sam_ref", "config"]))
    return extras + samples
175408
def _get_variant_file(x, key, suffix="", sample=None, ignore_do_upload=False):
    
    out = []
    fname = utils.get_in(x, key)
    upload_key = list(key)
    upload_key[-1] = "do_upload"
    do_upload = tz.get_in(tuple(upload_key), x, True)
    if fname and (ignore_do_upload or do_upload):
        if fname.endswith(".vcf.gz"):
            out.append({"path": fname,
                        "type": "vcf.gz",
                        "ext": "%s%s" % (x["variantcaller"], suffix),
                        "variantcaller": x["variantcaller"]})
            if utils.file_exists(fname + ".tbi"):
                out.append({"path": fname + ".tbi",
                            "type": "vcf.gz.tbi",
                            "index": True,
                            "ext": "%s%s" % (x["variantcaller"], suffix),
                            "variantcaller": x["variantcaller"]})
        elif fname.endswith((".vcf", ".bed", ".bedpe", ".bedgraph", ".cnr", ".cns", ".cnn", ".txt", ".tsv")):
            ftype = utils.splitext_plus(fname)[-1][1:]
            if ftype == "txt":
                extended_ftype = fname.split("-")[-1]
                if "/" not in extended_ftype:
                    ftype = extended_ftype
            out.append({"path": fname,
                        "type": ftype,
                        "ext": "%s%s" % (x["variantcaller"], suffix),
                        "variantcaller": x["variantcaller"]})
    if sample:
        out_sample = []
        for x in out:
            x["sample"] = sample
            out_sample.append(x)
        return out_sample
    else:
        return out
508087
def _variants_fills(fields, fills, info_types):
    """"""
    if fills is None:
        # no fills specified by user
        fills = dict()
    for f, vcf_type in zip(fields, info_types):
        if f == 'FILTER':
            fills[f] = False
        elif f not in fills:
            if f in config.STANDARD_VARIANT_FIELDS:
                fills[f] = config.DEFAULT_VARIANT_FILL[f]
            else:
                fills[f] = config.DEFAULT_FILL_MAP[vcf_type]
    # convert to tuple for zipping with fields
    fills = tuple(fills[f] for f in fields)
    return fills
375595
def _subset_to_variant_regions(callable_file, variant_regions, data):
    
    out_file = "%s-vrsubset.bed" % utils.splitext_plus(callable_file)[0]
    if not utils.file_uptodate(out_file, callable_file):
        with file_transaction(data, out_file) as tx_out_file:
            with utils.open_gzipsafe(callable_file) as in_handle:
                pybedtools.BedTool(in_handle).intersect(variant_regions).saveas(tx_out_file)
    return out_file
456220
private static List<String> getVariantKeys(String variantKeyPrefix, Collection<String> variants) {

		List<String> variantKeys = new ArrayList<>();
		for (String variant : variants) {
			if (variant == null) {
				variant = "";
			}
			if (variantKeyPrefix == null) {
				variantKeys.add(variant);
			} else {
				variantKeys.add(variantKeyPrefix + variant);
			}
		}

		return variantKeys;
	}
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
414317
private static String serializeVariantSets(Map<String, VariantSet> map) {
		StringBuilder result = new StringBuilder();

		for (Entry<String, VariantSet> entry : map.entrySet()) {
			result.append(entry.getKey()).append(":");
			VariantSet variantSet = (VariantSet) entry.getValue();
			result.append(variantSet.getDefaultVariant()).append(":");
			result.append(getCommaSeparatedString(variantSet));
			result.append(";");
		}

		return result.toString();
	}
661146
def _run_concat_variant_files_bcftools(in_list, out_file, config, naive=False):
    
    if not utils.file_exists(out_file):
        with file_transaction(config, out_file) as tx_out_file:
            bcftools = config_utils.get_program("bcftools", config)
            output_type = "z" if out_file.endswith(".gz") else "v"
            if naive:
                args = "--naive"
            else:
                args = "--allow-overlaps"
            cmd = "{bcftools} concat {args} -O {output_type} --file-list {in_list} -o {tx_out_file}"
            do.run(cmd.format(**locals()), "bcftools concat variants")
    if out_file.endswith(".gz"):
        bgzip_and_index(out_file, config)
    return out_file
9577
def _get_files_variantcall(sample):
    
    out = []
    algorithm = sample["config"]["algorithm"]
    out = _maybe_add_summary(algorithm, sample, out)
    out = _maybe_add_alignment(algorithm, sample, out)
    out = _maybe_add_callable(sample, out)
    out = _maybe_add_disambiguate(algorithm, sample, out)
    out = _maybe_add_variant_file(algorithm, sample, out)
    out = _maybe_add_sv(algorithm, sample, out)
    out = _maybe_add_hla(algorithm, sample, out)
    out = _maybe_add_heterogeneity(algorithm, sample, out)

    out = _maybe_add_validate(algorithm, sample, out)
    return _add_meta(out, sample)
1322353
def cyvcf2(context, vcf, include, exclude, chrom, start, end, loglevel, silent,
           individual, no_inds):
    """"""
    coloredlogs.install(log_level=loglevel)
    start_parsing = datetime.now()
    log.info("Running cyvcf2 version %s", __version__)

    if include and exclude:
        log.warning("Can not use include and exclude at the same time")
        context.abort()

    region = ''
    if (chrom or start or end):
        if not (chrom and start and end):
            log.warning("Please specify chromosome, start and end for region")
            context.abort()
        else:
            region = "{0}:{1}-{2}".format(chrom, start, end)

    vcf_obj = VCF(vcf)
    
    for inclusion in include:
        if vcf_obj.contains(inclusion):
            log.info("Including %s in output", inclusion)
        else:
            log.warning("%s does not exist in header", inclusion)
            context.abort()

    for exclusion in exclude:
        if vcf_obj.contains(exclusion):
            log.info("Excluding %s in output", exclusion)
        else:
            log.warning("%s does not exist in header", exclusion)
            context.abort()
    
    if individual:
        # Check if the choosen individuals exists in vcf
        test = True
        for ind_id in individual:
            if ind_id not in vcf_obj.samples:
                log.warning("Individual '%s' does not exist in vcf", ind_id)
                test = False
            if not test:
                context.abort()
        # Convert individuals to list for VCF.set_individuals
        individual = list(individual)
    else:
        individual = None
    
    # Set individual to be empty list to skip all genotypes
    if no_inds:
        individual = []

    if not silent:
        print_header(vcf_obj, include, exclude, individual)

    nr_variants = None
    try:
        for nr_variants, variant in enumerate(vcf_obj(region)):
            if not silent:
                print_variant(variant, include, exclude)
    except Exception as err:
        log.warning(err)
        context.abort()
    
    if nr_variants is None:
        log.info("No variants in vcf")
        return
    
    log.info("{0} variants parsed".format(nr_variants+1))
    log.info("Time to parse variants: {0}".format(datetime.now() - start_parsing))
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
1149192
def annotate_svs(adapter, vcf_obj):
    
    for nr_variants, variant in enumerate(vcf_obj, 1):
        variant_info = get_coords(variant)
        match = adapter.get_structural_variant(variant_info)
        if match:
            annotate_variant(variant, match)
        yield variant
409381
def generate_simple_vcf(filename, variant_collection):
    
    contigs = []
    positions = []
    refs = []
    alts = []
    for variant in variant_collection:
        contigs.append("chr" + variant.contig)
        positions.append(variant.start)
        refs.append(variant.ref)
        alts.append(variant.alt)
    df = pd.DataFrame()
    df["contig"] = contigs
    df["position"] = positions
    df["id"] = ["."] * len(variant_collection)
    df["ref"] = refs
    df["alt"] = alts
    df["qual"] = ["."] * len(variant_collection)
    df["filter"] = ["."] * len(variant_collection)
    df["info"] = ["."] * len(variant_collection)
    df["format"] = ["GT:AD:DP"] * len(variant_collection)
    normal_ref_depths = [randint(1, 10) for v in variant_collection]
    normal_alt_depths = [randint(1, 10) for v in variant_collection]
    df["n1"] = ["0:%d,%d:%d" % (normal_ref_depths[i], normal_alt_depths[i],
                                normal_ref_depths[i] + normal_alt_depths[i])
                for i in range(len(variant_collection))]
    tumor_ref_depths = [randint(1, 10) for v in variant_collection]
    tumor_alt_depths = [randint(1, 10) for v in variant_collection]
    df["t1"] = ["0/1:%d,%d:%d" % (tumor_ref_depths[i], tumor_alt_depths[i], tumor_ref_depths[i] + tumor_alt_depths[i])
                for i in range(len(variant_collection))]

    with open(filename, "w") as f:
        f.write("##fileformat=VCFv4.1\n")
        f.write("##reference=file:///projects/ngs/resources/gatk/2.3/ucsc.hg19.fasta\n")

    with open(filename, "a") as f:
        df.to_csv(f, sep="\t", index=None, header=None)
253701
def _run_concat_variant_files_gatk4(input_file_list, out_file, config):
    
    if not utils.file_exists(out_file):
        with file_transaction(config, out_file) as tx_out_file:
            params = ["-T", "GatherVcfs", "-I", input_file_list, "-O", tx_out_file]
            # Use GATK4 for merging, tools_off: [gatk4] applies to variant calling
            config = utils.deepish_copy(config)
            if "gatk4" in dd.get_tools_off({"config": config}):
                config["algorithm"]["tools_off"].remove("gatk4")
            # Allow specification of verbosity in the unique style this tool uses
            resources = config_utils.get_resources("gatk", config)
            opts = [str(x) for x in resources.get("options", [])]
            if "--verbosity" in opts:
                params += ["--VERBOSITY:%s" % opts[opts.index("--verbosity") + 1]]
            broad_runner = broad.runner_from_config(config)
            broad_runner.run_gatk(params)
    return out_file