Retrieve the variant region BED file from a population of items.

    If tumor/normal, return the tumor BED file. If a population, return
    the BED file covering the most bases.
10599
def population_variant_regions(items, merged=False):
    
    def _get_variant_regions(data):
        out = dd.get_variant_regions(data) or dd.get_sample_callable(data)
        # Only need to merge for variant region inputs, not callable BED regions which don't overlap
        if merged and dd.get_variant_regions(data):
            merged_out = dd.get_variant_regions_merged(data)
            if merged_out:
                out = merged_out
            else:
                out = merge_overlaps(out, data)
        return out
    import pybedtools
    if len(items) == 1:
        return _get_variant_regions(items[0])
    else:
        paired = vcfutils.get_paired(items)
        if paired:
            return _get_variant_regions(paired.tumor_data)
        else:
            vrs = []
            for data in items:
                vr_bed = _get_variant_regions(data)
                if vr_bed:
                    vrs.append((pybedtools.BedTool(vr_bed).total_coverage(), vr_bed))
            vrs.sort(reverse=True)
            if vrs:
                return vrs[0][1]
428434
def run(align_bams, items, ref_file, assoc_files, region=None, out_file=None):
    
    paired = vcfutils.get_paired_bams(align_bams, items)
    assert paired and not paired.normal_bam, ("smCounter2 supports tumor-only variant calling: %s" %
                                              (",".join([dd.get_sample_name(d) for d in items])))
    vrs = bedutils.population_variant_regions(items)
    target = shared.subset_variant_regions(vrs, region,
                                            out_file, items=items, do_merge=True)
    out_file = out_file.replace(".vcf.gz", ".vcf")
    out_prefix = utils.splitext_plus(os.path.basename(out_file))[0]
    if not utils.file_exists(out_file) and not utils.file_exists(out_file + ".gz"):
        with file_transaction(paired.tumor_data, out_file) as tx_out_file:
            cmd = ["smCounter2", "--runPath", os.path.dirname(tx_out_file),
                   "--outPrefix", out_prefix,
                   "--bedTarget", target, "--refGenome", ref_file,
                   "--bamFile", paired.tumor_bam, "--bamType", "consensus",
                   "--nCPU", dd.get_num_cores(paired.tumor_data)]
            do.run(cmd, "smcounter2 variant calling")
            for fname in glob.glob(os.path.join(os.path.dirname(tx_out_file), "*.smCounter*")):
                shutil.move(fname, os.path.join(os.path.dirname(out_file), os.path.basename(fname)))
            utils.symlink_plus(os.path.join(os.path.dirname(out_file),
                                            "%s.smCounter.cut.vcf" % out_prefix),
                               out_file)
    return vcfutils.bgzip_and_index(out_file, paired.tumor_data["config"], remove_orig=False,
                                    prep_cmd="sed 's#FORMAT\t%s#FORMAT\t%s#' | %s" %
                                    (out_prefix, dd.get_sample_name(paired.tumor_data),
                                     vcfutils.add_contig_to_header_cl(dd.get_ref_file(paired.tumor_data), out_file)))
660918
private Variant getGenomicVariant(Variant variant) {
        if (variant.getAlternate().equals(".")) {  // reference positions are not variants
            return null;
        } else {
            String ref;
            // large deletion
            if (variant.getAlternate().equals("<DEL>")) {
                // .get("_") because studyId and fileId are empty strings when VariantSource is initialized at readInputFile
                int end = Integer.valueOf(variant.getSourceEntries().get("_").getAttributes().get("END"));
                ref = StringUtils.repeat("N", end - variant.getStart());
                return new Variant(variant.getChromosome(), variant.getStart(),
                        ref, variant.getAlternate().equals("") ? "-" : variant.getAlternate());
                // TODO: structural variants are not yet properly handled. Implement and remove this patch asap
            } else if (variant.getAlternate().startsWith("<")
                    || (variant.getAlternate().length() > 1 && variant.getReference().length() > 1)) {
                return null;
            } else {
                ref = variant.getReference().equals("") ? "-" : variant.getReference();
                return new Variant(variant.getChromosome(), variant.getStart(),
                        ref, variant.getAlternate().equals("") ? "-" : variant.getAlternate());
            }
        }
    }
