Return DeepVariant calling on germline samples.

    region can be a single region or list of multiple regions for multicore calling.
766663
def run(align_bams, items, ref_file, assoc_files, region, out_file):
    
    assert not vcfutils.is_paired_analysis(align_bams, items), \
        ("DeepVariant currently only supports germline calling: %s" %
         (", ".join([dd.get_sample_name(d) for d in items])))
    assert len(items) == 1, \
        ("DeepVariant currently only supports single sample calling: %s" %
         (", ".join([dd.get_sample_name(d) for d in items])))
    out_file = _run_germline(align_bams[0], items[0], ref_file,
                             region, out_file)
    return vcfutils.bgzip_and_index(out_file, items[0]["config"])
1065473
def _get_calls(data, cnv_only=False):
    
    cnvs_supported = set(["cnvkit", "battenberg"])
    out = {}
    for sv in data.get("sv", []):
        if not cnv_only or sv["variantcaller"] in cnvs_supported:
            out[sv["variantcaller"]] = sv
    return out
496628
private Map<String, VariantSet> getVariants(String defaultSkin, Set<String> skinNames, String defaultLocaleName,
			Set<String> localeVariants) {

		Map<String, VariantSet> skinVariants = new HashMap<>();
		if (!skinNames.isEmpty()) {
			skinVariants.put(JawrConstant.SKIN_VARIANT_TYPE,
					new VariantSet(JawrConstant.SKIN_VARIANT_TYPE, defaultSkin, skinNames));
		}
		if (!localeVariants.isEmpty()) {

			skinVariants.put(JawrConstant.LOCALE_VARIANT_TYPE,
					new VariantSet(JawrConstant.LOCALE_VARIANT_TYPE, defaultLocaleName, localeVariants));
		}

		return skinVariants;
	}
1098679
def _get_variant_regions(items):
    
    return list(filter(lambda x: x is not None,
                       [tz.get_in(("config", "algorithm", "variant_regions"), data)
                        for data in items
                        if tz.get_in(["config", "algorithm", "coverage_interval"], data) != "genome"]))
902027
def variantcall_batch_region(items):
    
    items = [utils.to_single_data(x) for x in items]
    align_bams = [dd.get_align_bam(x) for x in items]
    variantcaller = _get_batch_variantcaller(items)
    region_blocks = list(set([tuple(x.get("region_block")) for x in items if "region_block" in x]))
    assert len(region_blocks) == 1, region_blocks
    region_block = region_blocks[0]
    # Pre-called input variant files
    if not variantcaller and all(d.get("vrn_file") for d in items):
        return {"vrn_file_region": None, "region_block": region_block}
    caller_fn = get_variantcallers()[variantcaller]
    assoc_files = tz.get_in(("genome_resources", "variation"), items[0], {})
    region = _region_to_coords(region_block[0])
    chrom, start, end = region
    region_str = "_".join(str(x) for x in region)
    batch_name = _get_batch_name(items)
    out_file = os.path.join(dd.get_work_dir(items[0]), variantcaller, chrom,
                            "%s-%s-block.vcf.gz" % (batch_name, region_str))
    utils.safe_makedir(os.path.dirname(out_file))
    with pshared.bedtools_tmpdir(items[0]):
        if variantcaller in SUPPORT_MULTICORE:
            call_file = caller_fn(align_bams, items, dd.get_ref_file(items[0]), assoc_files,
                                  [_region_to_coords(r) for r in region_block], out_file)
        else:
            call_file = _run_variantcall_batch_multicore(items, region_block, out_file)
    return {"vrn_file_region": call_file, "region_block": region_block}
