Run a variant comparison using the bcbio.variation toolkit, given an input configuration.
196829
def bcbio_variation_comparison(config_file, base_dir, data):
    
    tmp_dir = utils.safe_makedir(os.path.join(base_dir, "tmp"))
    resources = config_utils.get_resources("bcbio_variation", data["config"])
    jvm_opts = resources.get("jvm_opts", ["-Xms750m", "-Xmx2g"])
    cmd = ["bcbio-variation"] + jvm_opts + broad.get_default_jvm_opts(tmp_dir) + \
          ["variant-compare", config_file]
    do.run(cmd, "Comparing variant calls using bcbio.variation", data)
172604
def compare_config(self):
        
        # Set defaults
        base_file = "running-config"
        base_file_system = "system:"
        if self.config_replace:
            new_file = self.candidate_cfg
        else:
            new_file = self.merge_cfg
        new_file_system = self.dest_file_system

        base_file_full = self._gen_full_path(
            filename=base_file, file_system=base_file_system
        )
        new_file_full = self._gen_full_path(
            filename=new_file, file_system=new_file_system
        )

        if self.config_replace:
            cmd = "show archive config differences {} {}".format(
                base_file_full, new_file_full
            )
            diff = self.device.send_command_expect(cmd)
            diff = self._normalize_compare_config(diff)
        else:
            # merge
            cmd = "show archive config incremental-diffs {} ignorecase".format(
                new_file_full
            )
            diff = self.device.send_command_expect(cmd)
            if "error code 5" in diff or "returned error 5" in diff:
                diff = (
                    "You have encountered the obscure 'error 5' message. This generally "
                    "means you need to add an 'end' statement to the end of your merge changes."
                )
            elif "% Invalid" not in diff:
                diff = self._normalize_merge_diff_incr(diff)
            else:
                cmd = "more {}".format(new_file_full)
                diff = self.device.send_command_expect(cmd)
                diff = self._normalize_merge_diff(diff)

        return diff.strip()
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]]
1156483
def _cnvkit_prioritize(sample, genes, allele_file, metrics_file):
    
    mdf = pd.read_table(metrics_file)
    mdf.columns = [x.lower() for x in mdf.columns]
    if len(genes) > 0:
        mdf = mdf[mdf["gene"].str.contains("|".join(genes))]
    mdf = mdf[["chromosome", "start", "end", "gene", "log2", "ci_hi", "ci_lo"]]
    adf = pd.read_table(allele_file)
    if len(genes) > 0:
        adf = adf[adf["gene"].str.contains("|".join(genes))]
    if "cn1" in adf.columns and "cn2" in adf.columns:
        adf = adf[["chromosome", "start", "end", "cn", "cn1", "cn2"]]
    else:
        adf = adf[["chromosome", "start", "end", "cn"]]
    df = pd.merge(mdf, adf, on=["chromosome", "start", "end"])
    df = df[df["cn"] != 2]
    if len(df) > 0:
        def passes(row):
            spread = abs(row["ci_hi"] - row["ci_lo"])
            return spread < 0.25
        df["passes"] = df.apply(passes, axis=1)
    df.insert(0, "sample", [sample] * len(df))
    return df
996348
def run(self):
    """"""
    if self.exclude_from_doctest:
      for pattern in self.exclude_from_doctest:
        if fnmatch.fnmatch(self.filename, pattern):
          return

    checked_codes = pep257.conventions.pep257 | {'D998', 'D999'}
    for error in self._check_source():
      if isinstance(error, pep257.Error) and error.code in checked_codes:
        # NOTE(sigmavirus24): Fixes GitLab#3
        message = '%s %s' % (error.code, error.short_desc)
        yield (error.line, 0, message, type(self))
842874
def _check_variantcaller(item):
    
    allowed = set(list(genotype.get_variantcallers().keys()) + [None, False])
    vcs = item["algorithm"].get("variantcaller")
    if not isinstance(vcs, dict):
        vcs = {"variantcaller": vcs}
    for vc_set in vcs.values():
        if not isinstance(vc_set, (tuple, list)):
            vc_set = [vc_set]
        problem = [x for x in vc_set if x not in allowed]
        if len(problem) > 0:
            raise ValueError("Unexpected algorithm 'variantcaller' parameter: %s\n"
                             "Supported options: %s\n" % (problem, sorted(list(allowed))))
    # Ensure germline somatic calling only specified with tumor/normal samples
    if "germline" in vcs or "somatic" in vcs:
        paired = vcfutils.get_paired_phenotype(item)
        if not paired:
            raise ValueError("%s: somatic/germline calling in 'variantcaller' "
                             "but tumor/normal metadata phenotype not specified" % dd.get_sample_name(item))
