Run THetA analysis given output from CNV caller on a tumor/normal pair.
162027
def run(vrn_info, cnvs_by_name, somatic_info):
    
    cmd = _get_cmd("RunTHeTA.py")
    if not cmd:
        logger.info("THetA scripts not found in current PATH. Skipping.")
    else:
        from bcbio.structural import cnvkit
        work_dir = _sv_workdir(somatic_info.tumor_data)
        assert "cnvkit" in cnvs_by_name, "THetA requires CNVkit calls"
        cnv_info = cnvkit.export_theta(cnvs_by_name["cnvkit"], somatic_info.tumor_data)
        cnv_info["theta_input"] = subset_by_supported(cnv_info["theta_input"], _theta_to_coords,
                                                      cnvs_by_name, work_dir, somatic_info.tumor_data)
        return _run_theta(cnv_info, somatic_info.tumor_data, work_dir, run_n3=False)
814380
def _ready_for_het_analysis(items):
    
    paired = vcfutils.get_paired_bams([dd.get_align_bam(d) for d in items], items)
    has_het = any(dd.get_hetcaller(d) for d in items)
    if has_het and paired:
        return get_variants(paired.tumor_data) and _get_calls(paired.tumor_data, cnv_only=True)
1789
def run(self, args, pipeline_command):
        

        # output that must be moved but not renamed
        consistentNaming = ['alignments/normal_dna_fix_pg_sorted.bam',
                           'alignments/normal_dna_fix_pg_sorted.bam.bai',
                           'alignments/rna_genome_sorted.bam',
                           'alignments/rna_genome_sorted.bam.bai',
                           'alignments/rna_transcriptome.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam.bai',
                           'mutations/merged/all_merged.vcf',
                           'rankboost/mhcii_rankboost_concise_results.tsv',
                           'rankboost/mhci_rankboost_concise_results.tsv',
                           ]

        # output that must be renamed as well as moved
        # map of the original name to the final name
        renamingNeeded = {'binding_predictions': 'binding_predictions.tar',
                         'expression': 'expression.tar',
                         'haplotyping': 'haplotyping.tar',
                         'peptides': 'peptides.tar',
                         'rankboost': 'rankboost.tar',
                         'reports': 'reports.tar',
                         'mutations/snpeffed/mutations.vcf': 'all_snpeffed.vcf',
                         'mutations/transgened/mutations.vcf': 'all_transgened.vcf',
                         'mutations/merged': 'merged_perchrom.tar',
                         'mutations/muse': 'muse_perchrom.tar',
                         'mutations/mutect': 'mutect_perchrom.tar',
                         'mutations/radia': 'radia_perchrom.tar',
                         'mutations/somaticsniper': 'somaticsniper_perchrom.tar',
                         'mutations/strelka/snv': 'strelka_snv_perchrom.tar',
                         'mutations/strelka/indel': 'strelka_indel_perchrom.tar'}

        def make_output(output_dir, source_dir):
            """
            :param output_dir: dir to write the output to
            :param source_dir: dir containing the directory structure to be parsed
            :return:
            """
            def make_tar(dir, tar):
                with tarfile.open(tar, "w:gz") as tar:
                    tar.add(dir)
            # the output dir is where the real output directories are written
            protect_outputs = os.listdir(source_dir)
            for protectOut in protect_outputs:
                def getName(fileName):
                    return os.path.join(os.path.join(source_dir, protectOut), fileName)
                # move individual files out
                for fileName in consistentNaming:
                    shutil.copyfile(getName(fileName), os.path.join(output_dir, os.path.basename(fileName)))
                for src, dst in renamingNeeded.iteritems():
                    if dst.endswith('.tar'):
                        make_tar(getName(src), os.path.join(output_dir, dst))
                    else:
                        shutil.copyfile(getName(src), os.path.join(output_dir, dst))
            shutil.rmtree(source_dir)


        # prepare workdir
        mount = self._prepare_mount(args)
        self._workdir = os.path.join(mount, 'Toil-' + self._name)
        # insure the pairs are in the same directory, as protect expects
        # This is made more complicated by the fact CWLTool mounts inputs into random, read-only dirs
        # to get around this we copy all inputs into their own directories that we own
        tumor_dna_dir = os.path.expanduser('~/tumorDNA')
        tumor_rna_dir = os.path.expanduser('~/tumorRNA')
        normal_dna_dir = os.path.expanduser('~/normalDNA')
        os.mkdir(tumor_dna_dir)
        os.mkdir(tumor_rna_dir)
        os.mkdir(normal_dna_dir)
        shutil.copy(args.tumor_dna, tumor_dna_dir)
        shutil.copy(args.tumor_rna, tumor_rna_dir)
        shutil.copy(args.normal_dna, normal_dna_dir)
        shutil.copy(args.tumor_dna2, tumor_dna_dir)
        shutil.copy(args.tumor_rna2, tumor_rna_dir)
        shutil.copy(args.normal_dna2, normal_dna_dir)
        args.tumor_dna = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna))
        args.tumor_dna2 = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna2))
        args.tumor_rna = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna))
        args.tumor_rna2 = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna2))
        args.normal_dna = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna))
        args.normal_dna2 = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna2))

        # prepare config
        args_dict = vars(args)
        args_dict['output_dir'] = mount
        self._config = textwrap.dedent(self._config.format(**args_dict))
        self._sample_name = args_dict["sample_name"]
        config_path = os.path.join(self._workdir, 'config')
        command = self._make_prefix(os.path.join(self._workdir, 'jobStore'),
                                          config_path,
                                          self._workdir) + pipeline_command
        if self._resume and args.resume:
            command.append('--restart')
        self._create_workdir(args)
        with open(config_path, 'w') as f:
            f.write(self._config)

        try:
            subprocess.check_call(command)
        except subprocess.CalledProcessError as e:
            print(e, file=sys.stderr)
        finally:
            log.info('Pipeline terminated, changing ownership of output files from root to user.')
            stat = os.stat(self._mount)
            subprocess.check_call(['chown', '-R', '{}:{}'.format(stat.st_uid, stat.st_gid),
                                   self._mount])
            make_output(self._mount, os.path.join(self._mount, 'output'))
            if self._no_clean and args.no_clean:
                log.info('Flag "--no-clean" was used, therefore %s was not deleted.', self._workdir)
            else:
                log.info('Cleaning up temporary directory: %s', self._workdir)
                shutil.rmtree(self._workdir)