66493
def use_vqsr(algs, call_file=None):
    
    from bcbio.variation import vcfutils
    vqsr_callers = set(["gatk", "gatk-haplotype"])
    vqsr_sample_thresh = 50
    vqsr_supported = collections.defaultdict(int)
    coverage_intervals = set([])
    for alg in algs:
        callers = alg.get("variantcaller")
        if isinstance(callers, six.string_types):
            callers = [callers]
        if not callers:  # no variant calling, no VQSR
            continue
        if "vqsr" in (alg.get("tools_off") or []):  # VQSR turned off
            continue
        for c in callers:
            if c in vqsr_callers:
                if "vqsr" in (alg.get("tools_on") or []):  # VQSR turned on:
                    vqsr_supported[c] += 1
                    coverage_intervals.add("genome")
                # Do not try VQSR for gVCF inputs
                elif call_file and vcfutils.is_gvcf_file(call_file):
                    pass
                else:
                    coverage_intervals.add(alg.get("coverage_interval", "exome").lower())
                    vqsr_supported[c] += 1
    if len(vqsr_supported) > 0:
        num_samples = max(vqsr_supported.values())
        if "genome" in coverage_intervals or num_samples >= vqsr_sample_thresh:
            return True
    return False
1302920
def __getSequenceVariants(self, x1, polyStart, polyStop, listSequence) :
		""""""
		if polyStart < len(self.polymorphisms) and polyStart < polyStop: 
			sequence = copy.copy(listSequence)
			ret = []
			
			pk = self.polymorphisms[polyStart]
			posInSequence = pk[0]-x1
			if posInSequence < len(listSequence) : 
				for allele in pk[1] :
					sequence[posInSequence] = allele
					ret.extend(self.__getSequenceVariants(x1, polyStart+1, polyStop, sequence))
			
			return ret
		else :
			return [''.join(listSequence)]
1327429
def get_variants_in_region(self, chrom, start, end):
        """"""
        if self.chrom is not None and chrom == self.chrom:
            # We are going to search for 'NA' since the chromosome was set
            chrom = "NA"

        iterator = self._bgen.iter_variants_in_region(
            CHROM_STR_DECODE.get(chrom, chrom), start, end,
        )
        for info, dosage in iterator:
            yield Genotypes(
                Variant(
                    info.name, CHROM_STR_ENCODE.get(info.chrom, info.chrom),
                    info.pos, [info.a1, info.a2],
                ),
                dosage,
                reference=info.a1,
                coded=info.a2,
                multiallelic=True,
            )