192268
def _run_qsnp_paired(align_bams, items, ref_file, assoc_files,
                     region=None, out_file=None):
    
    config = items[0]["config"]
    if out_file is None:
        out_file = "%s-paired-variants.vcf" % os.path.splitext(align_bams[0])[0]
    if not utils.file_exists(out_file):
        out_file = out_file.replace(".gz", "")
        with file_transaction(config, out_file) as tx_out_file:
            with tx_tmpdir(config) as tmpdir:
                with utils.chdir(tmpdir):
                    paired = get_paired_bams(align_bams, items)
                    qsnp = config_utils.get_program("qsnp", config)
                    resources = config_utils.get_resources("qsnp", config)
                    mem = " ".join(resources.get("jvm_opts", ["-Xms750m -Xmx4g"]))
                    qsnp_log = os.path.join(tmpdir, "qsnp.log")
                    qsnp_init = os.path.join(tmpdir, "qsnp.ini")
                    if region:
                        paired = _create_bam_region(paired, region, tmpdir)
                    _create_input(paired, tx_out_file, ref_file, assoc_files['dbsnp'], qsnp_init)
                    cl = ("{qsnp} {mem} -i {qsnp_init} -log {qsnp_log}")
                    do.run(cl.format(**locals()), "Genotyping paired variants with Qsnp", {})
        out_file = _filter_vcf(out_file)
        out_file = bgzip_and_index(out_file, config)
    return out_file
632702
def _run_vardict_paired(align_bams, items, ref_file, assoc_files,
                          region=None, out_file=None):
    
    config = items[0]["config"]
    if out_file is None:
        out_file = "%s-paired-variants.vcf.gz" % os.path.splitext(align_bams[0])[0]
    if not utils.file_exists(out_file):
        with file_transaction(items[0], out_file) as tx_out_file:
            vrs = bedutils.population_variant_regions(items)
            target = shared.subset_variant_regions(vrs, region,
                                                   out_file, items=items, do_merge=True)
            paired = vcfutils.get_paired_bams(align_bams, items)
            if not _is_bed_file(target):
                vcfutils.write_empty_vcf(tx_out_file, config,
                                         samples=[x for x in [paired.tumor_name, paired.normal_name] if x])
            else:
                if not paired.normal_bam:
                    ann_file = _run_vardict_caller(align_bams, items, ref_file,
                                                   assoc_files, region, out_file)
                    return ann_file
                vardict = get_vardict_command(items[0])
                vcfstreamsort = config_utils.get_program("vcfstreamsort", config)
                compress_cmd = "| bgzip -c" if out_file.endswith("gz") else ""
                freq = float(utils.get_in(config, ("algorithm", "min_allele_fraction"), 10)) / 100.0
                # merge bed file regions as amplicon VarDict is only supported in single sample mode
                opts, var2vcf_opts = _vardict_options_from_config(items, config, out_file, target)
                fix_ambig_ref = vcfutils.fix_ambiguous_cl()
                fix_ambig_alt = vcfutils.fix_ambiguous_cl(5)
                remove_dup = vcfutils.remove_dup_cl()
                if any("vardict_somatic_filter" in tz.get_in(("config", "algorithm", "tools_off"), data, [])
                       for data in items):
                    somatic_filter = ""
                    freq_filter = ""
                else:
                    var2vcf_opts += " -M "  # this makes VarDict soft filter non-differential variants
                    somatic_filter = ("| sed 's/\\\\.*Somatic\\\\/Somatic/' "
                                      "| sed 's/REJECT,Description=\".*\">/REJECT,Description=\"Not Somatic via VarDict\">/' "
                                      """| %s -c 'from bcbio.variation import freebayes; freebayes.call_somatic("%s", "%s")' """
                                      % (sys.executable, paired.tumor_name, paired.normal_name))
                    freq_filter = ("| bcftools filter -m '+' -s 'REJECT' -e 'STATUS !~ \".*Somatic\"' 2> /dev/null "
                                   "| %s -x 'bcbio.variation.vardict.add_db_germline_flag(x)' "
                                   "| %s "
                                   "| %s -x 'bcbio.variation.vardict.depth_freq_filter(x, %s, \"%s\")'" %
                                   (os.path.join(os.path.dirname(sys.executable), "py"),
                                    _lowfreq_linear_filter(0, True),
                                    os.path.join(os.path.dirname(sys.executable), "py"),
                                    0, bam.aligner_from_header(paired.tumor_bam)))
                jvm_opts = _get_jvm_opts(items[0], tx_out_file)
                py_cl = os.path.join(utils.get_bcbio_bin(), "py")
                setup = ("%s && unset JAVA_HOME &&" % utils.get_R_exports())
                contig_cl = vcfutils.add_contig_to_header_cl(ref_file, tx_out_file)
                cmd = ("{setup}{jvm_opts}{vardict} -G {ref_file} "
                       "-N {paired.tumor_name} -b \"{paired.tumor_bam}|{paired.normal_bam}\" {opts} "
                       "| awk 'NF>=48' | testsomatic.R "
                       "| var2vcf_paired.pl -P 0.9 -m 4.25 {var2vcf_opts} "
                       "-N \"{paired.tumor_name}|{paired.normal_name}\" "
                       "| {contig_cl} {freq_filter} "
                       "| bcftools filter -i 'QUAL >= 0' "
                       "{somatic_filter} | {fix_ambig_ref} | {fix_ambig_alt} | {remove_dup} | {vcfstreamsort} "
                       "{compress_cmd} > {tx_out_file}")
                do.run(cmd.format(**locals()), "Genotyping with VarDict: Inference", {})
    return out_file