589967
def heterogzygote_counts(paired):
    
    work_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(paired.tumor_data), "structural", "counts"))
    key = "germline_het_pon"
    het_bed = tz.get_in(["genome_resources", "variation", key], paired.tumor_data)
    vr = bedutils.population_variant_regions([x for x in [paired.tumor_data, paired.normal_data] if x])
    cur_het_bed = bedutils.intersect_two(het_bed, vr, work_dir, paired.tumor_data)
    tumor_counts = _run_collect_allelic_counts(cur_het_bed, key, work_dir, paired.tumor_data)
    normal_counts = (_run_collect_allelic_counts(cur_het_bed, key, work_dir, paired.normal_data)
                     if paired.normal_data else None)
    if normal_counts:
        tumor_counts, normal_counts = _filter_by_normal(tumor_counts, normal_counts, paired.tumor_data)
    return tumor_counts, normal_counts
255278
def run_vardict(align_bams, items, ref_file, assoc_files, region=None,
                out_file=None):
    
    items = shared.add_highdepth_genome_exclusion(items)
    if vcfutils.is_paired_analysis(align_bams, items):
        call_file = _run_vardict_paired(align_bams, items, ref_file,
                                        assoc_files, region, out_file)
    else:
        vcfutils.check_paired_problems(items)
        call_file = _run_vardict_caller(align_bams, items, ref_file,
                                        assoc_files, region, out_file)
    return call_file
1001618
def get_analysis_intervals(data, vrn_file, base_dir):
    
    from bcbio.bam import callable
    if vrn_file and vcfutils.is_gvcf_file(vrn_file):
        callable_bed = _callable_from_gvcf(data, vrn_file, base_dir)
        if callable_bed:
            return callable_bed

    if data.get("ensemble_bed"):
        return data["ensemble_bed"]
    elif dd.get_sample_callable(data):
        return dd.get_sample_callable(data)
    elif data.get("align_bam"):
        return callable.sample_callable_bed(data["align_bam"], dd.get_ref_file(data), data)[0]
    elif data.get("work_bam"):
        return callable.sample_callable_bed(data["work_bam"], dd.get_ref_file(data), data)[0]
    elif data.get("work_bam_callable"):
        data = utils.deepish_copy(data)
        data["work_bam"] = data.pop("work_bam_callable")
        return callable.sample_callable_bed(data["work_bam"], dd.get_ref_file(data), data)[0]
    elif tz.get_in(["config", "algorithm", "callable_regions"], data):
        return tz.get_in(["config", "algorithm", "callable_regions"], data)
    elif tz.get_in(["config", "algorithm", "variant_regions"], data):
        return tz.get_in(["config", "algorithm", "variant_regions"], data)