550494
def haplotype_caller(align_bams, items, ref_file, assoc_files,
                       region=None, out_file=None):
    
    if out_file is None:
        out_file = "%s-variants.vcf.gz" % utils.splitext_plus(align_bams[0])[0]
    if not utils.file_exists(out_file):
        num_cores = dd.get_num_cores(items[0])
        broad_runner, params = \
            _shared_gatk_call_prep(align_bams, items, ref_file, region, out_file, num_cores)
        gatk_type = broad_runner.gatk_type()
        assert gatk_type in ["restricted", "gatk4"], \
            "Require full version of GATK 2.4+, or GATK4 for haplotype calling"
        with file_transaction(items[0], out_file) as tx_out_file:
            resources = config_utils.get_resources("gatk-spark", items[0]["config"])
            spark_opts = [str(x) for x in resources.get("options", [])]
            if _use_spark(num_cores, gatk_type, items, spark_opts):
                params += ["-T", "HaplotypeCallerSpark"]
                if spark_opts:
                    params += spark_opts
                else:
                    params += ["--spark-master", "local[%s]" % num_cores,
                               "--conf", "spark.local.dir=%s" % os.path.dirname(tx_out_file),
                               "--conf", "spark.driver.host=localhost", "--conf", "spark.network.timeout=800",
                               "--conf", "spark.executor.heartbeatInterval=100"]
            else:
                params += ["-T", "HaplotypeCaller"]
            params += ["--annotation", "ClippingRankSumTest",
                       "--annotation", "DepthPerSampleHC"]
            # Enable hardware based optimizations in GATK 3.1+
            if LooseVersion(broad_runner.gatk_major_version()) >= LooseVersion("3.1"):
                if _supports_avx():
                    # Scale down HMM thread default to avoid overuse of cores
                    # https://github.com/bcbio/bcbio-nextgen/issues/2442
                    if gatk_type == "gatk4":
                        params += ["--native-pair-hmm-threads", "1"]
                    # GATK4 selects the right HMM optimization automatically with FASTEST_AVAILABLE
                    # GATK3 needs to be explicitly set
                    else:
                        params += ["--pair_hmm_implementation", "VECTOR_LOGLESS_CACHING"]
            resources = config_utils.get_resources("gatk-haplotype", items[0]["config"])
            if "options" in resources:
                params += [str(x) for x in resources.get("options", [])]
            # Prepare gVCFs if doing joint calling
            is_joint = False
            if _joint_calling(items) or any("gvcf" in dd.get_tools_on(d) for d in items):
                is_joint = True
                # If joint calling parameters not set in user options
                if not any([x in ["--emit-ref-confidence", "-ERC", "--emitRefConfidence"] for x in params]):
                    if gatk_type == "gatk4":
                        params += ["--emit-ref-confidence", "GVCF"]
                    else:
                        params += ["--emitRefConfidence", "GVCF"]
                        params += ["--variant_index_type", "LINEAR", "--variant_index_parameter", "128000"]
                # Set GQ banding to not be single GQ resolution
                # No recommended default but try to balance resolution and size
                # http://gatkforums.broadinstitute.org/gatk/discussion/7051/recommendation-best-practices-gvcf-gq-bands

                if not any([x in ["-GQB"] for x in params]):
                    for boundary in [10, 20, 30, 40, 60, 80]:
                        params += ["-GQB", str(boundary)]
            # Enable non-diploid calling in GATK 3.3+
            if LooseVersion(broad_runner.gatk_major_version()) >= LooseVersion("3.3"):
                params += ["-ploidy", str(ploidy.get_ploidy(items, region))]
            if gatk_type == "gatk4":
                # GATK4 Spark calling does not support bgzipped output, use plain VCFs
                if is_joint and _use_spark(num_cores, gatk_type, items, spark_opts):
                    tx_out_file = tx_out_file.replace(".vcf.gz", ".vcf")
                params += ["--output", tx_out_file]
            else:
                params += ["-o", tx_out_file]
            broad_runner.new_resources("gatk-haplotype")
            memscale = {"magnitude": 0.9 * num_cores, "direction": "increase"} if num_cores > 1 else None
            try:
                broad_runner.run_gatk(params, os.path.dirname(tx_out_file), memscale=memscale,
                                      parallel_gc=_use_spark(num_cores, gatk_type, items, spark_opts))
            except subprocess.CalledProcessError as msg:
                # Spark failing on regions without any reads, write an empty VCF instead
                # https://github.com/broadinstitute/gatk/issues/4234
                if (_use_spark(num_cores, gatk_type, items, spark_opts) and
                      str(msg).find("java.lang.UnsupportedOperationException: empty collection") >= 0 and
                      str(msg).find("at org.apache.spark.rdd.RDD") >= 0):
                    vcfutils.write_empty_vcf(tx_out_file, samples=[dd.get_sample_name(d) for d in items])
                else:
                    raise
            if tx_out_file.endswith(".vcf"):
                vcfutils.bgzip_and_index(tx_out_file, items[0]["config"])


    # avoid bug in GATK where files can get output as non-compressed
    if out_file.endswith(".gz") and not os.path.exists(out_file + ".tbi"):
        with open(out_file, "r") as in_handle:
            is_plain_text = in_handle.readline().startswith("##fileformat")
        if is_plain_text:
            text_out_file = out_file
            out_file = out_file.replace(".vcf.gz", ".vcf")
            shutil.move(text_out_file, out_file)
    return vcfutils.bgzip_and_index(out_file, items[0]["config"])