1025750
def prepare_exclude_file(items, base_file, chrom=None):
    
    items = shared.add_highdepth_genome_exclusion(items)
    out_file = "%s-exclude%s.bed" % (utils.splitext_plus(base_file)[0], "-%s" % chrom if chrom else "")
    if not utils.file_exists(out_file) and not utils.file_exists(out_file + ".gz"):
        with shared.bedtools_tmpdir(items[0]):
            with file_transaction(items[0], out_file) as tx_out_file:
                # Get a bedtool for the full region if no variant regions
                want_bedtool = callable.get_ref_bedtool(tz.get_in(["reference", "fasta", "base"], items[0]),
                                                        items[0]["config"], chrom)
                want_bedtool = pybedtools.BedTool(shared.subset_variant_regions(want_bedtool.saveas().fn,
                                                                                chrom, tx_out_file, items))
                sv_exclude_bed = _get_sv_exclude_file(items)
                if sv_exclude_bed and len(want_bedtool) > 0:
                    want_bedtool = want_bedtool.subtract(sv_exclude_bed, nonamecheck=True).saveas()
                full_bedtool = callable.get_ref_bedtool(tz.get_in(["reference", "fasta", "base"], items[0]),
                                                        items[0]["config"])
                if len(want_bedtool) > 0:
                    full_bedtool.subtract(want_bedtool, nonamecheck=True).saveas(tx_out_file)
                else:
                    full_bedtool.saveas(tx_out_file)
    return out_file
1326332
def subset_bam_by_region(in_file, region, config, out_file_base=None):
    
    if out_file_base is not None:
        base, ext = os.path.splitext(out_file_base)
    else:
        base, ext = os.path.splitext(in_file)
    out_file = "%s-subset%s%s" % (base, region, ext)
    if not file_exists(out_file):
        with pysam.Samfile(in_file, "rb") as in_bam:
            target_tid = in_bam.gettid(region)
            assert region is not None, \
                   "Did not find reference region %s in %s" % \
                   (region, in_file)
            with file_transaction(config, out_file) as tx_out_file:
                with pysam.Samfile(tx_out_file, "wb", template=in_bam) as out_bam:
                    for read in in_bam:
                        if read.tid == target_tid:
                            out_bam.write(read)
    return out_file
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
225673
def has_variant_regions(items, base_file, chrom=None):
    
    if chrom:
        all_vrs = _get_variant_regions(items)
        if len(all_vrs) > 0:
            test = shared.subset_variant_regions(tz.first(all_vrs), chrom, base_file, items)
            if test == chrom:
                return False
    return True
857452
def median_vaf_purity(row, cohort, **kwargs):
    
    patient_id = row["patient_id"]
    patient = cohort.patient_from_id(patient_id)
    variants = cohort.load_variants(patients=[patient], filter_fn=no_filter)
    if patient_id in variants.keys():
        variants = variants[patient_id]
    else:
        return np.nan
    def grab_vaf(variant):
        filterable_variant = FilterableVariant(variant, variants, patient)
        return variant_stats_from_variant(variant, filterable_variant.variant_metadata).tumor_stats.variant_allele_frequency
    vafs = [grab_vaf(variant) for variant in variants]
    return 2 * pd.Series(vafs).median()
1086237
def sample_callable_bed(bam_file, ref_file, data):
    
    from bcbio.heterogeneity import chromhacks
    CovInfo = collections.namedtuple("CovInfo", "callable, raw_callable, depth_files")
    noalt_calling = "noalt_calling" in dd.get_tools_on(data) or "altcontigs" in dd.get_exclude_regions(data)
    def callable_chrom_filter(r):
        """Filter to callable region, potentially limiting by chromosomes.
        """
        return r.name == "CALLABLE" and (not noalt_calling or chromhacks.is_nonalt(r.chrom))
    out_file = "%s-callable_sample.bed" % os.path.splitext(bam_file)[0]
    with shared.bedtools_tmpdir(data):
        sv_bed = regions.get_sv_bed(data)
        callable_bed, depth_files = coverage.calculate(bam_file, data, sv_bed)
        input_regions_bed = dd.get_variant_regions(data)
        if not utils.file_uptodate(out_file, callable_bed):
            with file_transaction(data, out_file) as tx_out_file:
                callable_regions = pybedtools.BedTool(callable_bed)
                filter_regions = callable_regions.filter(callable_chrom_filter)
                if input_regions_bed:
                    if not utils.file_uptodate(out_file, input_regions_bed):
                        input_regions = pybedtools.BedTool(input_regions_bed)
                        filter_regions.intersect(input_regions, nonamecheck=True).saveas(tx_out_file)
                else:
                    filter_regions.saveas(tx_out_file)
    return CovInfo(out_file, callable_bed, depth_files)
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
805533
def _get_sv_exclude_file(items):
    
    sv_bed = utils.get_in(items[0], ("genome_resources", "variation", "sv_repeat"))
    if sv_bed and os.path.exists(sv_bed):
        return sv_bed