367946
def run_analysis(self, argv):
        """"""
        args = self._parser.parse_args(argv)

        if not HAVE_ST:
            raise RuntimeError(
                "Trying to run fermipy analysis, but don't have ST")

        workdir = os.path.dirname(args.config)
        _config_file = self._clone_config_and_srcmaps(args.config, args.seed)

        gta = GTAnalysis(_config_file, logging={'verbosity': 3},
                         fileio={'workdir_regex': '\.xml$|\.npy$'})
        gta.load_roi(args.roi_baseline)

        simfile = os.path.join(workdir, 'sim_%s_%s.yaml' %
                               (args.sim, args.sim_profile))

        mcube_file = "%s_%s_%06i" % (args.sim, args.sim_profile, args.seed)
        sim_config = utils.load_yaml(simfile)

        injected_source = sim_config.get('injected_source', None)
        if injected_source is not None:
            src_dict =  injected_source['source_model']
            src_dict['ra'] = gta.config['selection']['ra']
            src_dict['dec'] = gta.config['selection']['dec']
            injected_name = injected_source['name']
            gta.add_source(injected_name, src_dict)
            gta.write_model_map(mcube_file)
            mc_spec_dict = dict(true_counts=gta.model_counts_spectrum(injected_name),
                                energies=gta.energies,
                                model=src_dict)
            mcspec_file = os.path.join(workdir,
                                       "mcspec_%s_%06i.yaml" % (mcube_file, args.seed))
            utils.write_yaml(mc_spec_dict, mcspec_file)
        else:
            injected_name = None

        gta.write_roi('sim_baseline_%06i' % args.seed)

        test_sources = []
        for profile in args.profiles:
            profile_path = os.path.join(workdir, 'profile_%s.yaml' % profile)
            test_source = load_yaml(profile_path)
            test_sources.append(test_source)
            first = args.seed
            last = first + args.nsims
            for seed in range(first, last):
                self._run_simulation(gta, args.roi_baseline,
                                     injected_name, test_sources, first, seed,
                                     non_null_src=args.non_null_src)
1000316
def _calculate_sv_coverage_cnvkit(data, work_dir):
    
    from bcbio.variation import coverage
    from bcbio.structural import annotate
    out_target_file = os.path.join(work_dir, "%s-target-coverage.cnn" % dd.get_sample_name(data))
    out_anti_file = os.path.join(work_dir, "%s-antitarget-coverage.cnn" % dd.get_sample_name(data))
    if ((not utils.file_exists(out_target_file) or not utils.file_exists(out_anti_file)) and
          (dd.get_align_bam(data) or dd.get_work_bam(data))):
        target_cov = coverage.run_mosdepth(data, "target", tz.get_in(["regions", "bins", "target"], data))
        anti_cov = coverage.run_mosdepth(data, "antitarget", tz.get_in(["regions", "bins", "antitarget"], data))
        target_cov_genes = annotate.add_genes(target_cov.regions, data, max_distance=0)
        out_target_file = _add_log2_depth(target_cov_genes, out_target_file, data)
        out_anti_file = _add_log2_depth(anti_cov.regions, out_anti_file, data)
    return out_target_file, out_anti_file
435075
def _check_lods(parts, tumor_thresh, normal_thresh, indexes):
    
    try:
        gl_index = parts[8].split(":").index("GL")
    except ValueError:
        return True
    try:
        tumor_gls = [float(x) for x in parts[indexes["tumor"]].strip().split(":")[gl_index].split(",") if x != "."]
        if tumor_gls:
            tumor_lod = max(tumor_gls[i] - tumor_gls[0] for i in range(1, len(tumor_gls)))
        else:
            tumor_lod = -1.0
    # No GL information, no tumor call (so fail it)
    except IndexError:
        tumor_lod = -1.0
    try:
        normal_gls = [float(x) for x in parts[indexes["normal"]].strip().split(":")[gl_index].split(",") if x != "."]
        if normal_gls:
            normal_lod = min(normal_gls[0] - normal_gls[i] for i in range(1, len(normal_gls)))
        else:
            normal_lod = normal_thresh
    # No GL inofmration, no normal call (so pass it)
    except IndexError:
        normal_lod = normal_thresh
    return normal_lod >= normal_thresh and tumor_lod >= tumor_thresh