1123205
def _run_genotype_gvcfs_gatk3(data, region, vrn_files, ref_file, out_file):
    
    if not utils.file_exists(out_file):
        broad_runner = broad.runner_from_config(data["config"])
        with file_transaction(data, out_file) as tx_out_file:
            assoc_files = tz.get_in(("genome_resources", "variation"), data, {})
            if not assoc_files: assoc_files = {}
            params = ["-T", "GenotypeGVCFs",
                      "-R", ref_file, "-o", tx_out_file,
                      "-L", bamprep.region_to_gatk(region),
                      "--max_alternate_alleles", "4"]
            for vrn_file in vrn_files:
                params += ["--variant", vrn_file]
            if assoc_files.get("dbsnp"):
                params += ["--dbsnp", assoc_files["dbsnp"]]
            broad_runner.new_resources("gatk-haplotype")
            cores = dd.get_cores(data)
            if cores > 1:
                # GATK performs poorly with memory usage when parallelizing
                # with a large number of cores but makes use of extra memory,
                # so we cap at 6 cores.
                # See issue #1565 for discussion
                # Recent GATK 3.x versions also have race conditions with multiple
                # threads, so limit to 1 and keep memory available
                # https://gatkforums.broadinstitute.org/wdl/discussion/8718/concurrentmodificationexception-in-gatk-3-7-genotypegvcfs
                # params += ["-nt", str(min(6, cores))]
                memscale = {"magnitude": 0.9 * cores, "direction": "increase"}
            else:
                memscale = None
            broad_runner.run_gatk(params, memscale=memscale, parallel_gc=True)
    return vcfutils.bgzip_and_index(out_file, data["config"])
645239
private Map<String, String> getVariantMap(ResourceBundlesHandler bundlesHandler, String contextPath,
			String bundlePath) {

		JoinableResourceBundle bundle = bundlesHandler.resolveBundleForPath(bundlePath);
		Set<String> variantTypes = bundle.getVariants().keySet();

		String variantKey = getVariantKey(contextPath);
		String[] variantValues = new String[0];
		if (variantKey.length() > 0) {
			if (variantKey.length() == 1) {
				variantValues = new String[] { "", "" };
			} else {
				variantValues = variantKey.split(String.valueOf(JawrConstant.VARIANT_SEPARATOR_CHAR));
			}
		}

		Map<String, String> variants = new HashMap<>();
		if (variantTypes.size() != variantValues.length) {
			throw new BundlingProcessException("For the resource '" + contextPath
					+ "', the number variant types for the bundle don't match the variant values.");
		}
		int i = 0;
		for (String variantType : variantTypes) {
			variants.put(variantType, variantValues[i++]);
		}
		return variants;
	}
374067
def square_batch_region(data, region, bam_files, vrn_files, out_file):
    
    from bcbio.variation import sentieon, strelka2
    if not utils.file_exists(out_file):
        jointcaller = tz.get_in(("config", "algorithm", "jointcaller"), data)
        if jointcaller in ["%s-joint" % x for x in SUPPORTED["general"]]:
            _square_batch_bcbio_variation(data, region, bam_files, vrn_files, out_file, "square")
        elif jointcaller in ["%s-merge" % x for x in SUPPORTED["general"]]:
            _square_batch_bcbio_variation(data, region, bam_files, vrn_files, out_file, "merge")
        elif jointcaller in ["%s-joint" % x for x in SUPPORTED["gatk"]]:
            gatkjoint.run_region(data, region, vrn_files, out_file)
        elif jointcaller in ["%s-joint" % x for x in SUPPORTED["gvcf"]]:
            strelka2.run_gvcfgenotyper(data, region, vrn_files, out_file)
        elif jointcaller in ["%s-joint" % x for x in SUPPORTED["sentieon"]]:
            sentieon.run_gvcftyper(vrn_files, out_file, region, data)
        else:
            raise ValueError("Unexpected joint calling approach: %s." % jointcaller)
    if region:
        data["region"] = region
    data = _fix_orig_vcf_refs(data)
    data["vrn_file"] = out_file
    return [data]
