Determine if we should perform prioritization.

    Currently done on tumor-only input samples and feeding into PureCN
    which needs the germline annotations.
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
1281403
def _get_purecn_files(paired, work_dir, require_exist=False):
    
    out_base = os.path.join(work_dir, "%s-purecn" % (dd.get_sample_name(paired.tumor_data)))
    out = {"plot": {}}
    all_files = []
    for plot in ["chromosomes", "local_optima", "segmentation", "summary"]:
        if plot == "summary":
            cur_file = "%s.pdf" % out_base
        else:
            cur_file = "%s_%s.pdf" % (out_base, plot)
        if not require_exist or os.path.exists(cur_file):
            out["plot"][plot] = cur_file
            all_files.append(os.path.basename(cur_file))
    for key, ext in [["hetsummary", ".csv"], ["dnacopy", "_dnacopy.seg"], ["genes", "_genes.csv"],
                     ["log", ".log"], ["loh", "_loh.csv"], ["rds", ".rds"],
                     ["variants", "_variants.csv"]]:
        cur_file = "%s%s" % (out_base, ext)
        if not require_exist or os.path.exists(cur_file):
            out[key] = cur_file
            all_files.append(os.path.basename(cur_file))
    return out_base, out, all_files
1167150
def _init_objaartall(self):
        """"""
        kws = {
            'sortgo':lambda nt: [nt.NS, nt.dcnt],
            # fmtgo=('{p_fdr_bh:8.2e} {GO} '
            # Formatting for GO terms in grouped GO list
            'fmtgo':('{hdr1usr01:2} {NS} {GO} {s_fdr_bh:8} '
                     '{dcnt:5} {childcnt:3} R{reldepth:02} '
                     '{D1:5} {GO_name} ({study_count} study genes)\n'),
            # Formatting for GO terms listed under each gene
            'fmtgo2':('{hdr1usr01:2} {NS} {GO} {s_fdr_bh:8} '
                      '{dcnt:5} R{reldepth:02} '
                      '{GO_name} ({study_count} study genes)\n'),
            # itemid2name=ensmusg2symbol}
            }
        return AArtGeneProductSetsAll(self.grprdflt, self.hdrobj, **kws)
960489
def prep_seq2c_bed(data):
    
    if dd.get_background_cnv_reference(data, "seq2c"):
        bed_file = _background_to_bed(dd.get_background_cnv_reference(data, "seq2c"), data)
    else:
        bed_file = regions.get_sv_bed(data)
    if bed_file:
        bed_file = bedutils.clean_file(bed_file, data, prefix="svregions-")
    else:
        bed_file = bedutils.clean_file(dd.get_variant_regions(data), data)
    if not bed_file:
        return None

    col_num = bt.BedTool(bed_file).field_count()
    if col_num < 4:
        annotated_file = annotate.add_genes(bed_file, data, max_distance=0)
        if annotated_file == bed_file:
            raise ValueError("BED file for Seq2C must be annotated with gene names, "
                             "however the input BED is 3-columns and we have no transcript "
                             "data to annotate with " + bed_file)
        annotated_file = annotate.gene_one_per_line(annotated_file, data)
    else:
        annotated_file = bed_file

    ready_file = "%s-seq2cclean.bed" % (utils.splitext_plus(annotated_file)[0])
    if not utils.file_uptodate(ready_file, annotated_file):
        bed = bt.BedTool(annotated_file)
        if col_num > 4 and col_num != 8:
            bed = bed.cut(range(4))
        bed = bed.filter(lambda x: x.name not in ["", ".", "-"])
        with file_transaction(data, ready_file) as tx_out_file:
            bed.saveas(tx_out_file)
        logger.debug("Saved Seq2C clean annotated ready input BED into " + ready_file)

    return ready_file
1016037
def max_normal_germline_depth(in_file, params, somatic_info):
    
    bcf_in = pysam.VariantFile(in_file)
    depths = []
    for rec in bcf_in:
        stats = _is_possible_loh(rec, bcf_in, params, somatic_info)
        if tz.get_in(["normal", "depth"], stats):
            depths.append(tz.get_in(["normal", "depth"], stats))
    if depths:
        return np.median(depths) * NORMAL_FILTER_PARAMS["max_depth_percent"]