518387
def _cnvkit_segment(cnr_file, cov_interval, data, items, out_file=None, detailed=False):
    
    if not out_file:
        out_file = "%s.cns" % os.path.splitext(cnr_file)[0]
    if not utils.file_uptodate(out_file, cnr_file):
        with file_transaction(data, out_file) as tx_out_file:
            if not _cna_has_values(cnr_file):
                with open(tx_out_file, "w") as out_handle:
                    out_handle.write("chromosome\tstart\tend\tgene\tlog2\tprobes\tCN1\tCN2\tbaf\tweight\n")
            else:
                # Scale cores to avoid memory issues with segmentation
                # https://github.com/etal/cnvkit/issues/346
                if cov_interval == "genome":
                    cores = max(1, dd.get_cores(data) // 2)
                else:
                    cores = dd.get_cores(data)
                cmd = [_get_cmd(), "segment", "-p", str(cores), "-o", tx_out_file, cnr_file]
                small_vrn_files = _compatible_small_variants(data, items)
                if len(small_vrn_files) > 0 and _cna_has_values(cnr_file) and cov_interval != "genome":
                    cmd += ["--vcf", small_vrn_files[0].name, "--sample-id", small_vrn_files[0].sample]
                    if small_vrn_files[0].normal:
                        cmd += ["--normal-id", small_vrn_files[0].normal]
                resources = config_utils.get_resources("cnvkit_segment", data["config"])
                user_options = resources.get("options", [])
                cmd += [str(x) for x in user_options]
                if cov_interval == "genome" and "--threshold" not in user_options:
                    cmd += ["--threshold", "0.00001"]
                # For tumors, remove very low normalized regions, avoiding upcaptured noise
                # https://github.com/bcbio/bcbio-nextgen/issues/2171#issuecomment-348333650
                # unless we want detailed segmentation for downstream tools
                paired = vcfutils.get_paired(items)
                if paired:
                    #if detailed:
                    #    cmd += ["-m", "hmm-tumor"]
                    if "--drop-low-coverage" not in user_options:
                        cmd += ["--drop-low-coverage"]
                # preferentially use conda installed Rscript
                export_cmd = ("%s && export TMPDIR=%s && "
                              % (utils.get_R_exports(), os.path.dirname(tx_out_file)))
                do.run(export_cmd + " ".join(cmd), "CNVkit segment")
    return out_file
1225301
def run_analysis(self, argv):
        """"""
        args = self._parser.parse_args(argv)

        name_keys = dict(target_type=args.ttype,
                         target_name=args.target,
                         sim_name=args.sim,
                         fullpath=True)

        orig_dir = NAME_FACTORY.targetdir(**name_keys)
        dest_dir = NAME_FACTORY.sim_targetdir(**name_keys)
        self.copy_target_dir(orig_dir, dest_dir,
                             args.roi_baseline, args.extracopy)
677102
def scan(self, seqs, nreport=100, scan_rc=True, normalize=False):
        

        if not self.threshold:
            sys.stderr.write(
                "Using default threshold of 0.95. "
                "This is likely not optimal!\n"
                )
            self.set_threshold(threshold=0.95)

        seqs = as_fasta(seqs, genome=self.genome)
           
        it = self._scan_sequences(seqs.seqs, 
                    nreport, scan_rc)
       
        if normalize:
            if len(self.meanstd) == 0:
                self.set_meanstd()
            mean_std = [self.meanstd.get(m_id) for m_id in self.motif_ids]
            means = [x[0] for x in  mean_std]
            stds = [x[1] for x in  mean_std]


        for result in it:
            if normalize:
                zresult = [] 
                for i,mrow in enumerate(result):
                    mrow = [((x[0] - means[i]) / stds[i], x[1], x[2]) for x in mrow]
                    zresult.append(mrow)
                yield zresult
            else:
                yield result
692744
def _run_bcbio_variation(vrn_file, rm_file, rm_interval_file, base_dir, sample, caller, data):
    
    val_config_file = _create_validate_config_file(vrn_file, rm_file, rm_interval_file,
                                                   base_dir, data)
    work_dir = os.path.join(base_dir, "work")
    out = {"summary": os.path.join(work_dir, "validate-summary.csv"),
           "grading": os.path.join(work_dir, "validate-grading.yaml"),
           "discordant": os.path.join(work_dir, "%s-eval-ref-discordance-annotate.vcf" % sample)}
    if not utils.file_exists(out["discordant"]) or not utils.file_exists(out["grading"]):
        bcbio_variation_comparison(val_config_file, base_dir, data)
    out["concordant"] = filter(os.path.exists,
                                [os.path.join(work_dir, "%s-%s-concordance.vcf" % (sample, x))
                                 for x in ["eval-ref", "ref-eval"]])[0]
    return out
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
1356696
def _tumor_normal_stats(rec, somatic_info, vcf_rec):
    
    out = {"normal": {"alt": None, "depth": None, "freq": None},
           "tumor": {"alt": 0, "depth": 0, "freq": None}}
    if hasattr(vcf_rec, "samples"):
        samples = [(s, {}) for s in vcf_rec.samples]
        for fkey in ["AD", "AO", "RO", "AF", "DP"]:
            try:
                for i, v in enumerate(rec.format(fkey)):
                    samples[i][1][fkey] = v
            except KeyError:
                pass
    # Handle INFO only inputs
    elif len(rec.samples) == 0:
        samples = [(somatic_info.tumor_name, None)]
    else:
        samples = rec.samples.items()
    for name, sample in samples:
        alt, depth, freq = sample_alt_and_depth(rec, sample)
        if depth is not None and freq is not None:
            if name == somatic_info.normal_name:
                key = "normal"
            elif name == somatic_info.tumor_name:
                key = "tumor"
            out[key]["freq"] = freq
            out[key]["depth"] = depth
            out[key]["alt"] = alt
    return out
566947
def _af_filter(data, in_file, out_file):
    
    min_freq = float(utils.get_in(data["config"], ("algorithm", "min_allele_fraction"), 10)) / 100.0
    logger.debug("Filtering MuTect2 calls with allele fraction threshold of %s" % min_freq)
    ungz_out_file = "%s.vcf" % utils.splitext_plus(out_file)[0]
    if not utils.file_exists(ungz_out_file) and not utils.file_exists(ungz_out_file + ".gz"):
        with file_transaction(data, ungz_out_file) as tx_out_file:
            vcf = cyvcf2.VCF(in_file)
            vcf.add_filter_to_header({
                'ID': 'MinAF',
                'Description': 'Allele frequency is lower than %s%% ' % (min_freq*100) + (
                    '(configured in bcbio as min_allele_fraction)'
                    if utils.get_in(data["config"], ("algorithm", "min_allele_fraction"))
                    else '(default threshold in bcbio; override with min_allele_fraction in the algorithm section)')})
            w = cyvcf2.Writer(tx_out_file, vcf)
            # GATK 3.x can produce VCFs without sample names for empty VCFs
            try:
                tumor_index = vcf.samples.index(dd.get_sample_name(data))
            except ValueError:
                tumor_index = None
            for rec in vcf:
                if tumor_index is not None and np.all(rec.format('AF')[tumor_index] < min_freq):
                    vcfutils.cyvcf_add_filter(rec, 'MinAF')
                w.write_record(rec)
            w.close()
    return vcfutils.bgzip_and_index(ungz_out_file, data["config"])
535424
def run(vrn_info, cnvs_by_name, somatic_info):
    
    config = {"sample_size": 5000}
    work_dir = _cur_workdir(somatic_info.tumor_data)
    if "battenberg" not in cnvs_by_name:
        logger.warn("PhyloWGS requires Battenberg CNV calls, skipping %s"
                    % dd.get_sample_name(somatic_info.tumor_data))
    else:
        ssm_file, cnv_file = _prep_inputs(vrn_info, cnvs_by_name["battenberg"], somatic_info, work_dir, config)
        evolve_file = _run_evolve(ssm_file, cnv_file, work_dir, somatic_info.tumor_data)
        summary_file = _prepare_summary(evolve_file, ssm_file, cnv_file, work_dir, somatic_info)
        print(summary_file)