309815
def _collapse_by_bam_variantcaller(samples):
    
    by_bam = collections.OrderedDict()
    for data in (x[0] for x in samples):
        work_bam = utils.get_in(data, ("combine", "work_bam", "out"), data.get("align_bam"))
        variantcaller = get_variantcaller(data)
        if isinstance(work_bam, list):
            work_bam = tuple(work_bam)
        key = (multi.get_batch_for_key(data), work_bam, variantcaller)
        try:
            by_bam[key].append(data)
        except KeyError:
            by_bam[key] = [data]
    out = []
    for grouped_data in by_bam.values():
        cur = grouped_data[0]
        cur.pop("region", None)
        region_bams = cur.pop("region_bams", None)
        if region_bams and len(region_bams[0]) > 1:
            cur.pop("work_bam", None)
        out.append([cur])
    return out
1087975
private String getAvailableVariant(Map<String, String> curVariants) {

		String variantKey = null;
		if (variants != null) {
			Map<String, String> availableVariants = generatorRegistry.getAvailableVariantMap(variants, curVariants);
			variantKey = VariantUtils.getVariantKey(availableVariants);
		}

		return variantKey;
	}
605506
def _run_cortex_on_region(region, align_bam, ref_file, work_dir, out_file_base, config):
    
    kmers = [31, 51, 71]
    min_reads = 1750
    cortex_dir = config_utils.get_program("cortex", config, "dir")
    stampy_dir = config_utils.get_program("stampy", config, "dir")
    vcftools_dir = config_utils.get_program("vcftools", config, "dir")
    if cortex_dir is None or stampy_dir is None:
        raise ValueError("cortex_var requires path to pre-built cortex and stampy")
    region_str = "{0}-{1}-{2}".format(*region)
    base_dir = safe_makedir(os.path.join(work_dir, region_str))
    try:
        out_vcf_base = os.path.join(base_dir, "{0}-{1}".format(
                    os.path.splitext(os.path.basename(out_file_base))[0], region_str))
        out_file = os.path.join(work_dir, os.path.basename("{0}.vcf".format(out_vcf_base)))
        if not file_exists(out_file):
            fastq = _get_fastq_in_region(region, align_bam, out_vcf_base)
            if _count_fastq_reads(fastq, min_reads) < min_reads:
                vcfutils.write_empty_vcf(out_file)
            else:
                local_ref, genome_size = _get_local_ref(region, ref_file, out_vcf_base)
                indexes = _index_local_ref(local_ref, cortex_dir, stampy_dir, kmers)
                cortex_out = _run_cortex(fastq, indexes, {"kmers": kmers, "genome_size": genome_size,
                                                          "sample": get_sample_name(align_bam)},
                                         out_vcf_base, {"cortex": cortex_dir, "stampy": stampy_dir,
                                                        "vcftools": vcftools_dir},
                                         config)
                if cortex_out:
                    _remap_cortex_out(cortex_out, region, out_file)
                else:
                    vcfutils.write_empty_vcf(out_file)
    finally:
        if os.path.exists(base_dir):
            shutil.rmtree(base_dir)
    return [region[0], int(region[1]), int(region[2]), out_file]
356759
def get_genomic_seq_for_region(self, chrom, start_pos, end_pos):
        
        
        headers = {"content-type": "text/plain"}
        
        self.attempt = 0
        ext = "/sequence/region/human/{}:{}..{}:1".format(chrom, start_pos, end_pos)
        
        return self.ensembl_request(ext, headers)
973044
def _create_genome_regions(data):
    
    work_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "coverage", dd.get_sample_name(data)))
    variant_regions = os.path.join(work_dir, "target-genome.bed")
    with file_transaction(data, variant_regions) as tx_variant_regions:
        with open(tx_variant_regions, "w") as out_handle:
            for c in shared.get_noalt_contigs(data):
                out_handle.write("%s\t%s\t%s\n" % (c.name, 0, c.size))
    return variant_regions