164584
def get_variant_by_name(self, name):
        
        results = []

        try:
            for info, dosage in self._bgen.get_variant(name):
                results.append(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=False,
                ))

        except ValueError:
            logging.variant_name_not_found(name)

        return results
165371
def main():
    ''''''
    args = docopt(__doc__, version=__version__)

    if 'bam_coverage' in args:
        bam_coverage(args['<reference>'],
                     args['<alignments>'],
                     int(args['<minmatch>']),
                     min_mapq=int(args['--mapq']),
                     min_len=float(args['--minlen']))
458703
def variant(samples):
    
    checkpoints = _variant_checkpoints(samples)
    if checkpoints["align"]:
        align_wf = _alignment(checkpoints)
        alignin = [["files"], ["analysis"],
                   ["config", "algorithm", "align_split_size"],
                   ["reference", "fasta", "base"],
                   ["rgnames", "pl"], ["rgnames", "sample"], ["rgnames", "pu"],
                   ["rgnames", "lane"], ["rgnames", "rg"], ["rgnames", "lb"],
                   ["reference", "aligner", "indexes"],
                   ["config", "algorithm", "aligner"],
                   ["config", "algorithm", "trim_reads"],
                   ["config", "algorithm", "adapters"],
                   ["config", "algorithm", "bam_clean"],
                   ["config", "algorithm", "variant_regions"],
                   ["config", "algorithm", "mark_duplicates"]]
        if checkpoints["hla"]:
            alignin.append(["config", "algorithm", "hlacaller"])
        if checkpoints["umi"]:
            alignin.append(["config", "algorithm", "umi_type"])
        align = [s("alignment_to_rec", "multi-combined", alignin,
                   [cwlout("alignment_rec", "record")],
                   "bcbio-vc",
                   disk={"files": 1.5}, cores=1, no_files=True),
                 w("alignment", "multi-parallel", align_wf,
                   [["align_split"], ["process_alignment_rec"],
                    ["work_bam"], ["config", "algorithm", "quality_format"]])]
    else:
        align = [s("organize_noalign", "multi-parallel",
                   ["files"],
                   [cwlout(["align_bam"], ["File", "null"], [".bai"]),
                    cwlout(["work_bam_plus", "disc"], ["File", "null"]),
                    cwlout(["work_bam_plus", "sr"], ["File", "null"]),
                    cwlout(["hla", "fastq"], ["File", "null"])],
                   "bcbio-vc", cores=1)]
    align_out = [["rgnames", "sample"], ["align_bam"]]
    pp_align, pp_align_out = _postprocess_alignment(checkpoints)
    if checkpoints["umi"]:
        align_out += [["umi_bam"]]
    vc, vc_out = _variant_vc(checkpoints)
    sv, sv_out = _variant_sv(checkpoints)
    hla, hla_out = _variant_hla(checkpoints)
    qc, qc_out = _qc_workflow(checkpoints)
    steps = align + pp_align + hla + vc + sv + qc
    final_outputs = align_out + pp_align_out + vc_out + hla_out + sv_out + qc_out
    return steps, final_outputs