667731
private void handleSequenceAnnotation(String seqName, String featureName, String value) {
		if (featureName.equals(GS_ACCESSION_NUMBER)) {
			stockholmStructure.addGSAccessionNumber(seqName, value);
		} else if (featureName.equals(GS_DESCRIPTION)) {
			stockholmStructure.addGSDescription(seqName, value);
		} else if (featureName.equals(GS_DATABASE_REFERENCE)) {
			stockholmStructure.addGSdbReference(seqName, value);
		} else if (featureName.equals(GS_ORGANISM_SPECIES)) {
			stockholmStructure.addGSOrganismSpecies(seqName, value);
		} else if (featureName.equals(GS_ORGANISM_CLASSIFICATION)) {
			stockholmStructure.addGSOrganismClassification(seqName, value);
		} else if (featureName.equals(GS_LOOK)) {
			stockholmStructure.addGSLook(seqName, value);
		} else {
			// unknown feature
			logger.warn("Unknown Sequence Feature [{}].\nPlease contact the Biojava team.", featureName);
		}
	}
1158816
def extract(data, items, out_dir=None):
    
    if vcfutils.get_paired_phenotype(data):
        if len(items) == 1:
            germline_vcf = _remove_prioritization(data["vrn_file"], data, out_dir)
            germline_vcf = vcfutils.bgzip_and_index(germline_vcf, data["config"])
            data["vrn_file_plus"] = {"germline": germline_vcf}
    return data
257282
protected ImmutableList<Variable> prioritizeVariables(IntermediateQuery query, JoinLikeNode joinLikeNode) {

        /*
         * Sequential (order matters)
         */
        return extractAncestors(query, joinLikeNode).stream()
                .sequential()
                .flatMap(p -> extractPriorityVariables(query, p.parent, p.position))
                .distinct()
                .collect(ImmutableCollectors.toList());
    }
339699
def top_priority_effect_per_gene_id(self):
        """"""
        return OrderedDict(
            (gene_id, top_priority_effect(variant_effects))
            for (gene_id, variant_effects)
            in self.groupby_gene_id().items())
670565
def create_dexseq_annotation(gff, count_file):
    
    out_file = count_file + ".ann"
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out:
        with open(tx_out, 'w') as out_handle:
            with open(gff) as in_handle:
                for line in in_handle:
                    cols = line.strip().split("\t")
                    if cols[2] == "exonic_part":
                        exon = [f for f in cols[8].split(";") if f.strip().startswith("exonic_part_number")]
                        gene = [f for f in cols[8].split(";") if f.strip().startswith("gene_id")]
                        exon = exon[0].replace("\"", "").split()[1]
                        gene = gene[0].replace("\"", "").split()[1]
                        length = int(cols[4]) - int(cols[3]) + 1
                        line = "%s:%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (gene, exon, gene,
                                                                    cols[0], cols[3],
                                                                    cols[4],
                                                                    length, cols[6])
                        out_handle.write(line)
374164
def _has_population_germline(rec):
    
    for k in population_keys:
        if k in rec.header.info:
            return True
    return False
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)))
258465
def _run_purecn(paired, work_dir):
    
    segfns = {"cnvkit": _segment_normalized_cnvkit, "gatk-cnv": _segment_normalized_gatk}
    out_base, out, all_files = _get_purecn_files(paired, work_dir)
    failed_file = out_base + "-failed.log"
    cnr_file = tz.get_in(["depth", "bins", "normalized"], paired.tumor_data)
    if not utils.file_uptodate(out["rds"], cnr_file) and not utils.file_exists(failed_file):
        cnr_file, seg_file = segfns[cnvkit.bin_approach(paired.tumor_data)](cnr_file, work_dir, paired)
        from bcbio import heterogeneity
        vcf_file = heterogeneity.get_variants(paired.tumor_data, include_germline=False)[0]["vrn_file"]
        vcf_file = germline.filter_to_pass_and_reject(vcf_file, paired, out_dir=work_dir)
        with file_transaction(paired.tumor_data, out_base) as tx_out_base:
            # Use UCSC style naming for human builds to support BSgenome
            genome = ("hg19" if dd.get_genome_build(paired.tumor_data) in ["GRCh37", "hg19"]
                      else dd.get_genome_build(paired.tumor_data))
            cmd = ["PureCN.R", "--seed", "42", "--out", tx_out_base, "--rds", "%s.rds" % tx_out_base,
                   "--sampleid", dd.get_sample_name(paired.tumor_data),
                   "--genome", genome,
                   "--vcf", vcf_file, "--tumor", cnr_file,
                   "--segfile", seg_file, "--funsegmentation", "Hclust", "--maxnonclonal", "0.3"]
            if dd.get_num_cores(paired.tumor_data) > 1:
                cmd += ["--cores", str(dd.get_num_cores(paired.tumor_data))]
            try:
                cmd = "export R_LIBS_USER=%s && %s && %s" % (utils.R_sitelib(), utils.get_R_exports(),
                                                             " ".join([str(x) for x in cmd]))
                do.run(cmd, "PureCN copy number calling")
            except subprocess.CalledProcessError as msg:
                if _allowed_errors(str(msg)):
                    logger.info("PureCN failed to find solution for %s: skipping" %
                                dd.get_sample_name(paired.tumor_data))
                    with open(failed_file, "w") as out_handle:
                        out_handle.write(str(msg))
                else:
                    logger.exception()
                    raise
            for f in all_files:
                if os.path.exists(os.path.join(os.path.dirname(tx_out_base), f)):
                    shutil.move(os.path.join(os.path.dirname(tx_out_base), f),
                                os.path.join(os.path.dirname(out_base), f))
    out = _get_purecn_files(paired, work_dir, require_exist=True)[1]
    return out if (out.get("rds") and os.path.exists(out["rds"])) else None