171599
def minimize(bed_file):
    
    if not bed_file:
        return bed_file
    else:
        sorted_bed = bt.BedTool(bed_file).cut(range(3)).sort()
        if not sorted_bed.fn.endswith(".bed"):
            return sorted_bed.moveto(sorted_bed.fn + ".bed")
        else:
            return sorted_bed
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}
518077
def _do_prioritize(items):
    
    if not any("tumoronly-prioritization" in dd.get_tools_off(d) for d in items):
        if vcfutils.get_paired_phenotype(items[0]):
            has_tumor = False
            has_normal = False
            for sub_data in items:
                if vcfutils.get_paired_phenotype(sub_data) == "tumor":
                    has_tumor = True
                elif vcfutils.get_paired_phenotype(sub_data) == "normal":
                    has_normal = True
            return has_tumor and not has_normal
629415
def get_split_regions(bed_file, data):
    
    out_file = "%s-analysis_blocks.bed" % utils.splitext_plus(bed_file)[0]
    with shared.bedtools_tmpdir(data):
        if not utils.file_uptodate(out_file, bed_file):
            ref_regions = get_ref_bedtool(dd.get_ref_file(data), data["config"])
            nblock_regions = ref_regions.subtract(pybedtools.BedTool(bed_file)).saveas()
            min_n_size = int(tz.get_in(["config", "algorithm", "nomap_split_size"], data, 250))
            block_filter = NBlockRegionPicker(ref_regions, data["config"], min_n_size)
            final_nblock_regions = nblock_regions.filter(
                block_filter.include_block).saveas().each(block_filter.expand_block).saveas()
            with file_transaction(data, out_file) as tx_out_file:
                final_regions = ref_regions.subtract(final_nblock_regions, nonamecheck=True).\
                                saveas().merge(d=min_n_size).saveas(tx_out_file)
        chroms = set([])
        with shared.bedtools_tmpdir(data):
            for r in pybedtools.BedTool(bed_file):
                chroms.add(r.chrom)
        out = []
        for r in pybedtools.BedTool(out_file):
            if r.chrom in chroms:
                out.append((r.chrom, r.start, r.stop))
        return out
820509
def _call_variants(example_dir, region_bed, data, out_file):
    
    tf_out_file = "%s-tfrecord.gz" % utils.splitext_plus(out_file)[0]
    if not utils.file_exists(tf_out_file):
        with file_transaction(data, tf_out_file) as tx_out_file:
            model = "wes" if strelka2.coverage_interval_from_bed(region_bed) == "targeted" else "wgs"
            cmd = ["dv_call_variants.py", "--cores", dd.get_num_cores(data),
                   "--outfile", tx_out_file, "--examples", example_dir,
                   "--sample", dd.get_sample_name(data), "--model", model]
            do.run(cmd, "DeepVariant call_variants %s" % dd.get_sample_name(data))
    return tf_out_file
619722
def detect_regions(bam_in, bed_file, out_dir, prefix):
    
    bed_file = _reorder_columns(bed_file)
    counts_reads_cmd = ("coverageBed -s -counts -b {bam_in} "
                        "-a {bed_file} | sort -k4,4 "
                        "> {out_dir}/loci.cov")
    # with tx_tmpdir() as temp_dir:
    with utils.chdir(out_dir):
        run(counts_reads_cmd.format(min_trimmed_read_len=min_trimmed_read_len, max_trimmed_read_len=max_trimmed_read_len, **locals()), "Run counts_reads")
        loci_file = _fix_score_column(op.join(out_dir, "loci.cov"))
        return loci_file