766311
def run(config, path=None, stop_on_error=True, just_tests=False):
    """"""
    if path is None:
        path = os.getcwd()
    prep(config, path=path)

    models = __import__('models')
    tests = __import__('tests')
    suites = __import__('suites')

    print('\n')
    for x in ['models', 'tests', 'suites']:
        module = __import__(x)
        assert hasattr(module, x), "'%s' module requires attribute '%s'" %\
                                   (x, x)

    if just_tests:
        for test in tests.tests:
            _run(test, models, stop_on_error)

    else:
        for suite in suites.suites:
            _run(suite, models, stop_on_error)
780305
def run(config_file):
	""""""
	# Load configuration.
	config = neat.Config(neat.DefaultGenome, neat.DefaultReproduction,
						 neat.DefaultSpeciesSet, neat.DefaultStagnation,
						 config_file)

	# Create the population, which is the top-level object for a NEAT run.
	p = neat.Population(config)

	# Add a stdout reporter to show progress in the terminal.
	p.add_reporter(neat.StdOutReporter(True))
	stats = neat.StatisticsReporter()
	p.add_reporter(stats)

	# Run for up to 300 generations.
	pe = neat.ThreadedEvaluator(4, eval_genome)
	winner = p.run(pe.evaluate, 300)
	pe.stop()

	# Display the winning genome.
	print('\nBest genome:\n{!s}'.format(winner))

	# Show output of the most fit genome against training data.
	print('\nOutput:')
	winner_net = neat.nn.FeedForwardNetwork.create(winner, config)
	for xi, xo in zip(xor_inputs, xor_outputs):
		output = winner_net.activate(xi)
		print(
			"input {!r}, expected output {!r}, got {!r}".format(xi, xo, output)
			)

	if visualize is not None:
		node_names = {-1: 'A', -2: 'B', 0: 'A XOR B'}
		visualize.draw_net(config, winner, True, node_names=node_names)
		visualize.plot_stats(stats, ylog=False, view=True)
		visualize.plot_species(stats, view=True)
661146
def _run_concat_variant_files_bcftools(in_list, out_file, config, naive=False):
    
    if not utils.file_exists(out_file):
        with file_transaction(config, out_file) as tx_out_file:
            bcftools = config_utils.get_program("bcftools", config)
            output_type = "z" if out_file.endswith(".gz") else "v"
            if naive:
                args = "--naive"
            else:
                args = "--allow-overlaps"
            cmd = "{bcftools} concat {args} -O {output_type} --file-list {in_list} -o {tx_out_file}"
            do.run(cmd.format(**locals()), "bcftools concat variants")
    if out_file.endswith(".gz"):
        bgzip_and_index(out_file, config)
    return out_file
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
1150832
def run(configobj=None,editpars=False):
    

    if configobj is None:
        configobj =teal.teal(__taskname__,loadOnly=(not editpars))

    update(configobj['input'],configobj['refdir'],
        local=configobj['local'],interactive=configobj['interactive'],
        wcsupdate=configobj['wcsupdate'])
1344741
def run_haplotyper(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):
        variant_regions = bedutils.population_variant_regions(items, merged=True)
        interval = _get_interval(variant_regions, region, out_file, items)
        with file_transaction(items[0], out_file) as tx_out_file:
            dbsnp = "--dbsnp %s" % (assoc_files.get("dbsnp")) if "dbsnp" in assoc_files else ""
            bams = " ".join(["-i %s" % x for x in align_bams])
            license = license_export(items[0])
            cores = dd.get_num_cores(items[0])
            out_mode = "--emit_mode gvcf" if joint.want_gvcf(items) else ""
            cmd = ("{license}sentieon driver -t {cores} -r {ref_file} "
                   "{bams} {interval} --algo Haplotyper {out_mode} {dbsnp} {tx_out_file}")
            do.run(cmd.format(**locals()), "Sentieon Haplotyper")
    return out_file