1233630
def _is_possible_loh(rec, vcf_rec, params, somatic_info, use_status=False, max_normal_depth=None):
    
    if _is_biallelic_snp(rec) and _passes_plus_germline(rec, use_status=use_status):
        stats = _tumor_normal_stats(rec, somatic_info, vcf_rec)
        depths = [tz.get_in([x, "depth"], stats) for x in ["normal", "tumor"]]
        depths = [d for d in depths if d is not None]
        normal_freq = tz.get_in(["normal", "freq"], stats)
        tumor_freq = tz.get_in(["tumor", "freq"], stats)
        if all([d > params["min_depth"] for d in depths]):
            if max_normal_depth and tz.get_in(["normal", "depth"], stats, 0) > max_normal_depth:
                return None
            if normal_freq is not None:
                if normal_freq >= params["min_freq"] and normal_freq <= params["max_freq"]:
                    return stats
            elif (tumor_freq >= params["tumor_only"]["min_freq"] and
                    tumor_freq <= params["tumor_only"]["max_freq"]):
                if (vcf_rec and not _has_population_germline(vcf_rec)) or is_population_germline(rec):
                    return stats
116658
def best_precursor(clus, loci):
    
    data_loci = sort_precursor(clus, loci)
    current_size = data_loci[0][5]
    best = 0
    for item, locus in enumerate(data_loci):
        if locus[3] - locus[2] > 70:
            if locus[5] > current_size * 0.8:
                best = item
                break
    best_loci = data_loci[best]
    del data_loci[best]
    data_loci.insert(0, best_loci)
    return data_loci
1045158
def split_variants_by_sample(data):
    
    # not split, do nothing
    if "group_orig" not in data:
        return [[data]]
    # cancer tumor/normal
    elif (vcfutils.get_paired_phenotype(data)
            and "tumor" in [vcfutils.get_paired_phenotype(d) for d in get_orig_items(data)]):
        out = []
        for i, sub_data in enumerate(get_orig_items(data)):
            if vcfutils.get_paired_phenotype(sub_data) == "tumor":
                cur_batch = tz.get_in(["metadata", "batch"], data)
                if cur_batch:
                    sub_data["metadata"]["batch"] = cur_batch
                sub_data["vrn_file"] = data["vrn_file"]
            else:
                sub_data.pop("vrn_file", None)
            out.append([sub_data])
        return out
    # joint calling or population runs, do not split back up and keep in batches
    else:
        out = []
        for sub_data in get_orig_items(data):
            cur_batch = tz.get_in(["metadata", "batch"], data)
            if cur_batch:
                sub_data["metadata"]["batch"] = cur_batch
            sub_data["vrn_file_batch"] = data["vrn_file"]
            sub_data["vrn_file"] = data["vrn_file"]
            out.append([sub_data])
        return out
1011812
def _run_cnvkit_cancer(items, background):
    
    paired = vcfutils.get_paired_bams([x["align_bam"] for x in items], items)
    normal_data = [x for x in items if dd.get_sample_name(x) != paired.tumor_name]
    tumor_ready, normal_ready = _match_batches(paired.tumor_data, normal_data[0] if normal_data else None)
    ckouts = _run_cnvkit_shared([tumor_ready], [normal_ready] if normal_ready else [])
    if not ckouts:
        return items
    assert len(ckouts) == 1
    tumor_data = _associate_cnvkit_out(ckouts, [paired.tumor_data], is_somatic=True)
    return tumor_data + normal_data