Run Scalpel indel calling, either paired tumor/normal or germline calling.
321493
def run_scalpel(align_bams, items, ref_file, assoc_files, region=None,
                  out_file=None):
    
    if region is None:
        message = ("A region must be provided for Scalpel")
        raise ValueError(message)
    if is_paired_analysis(align_bams, items):
        call_file = _run_scalpel_paired(align_bams, items, ref_file,
                                          assoc_files, region, out_file)
    else:
        call_file = _run_scalpel_caller(align_bams, items, ref_file,
                                          assoc_files, region, out_file)
    return call_file
1360704
def run
      raise unless prediction.is_a?(Query)

      start  = Time.now

      number = '-?\d*\.?\d*'
      match  = @prediction.definition.match(/QI:#{number}\|(#{number})\|
                                             (#{number})\|#{number}\|
                                             #{number}\|#{number}\|#{number}\|
                                             #{number}\|#{number}/x)

      raise NotEnoughEvidence if match.nil?

      # % of splice sites confirmed by EST/mRNA-seq alignment
      splice_sites = (match[1].to_f * 100).round
      # % of exons that match an EST/mRNA-seq alignment
      exons = (match[2].to_f * 100).round

      @validation_report = MakerQIValidationOutput.new(@short_header, @header,
                                                       @description,
                                                       splice_sites, exons)
      @validation_report.run_time = Time.now - start
      @validation_report
    rescue NotEnoughEvidence
      @validation_report = ValidationReport.new('No MAKER Quality Index',
                                                :warning, @short_header,
                                                @header, @description)
    rescue StandardError
      @validation_report = ValidationReport.new('Unexpected error', :error,
                                                @short_header, @header,
                                                @description)
      @validation_report.errors.push 'Unexpected Error'
    end
750060
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, ("Pisces 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)
    min_af = float(dd.get_min_allele_fraction(paired.tumor_data)) / 100.0
    if not utils.file_exists(out_file):
        base_out_name = utils.splitext_plus(os.path.basename(paired.tumor_bam))[0]
        raw_file = "%s.vcf" % utils.splitext_plus(out_file)[0]
        with file_transaction(paired.tumor_data, raw_file) as tx_out_file:
            ref_dir = _prep_genome(os.path.dirname(tx_out_file), paired.tumor_data)
            out_dir = os.path.dirname(tx_out_file)
            cores = dd.get_num_cores(paired.tumor_data)
            emit_min_af = min_af / 10.0
            cmd = ("pisces --bampaths {paired.tumor_bam} --genomepaths {ref_dir} --intervalpaths {target} "
                   "--maxthreads {cores} --minvf {emit_min_af} --vffilter {min_af} "
                   "--ploidy somatic --gvcf false -o {out_dir}")
            # Recommended filtering for low frequency indels
            # https://github.com/bcbio/bcbio-nextgen/commit/49d0cbb1f6dcbea629c63749e2f9813bd06dcee3#commitcomment-29765373
            cmd += " -RMxNFilter 5,9,0.35"
            # For low frequency UMI tagged variants, set higher variant thresholds
            # https://github.com/Illumina/Pisces/issues/14#issuecomment-399756862
            if min_af < (1.0 / 100.0):
                cmd += " --minbasecallquality 30"
            do.run(cmd.format(**locals()), "Pisces tumor-only somatic calling")
            shutil.move(os.path.join(out_dir, "%s.vcf" % base_out_name),
                        tx_out_file)
        vcfutils.bgzip_and_index(raw_file, paired.tumor_data["config"],
                                 prep_cmd="sed 's#%s.bam#%s#' | %s" %
                                 (base_out_name, dd.get_sample_name(paired.tumor_data),
                                  vcfutils.add_contig_to_header_cl(dd.get_ref_file(paired.tumor_data), out_file)))
    return vcfutils.bgzip_and_index(out_file, paired.tumor_data["config"])
316841
def _create_input(paired, out_file, ref_file, snp_file, qsnp_file):
    """"""
    ini_file["[inputFiles]"]["dbSNP"] = snp_file
    ini_file["[inputFiles]"]["ref"] = ref_file
    ini_file["[inputFiles]"]["normalBam"] = paired.normal_bam
    ini_file["[inputFiles]"]["tumourBam"] = paired.tumor_bam
    ini_file["[ids]"]["normalSample"] = paired.normal_name
    ini_file["[ids]"]["tumourSample"] = paired.tumor_name
    ini_file["[ids]"]["donor"] = paired.tumor_name
    ini_file["[outputFiles]"]["vcf"] = out_file
    with open(qsnp_file, "w") as out_handle:
        for k, v in ini_file.items():
            out_handle.write("%s\n" % k)
            for opt, value in v.items():
                if value != "":
                    out_handle.write("%s = %s\n" % (opt, value))
836131
def _run_amber(paired, work_dir, lenient=False):
    
    amber_dir = utils.safe_makedir(os.path.join(work_dir, "amber"))
    out_file = os.path.join(amber_dir, "%s.amber.baf" % dd.get_sample_name(paired.tumor_data))
    if not utils.file_exists(out_file) or not utils.file_exists(out_file + ".pcf"):
        with file_transaction(paired.tumor_data, out_file) as tx_out_file:
            key = "germline_het_pon"
            het_bed = tz.get_in(["genome_resources", "variation", key], paired.tumor_data)
            cmd = ["AMBER"] + _get_jvm_opts(tx_out_file, paired.tumor_data) + \
                  ["-threads", dd.get_num_cores(paired.tumor_data),
                   "-tumor", dd.get_sample_name(paired.tumor_data),
                   "-tumor_bam", dd.get_align_bam(paired.tumor_data),
                   "-reference", dd.get_sample_name(paired.normal_data),
                   "-reference_bam", dd.get_align_bam(paired.normal_data),
                   "-ref_genome", dd.get_ref_file(paired.tumor_data),
                   "-bed", het_bed,
                   "-output_dir", os.path.dirname(tx_out_file)]
            if lenient:
                cmd += ["-max_het_af_percent", "1.0"]
            try:
                do.run(cmd, "PURPLE: AMBER baf generation")
            except subprocess.CalledProcessError as msg:
                if not lenient and _amber_allowed_errors(str(msg)):
                    return _run_amber(paired, work_dir, True)
            for f in os.listdir(os.path.dirname(tx_out_file)):
                if f != os.path.basename(tx_out_file):
                    shutil.move(os.path.join(os.path.dirname(tx_out_file), f),
                                os.path.join(amber_dir, f))
    return out_file
259234
def run(vrn_info, calls_by_name, somatic_info, do_plots=True, handle_failures=True):
    
    if "seq2c" in calls_by_name:
        cnv_info = calls_by_name["seq2c"]
    elif "cnvkit" in calls_by_name:
        cnv_info = calls_by_name["cnvkit"]
    else:
        raise ValueError("BubbleTree only currently support CNVkit and Seq2c: %s" % ", ".join(calls_by_name.keys()))
    work_dir = _cur_workdir(somatic_info.tumor_data)
    class OutWriter:
        def __init__(self, out_handle):
            self.writer = csv.writer(out_handle)
        def write_header(self):
            self.writer.writerow(["chrom", "start", "end", "freq"])
        def write_row(self, rec, stats):
            self.writer.writerow([_to_ucsc_style(rec.chrom), rec.start, rec.stop, stats["tumor"]["freq"]])
    vcf_csv = prep_vrn_file(vrn_info["vrn_file"], vrn_info["variantcaller"],
                            work_dir, somatic_info, OutWriter, cnv_info["cns"])
    cnv_csv = _prep_cnv_file(cnv_info["cns"], cnv_info["variantcaller"], work_dir,
                             somatic_info.tumor_data)
    wide_lrr = cnv_info["variantcaller"] == "cnvkit" and somatic_info.normal_bam is None
    return _run_bubbletree(vcf_csv, cnv_csv, somatic_info.tumor_data, wide_lrr, do_plots,
                           handle_failures)
1365005
def run_rnaseq_variant_calling(data):
    
    variantcaller = dd.get_variantcaller(data)
    if isinstance(variantcaller, list) and len(variantcaller) > 1:
        logger.error("Only one variantcaller can be run for RNA-seq at "
                     "this time. Post an issue here "
                     "(https://github.com/bcbio/bcbio-nextgen/issues) "
                     "if this is something you need to do.")
        sys.exit(1)

    if variantcaller:
        if "gatk-haplotype" in variantcaller:
            data = variation.rnaseq_gatk_variant_calling(data)
        if vardict.get_vardict_command(data):
            data = variation.rnaseq_vardict_variant_calling(data)
        vrn_file = dd.get_vrn_file(data)
    return [[data]]
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
117022
def _run_purple(paired, het_file, depth_file, vrn_files, work_dir):
    
    purple_dir = utils.safe_makedir(os.path.join(work_dir, "purple"))
    out_file = os.path.join(purple_dir, "%s.purple.cnv" % dd.get_sample_name(paired.tumor_data))
    if not utils.file_exists(out_file):
        with file_transaction(paired.tumor_data, out_file) as tx_out_file:
            cmd = ["PURPLE"] + _get_jvm_opts(tx_out_file, paired.tumor_data) + \
                  ["-amber", os.path.dirname(het_file), "-baf", het_file,
                   "-cobalt", os.path.dirname(depth_file),
                   "-gc_profile", dd.get_variation_resources(paired.tumor_data)["gc_profile"],
                   "-output_dir", os.path.dirname(tx_out_file),
                   "-ref_genome", "hg38" if dd.get_genome_build(paired.tumor_data) == "hg38" else "hg19",
                   "-run_dir", work_dir,
                   "-threads", dd.get_num_cores(paired.tumor_data),
                   "-tumor_sample", dd.get_sample_name(paired.tumor_data),
                   "-ref_sample", dd.get_sample_name(paired.normal_data)]
            if vrn_files:
                cmd += ["-somatic_vcf", vrn_files[0]["vrn_file"]]
            # Avoid X11 display errors when writing plots
            cmd = "unset DISPLAY && %s" % " ".join([str(x) for x in cmd])
            do.run(cmd, "PURPLE: purity and ploidy estimation")
            for f in os.listdir(os.path.dirname(tx_out_file)):
                if f != os.path.basename(tx_out_file):
                    shutil.move(os.path.join(os.path.dirname(tx_out_file), f),
                                os.path.join(purple_dir, f))
    out_file_export = os.path.join(purple_dir, "%s-purple-cnv.tsv" % (dd.get_sample_name(paired.tumor_data)))
    if not utils.file_exists(out_file_export):
        utils.symlink_plus(out_file, out_file_export)
    out = {"variantcaller": "purple", "call_file": out_file_export,
           "vrn_file": titancna.to_vcf(out_file_export, "PURPLE", _get_header, _export_to_vcf,
                                       paired.tumor_data),
           "plot": {}, "metrics": {}}
    for name, ext in [("copy_number", "copyNumber"), ("minor_allele", "minor_allele"), ("variant", "variant")]:
        plot_file = os.path.join(purple_dir, "plot", "%s.%s.png" % (dd.get_sample_name(paired.tumor_data), ext))
        if os.path.exists(plot_file):
            out["plot"][name] = plot_file
    purity_file = os.path.join(purple_dir, "%s.purple.purity" % dd.get_sample_name(paired.tumor_data))
    with open(purity_file) as in_handle:
        header = in_handle.readline().replace("#", "").split("\t")
        vals = in_handle.readline().split("\t")
        for h, v in zip(header, vals):
            try:
                v = float(v)
            except ValueError:
                pass
            out["metrics"][h] = v
    return out
367471
def _pick_lead_item(items):
    
    if vcfutils.is_paired_analysis([dd.get_align_bam(x) for x in items], items):
        for data in items:
            if vcfutils.get_paired_phenotype(data) == "tumor":
                return data
        raise ValueError("Did not find tumor sample in paired tumor/normal calling")
    else:
        return items[0]
1018157
def mutect2_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):
        paired = vcfutils.get_paired_bams(align_bams, items)
        broad_runner = broad.runner_from_config(items[0]["config"])
        gatk_type = broad_runner.gatk_type()
        _prep_inputs(align_bams, ref_file, items)
        with file_transaction(items[0], out_file) as tx_out_file:
            params = ["-T", "Mutect2" if gatk_type == "gatk4" else "MuTect2",
                      "--annotation", "ClippingRankSumTest",
                      "--annotation", "DepthPerSampleHC"]
            if gatk_type == "gatk4":
                params += ["--reference", ref_file]
            else:
                params += ["-R", ref_file]
            for a in annotation.get_gatk_annotations(items[0]["config"], include_baseqranksum=False):
                params += ["--annotation", a]
            # Avoid issues with BAM CIGAR reads that GATK doesn't like
            if gatk_type == "gatk4":
                params += ["--read-validation-stringency", "LENIENT"]
            params += _add_tumor_params(paired, items, gatk_type)
            params += _add_region_params(region, out_file, items, gatk_type)
            # Avoid adding dbSNP/Cosmic so they do not get fed to variant filtering algorithm
            # Not yet clear how this helps or hurts in a general case.
            #params += _add_assoc_params(assoc_files)
            resources = config_utils.get_resources("mutect2", items[0]["config"])
            if "options" in resources:
                params += [str(x) for x in resources.get("options", [])]
            assert LooseVersion(broad_runner.gatk_major_version()) >= LooseVersion("3.5"), \
                "Require full version of GATK 3.5+ for mutect2 calling"
            broad_runner.new_resources("mutect2")
            gatk_cmd = broad_runner.cl_gatk(params, os.path.dirname(tx_out_file))
            if gatk_type == "gatk4":
                tx_raw_prefilt_file = "%s-raw%s" % utils.splitext_plus(tx_out_file)
                tx_raw_file = "%s-raw-filt%s" % utils.splitext_plus(tx_out_file)
                filter_cmd = _mutect2_filter(broad_runner, tx_raw_prefilt_file, tx_raw_file, ref_file)
                cmd = "{gatk_cmd} -O {tx_raw_prefilt_file} && {filter_cmd}"
            else:
                tx_raw_file = "%s-raw%s" % utils.splitext_plus(tx_out_file)
                cmd = "{gatk_cmd} > {tx_raw_file}"
            do.run(cmd.format(**locals()), "MuTect2")
            out_file = _af_filter(paired.tumor_data, tx_raw_file, out_file)
    return vcfutils.bgzip_and_index(out_file, items[0]["config"])
228064
def _fix_generic_tn_names(paired):
    
    def run(line):
        parts = line.rstrip("\n\r").split("\t")
        if "TUMOR" in parts:
            parts[parts.index("TUMOR")] = paired.tumor_name
        if "TUMOUR" in parts:
            parts[parts.index("TUMOUR")] = paired.tumor_name
        if "NORMAL" in parts:
            assert paired.normal_name
            parts[parts.index("NORMAL")] = paired.normal_name
        return "\t".join(parts) + "\n"
    return run
451861
def _run_germline(bam_file, data, ref_file, region, out_file):
    
    work_dir = utils.safe_makedir("%s-work" % utils.splitext_plus(out_file)[0])
    region_bed = strelka2.get_region_bed(region, [data], out_file, want_gzip=False)
    example_dir = _make_examples(bam_file, data, ref_file, region_bed, out_file, work_dir)
    if _has_candidate_variants(example_dir):
        tfrecord_file = _call_variants(example_dir, region_bed, data, out_file)
        return _postprocess_variants(tfrecord_file, data, ref_file, out_file)
    else:
        return vcfutils.write_empty_vcf(out_file, data["config"], [dd.get_sample_name(data)])
1192181
def rnaseq_variant_calling(samples, run_parallel):
    
    samples = run_parallel("run_rnaseq_variant_calling", samples)
    variantcaller = dd.get_variantcaller(to_single_data(samples[0]))
    if variantcaller and ("gatk-haplotype" in variantcaller):
        out = []
        for d in joint.square_off(samples, run_parallel):
            out.extend([[to_single_data(xs)] for xs in multi.split_variants_by_sample(to_single_data(d))])
        samples = out
    if variantcaller:
        samples = run_parallel("run_rnaseq_ann_filter", samples)
    if variantcaller and ("gatk-haplotype" in variantcaller):
        out = []
        for data in (to_single_data(xs) for xs in samples):
            if "variants" not in data:
                data["variants"] = []
            data["variants"].append({"variantcaller": "gatk-haplotype", "vcf": data["vrn_file_orig"],
                                     "population": {"vcf": data["vrn_file"]}})
            data["vrn_file"] = data.pop("vrn_file_orig")
            out.append([data])
        samples = out
    return samples
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)
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
1036069
def reporter(self):
        
        logging.info('Preparing reports')
        # Populate self.plusdict in order to reuse parsing code from an assembly-based method
        for sample in self.runmetadata.samples:
            self.plusdict[sample.name] = dict()
            self.matchdict[sample.name] = dict()
            if sample.general.bestassemblyfile != 'NA':
                for gene in sample[self.analysistype].allelenames:
                    self.plusdict[sample.name][gene] = dict()
                    for allele, percentidentity in sample[self.analysistype].results.items():

                        if gene in allele:
                            # Split the allele number from the gene name using the appropriate delimiter
                            if '_' in allele:
                                splitter = '_'
                            elif '-' in allele:
                                splitter = '-'
                            else:
                                splitter = ''
                            self.matchdict[sample.name].update({gene: allele.split(splitter)[-1]})
                            # Create the plusdict dictionary as in the assembly-based (r)MLST method. Allows all the
                            # parsing and sequence typing code to be reused.
                            try:
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]][percentidentity] \
                                    = sample[self.analysistype].avgdepth[allele]
                            except KeyError:
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]] = dict()
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]][percentidentity] \
                                    = sample[self.analysistype].avgdepth[allele]
                    if gene not in self.matchdict[sample.name]:
                        self.matchdict[sample.name].update({gene: 'N'})
        self.profiler()
        self.sequencetyper()
        self.mlstreporter()