Generalized vcfanno/vcf2db workflow for loading variants into a GEMINI database.
965979
def create_gemini_db(gemini_vcf, data, gemini_db=None, ped_file=None):
    
    if not gemini_db:
        gemini_db = "%s.db" % utils.splitext_plus(gemini_vcf)[0]
    if not vcfutils.vcf_has_variants(gemini_vcf):
        return None
    if not utils.file_exists(gemini_db):
        with file_transaction(data, gemini_db) as tx_gemini_db:
            vcf2db = config_utils.get_program("vcf2db.py", data)
            if "vcf2db_expand" in dd.get_tools_on(data):
                vcf2db_args = ["--expand", "gt_types", "--expand", "gt_ref_depths", "--expand", "gt_alt_depths"]
            else:
                vcf2db_args = []
            cmd = [vcf2db, gemini_vcf, ped_file, tx_gemini_db] + vcf2db_args
            do.run(cmd, "GEMINI: create database with vcf2db")
    return gemini_db
579009
def variation(fastas, sams, s2b, annots, qual, min_cov, min_per, no_sub, bp_table, print_con):
    
    genomes = {}

    # save base-pair data structures for each genome-sample pair
    for sam in sams:
        genomes, s2b = save_refs(sam, fastas, genomes, s2b)

    # get base-pair frequencies for each genome-sample pair
    for sam in sams:
        genomes = get_bp_stats(sam, genomes, s2b, qual)

    # calculate base-pair frequencies
    genomes = calc_frequencies(genomes, bp_table, min_cov, min_per)

    # print consensus genome
    if print_con is True:
        print('# saving consensus sequences', file = sys.stderr)
        genomes = print_consensus(genomes)

    if no_sub is True:
        return genomes

    # calculate substitution rates
    for genome, contigs in list(genomes.items()):
        for contig, samples in list(contigs.items()):
            for sample in samples:
#                print(genomes.keys())
#                print(genomes[genome].keys())
#                print(genomes[genome][contig].keys())
#                print(annots.keys())
                genomes[genome][contig][sample]['sub_rates'] = \
                    sub_rates(genomes[genome][contig][sample]['bp_stats'], annots[contig])

    return genomes
809612
def _get_build_type(fnames, samples, caller):
    
    build_type = set()
    if any(vcfutils.vcf_has_variants(f) for f in fnames) and caller not in NO_DB_CALLERS:
        for data in samples:
            if any([x in dd.get_tools_on(data)
                    for x in ["gemini", "gemini_orig", "gemini_allvariants", "vcf2db_expand"]]):
                if vcfanno.annotate_gemini(data):
                    build_type.add("gemini_orig" if "gemini_orig" in dd.get_tools_on(data) else "gemini")
                else:
                    logger.info("Not running gemini, input data not found: %s" % dd.get_sample_name(data))
            else:
                logger.info("Not running gemini, not configured in tools_on: %s" % dd.get_sample_name(data))
    else:
        logger.info("Not running gemini, no samples with variants found: %s" %
                    (", ".join([dd.get_sample_name(d) for d in samples])))
    return build_type
774878
def finalize_vcf(in_file, variantcaller, items):
    
    out_file = "%s-annotated%s" % utils.splitext_plus(in_file)
    if not utils.file_uptodate(out_file, in_file):
        header_cl = _add_vcf_header_sample_cl(in_file, items, out_file)
        contig_cl = _add_contig_cl(in_file, items, out_file)
        cls = [x for x in (contig_cl, header_cl) if x]
        if cls:
            post_cl = " | ".join(cls) + " | "
        else:
            post_cl = None
        dbsnp_file = tz.get_in(("genome_resources", "variation", "dbsnp"), items[0])
        if dbsnp_file:
            out_file = _add_dbsnp(in_file, dbsnp_file, items[0], out_file, post_cl)
    if utils.file_exists(out_file):
        return vcfutils.bgzip_and_index(out_file, items[0]["config"])
    else:
        return in_file
368479
def _add_transcripts(self, variant_obj, gemini_variant):
        
        query = "SELECT * from variant_impacts WHERE variant_id = {0}".format(
            gemini_variant['variant_id']
        )

        gq = GeminiQuery(self.db)
        gq.run(query)
        
        for gemini_transcript in gq:
            transcript = Transcript(
                hgnc_symbol=gemini_transcript['gene'],
                transcript_id=gemini_transcript['transcript'],
                consequence=gemini_transcript['impact_so'],
                biotype=gemini_transcript['biotype'],
                polyphen=gemini_transcript['polyphen_pred'],
                sift=gemini_transcript['sift_pred'],
                HGVSc=gemini_transcript['codon_change'],
                HGVSp=', '.join([gemini_transcript['aa_change'] or '', gemini_transcript['aa_length'] or ''])
                )
            variant_obj.add_transcript(transcript)
137475
def compilevcf(args):
    
    p = OptionParser(compilevcf.__doc__)
    p.add_option("--db", default="hg38", help="Use these lobSTR db")
    p.add_option("--nofilter", default=False, action="store_true",
                 help="Do not filter the variants")
    p.set_home("lobstr")
    p.set_cpus()
    p.set_aws_opts(store="hli-mv-data-science/htang/str-data")
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    samples, = args
    workdir = opts.workdir
    store = opts.output_path
    cleanup = not opts.nocleanup
    filtered = not opts.nofilter
    dbs = opts.db.split(",")
    cwd = os.getcwd()
    mkdir(workdir)
    os.chdir(workdir)
    samples = op.join(cwd, samples)

    stridsfile = "STR.ids"
    if samples.endswith((".vcf", ".vcf.gz")):
        vcffiles = [samples]
    else:
        vcffiles = [x.strip() for x in must_open(samples)]
    if not op.exists(stridsfile):
        ids = []
        for db in dbs:
            ids.extend(STRFile(opts.lobstr_home, db=db).ids)
        uids = uniqify(ids)
        logging.debug("Combined: {} Unique: {}".format(len(ids), len(uids)))

        fw = open(stridsfile, "w")
        print("\n".join(uids), file=fw)
        fw.close()

    run_args = [(x, filtered, cleanup, store) for x in vcffiles]
    cpus = min(opts.cpus, len(run_args))
    p = Pool(processes=cpus)
    for res in p.map_async(run_compile, run_args).get():
        continue
892685
def reads_to_dataframe(variants_and_allele_reads):
    
    df_builder = DataFrameBuilder(
        AlleleRead,
        extra_column_fns={
            "gene": lambda variant, _: ";".join(variant.gene_names),
        })
    for variant, allele_reads in variants_and_allele_reads:
        df_builder.add_many(variant, allele_reads)
    return df_builder.to_dataframe()
92098
def _get_vcf_breakends(hydra_file, genome_2bit, options=None):
    
    if options is None: options = {}
    for features in group_hydra_breakends(hydra_parser(hydra_file, options)):
        if len(features) == 1 and is_deletion(features[0], options):
            yield build_vcf_deletion(features[0], genome_2bit)
        elif len(features) == 1 and is_tandem_dup(features[0], options):
            yield build_tandem_deletion(features[0], genome_2bit)
        elif len(features) == 2 and is_inversion(*features):
            yield build_vcf_inversion(features[0], features[1], genome_2bit)
        elif len(features) == 2 and is_translocation(*features):
            info = get_translocation_info(features[0], features[1])
            for feature in features:
                for brend in build_vcf_parts(feature, genome_2bit, info):
                    yield brend
        else:
            for feature in features:
                for brend in build_vcf_parts(feature, genome_2bit):
                    yield brend
459244
def variant_sequences_generator_to_dataframe(variant_sequences_generator):
    
    # TODO: Change VariantSequence.alt to VariantSequence.alt_nucleotides
    # or something else that doesn't clash with a variant's `alt` field
    return dataframe_from_generator(
        VariantSequence,
        variant_sequences_generator,
        rename_dict={"alt": "allele"},
        extra_column_fns={
            "gene": lambda variant, _: ";".join(variant.gene_names),
        })
95450
def postprocess_variants(items):
    
    vrn_key = "vrn_file"
    if not isinstance(items, dict):
        items = [utils.to_single_data(x) for x in items]
        if "vrn_file_joint" in items[0]:
            vrn_key = "vrn_file_joint"
    data, items = _get_batch_representative(items, vrn_key)
    items = cwlutils.unpack_tarballs(items, data)
    data = cwlutils.unpack_tarballs(data, data)
    cur_name = "%s, %s" % (dd.get_sample_name(data), get_variantcaller(data, require_bam=False))
    logger.info("Finalizing variant calls: %s" % cur_name)
    orig_vrn_file = data.get(vrn_key)
    data = _symlink_to_workdir(data, [vrn_key])
    data = _symlink_to_workdir(data, ["config", "algorithm", "variant_regions"])
    if data.get(vrn_key):
        logger.info("Calculating variation effects for %s" % cur_name)
        ann_vrn_file, vrn_stats = effects.add_to_vcf(data[vrn_key], data)
        if ann_vrn_file:
            data[vrn_key] = ann_vrn_file
        if vrn_stats:
            data["vrn_stats"] = vrn_stats
        orig_items = _get_orig_items(items)
        logger.info("Annotate VCF file: %s" % cur_name)
        data[vrn_key] = annotation.finalize_vcf(data[vrn_key], get_variantcaller(data, require_bam=False),
                                                orig_items)
        if cwlutils.is_cwl_run(data):
            logger.info("Annotate with population level variation data")
            ann_file = population.run_vcfanno(data[vrn_key], data)
            if ann_file:
                data[vrn_key] = ann_file
        logger.info("Filtering for %s" % cur_name)
        data[vrn_key] = variant_filtration(data[vrn_key], dd.get_ref_file(data),
                                           tz.get_in(("genome_resources", "variation"), data, {}),
                                           data, orig_items)
        logger.info("Prioritization for %s" % cur_name)
        prio_vrn_file = prioritize.handle_vcf_calls(data[vrn_key], data, orig_items)
        if prio_vrn_file != data[vrn_key]:
            data[vrn_key] = prio_vrn_file
            logger.info("Germline extraction for %s" % cur_name)
            data = germline.extract(data, orig_items)

        if dd.get_align_bam(data):
            data = damage.run_filter(data[vrn_key], dd.get_align_bam(data), dd.get_ref_file(data),
                                     data, orig_items)
    if orig_vrn_file and os.path.samefile(data[vrn_key], orig_vrn_file):
        data[vrn_key] = orig_vrn_file
    return [[data]]
89495
def get_variants_in_region(self, chrom, start, end):
        """"""
        bim = self.bim.loc[
            (self.bim["chrom"] == CHROM_STR_TO_INT[chrom]) &
            (start <= self.bim["pos"]) &
            (self.bim["pos"] <= end)
        ]
        for i, g in enumerate(self.bed.iter_geno_marker(bim.index)):
            info = bim.iloc[i, :]
            name, geno = g
            yield Genotypes(
                Variant(info.name, CHROM_INT_TO_STR[info.chrom],
                        info.pos, [info.a1, info.a2]),
                self._normalize_missing(geno),
                reference=info.a2,
                coded=info.a1,
                multiallelic=info.multiallelic
            )
857431
def build_vcf_parts(feature, genome_2bit, info=None):
    
    base1 = genome_2bit[feature.chrom1].get(
        feature.start1, feature.start1 + 1).upper()
    id1 = "hydra{0}a".format(feature.name)
    base2 = genome_2bit[feature.chrom2].get(
        feature.start2, feature.start2 + 1).upper()
    id2 = "hydra{0}b".format(feature.name)
    orientation = _breakend_orientation(feature.strand1, feature.strand2)
    return (VcfLine(feature.chrom1, feature.start1, id1, base1,
                    _vcf_alt(base1, feature.chrom2, feature.start2,
                             orientation.is_rc1, orientation.is_first1),
                    _vcf_info(feature.start1, feature.end1, id2, info)),
            VcfLine(feature.chrom2, feature.start2, id2, base2,
                    _vcf_alt(base2, feature.chrom1, feature.start1,
                             orientation.is_rc2, orientation.is_first2),
                    _vcf_info(feature.start2, feature.end2, id1, info)))
270084
def variants(ctx, variant_id, chromosome, end_chromosome, start, end, variant_type,
             sv_type):
    """"""
    if sv_type:
        variant_type = 'sv'

    adapter = ctx.obj['adapter']

    if (start or end):
        if not (chromosome and start and end):
            LOG.warning("Regions must be specified with chromosome, start and end")
            return

    if variant_id:
        variant = adapter.get_variant({'_id':variant_id})
        if variant:
            click.echo(variant)
        else:
            LOG.info("Variant {0} does not exist in database".format(variant_id))
        return

    if variant_type == 'snv':
        result = adapter.get_variants(
            chromosome=chromosome,
            start=start,
            end=end
        )
    else:
        LOG.info("Search for svs")
        result = adapter.get_sv_variants(
            chromosome=chromosome,
            end_chromosome=end_chromosome,
            sv_type=sv_type,
            pos=start,
            end=end
        )

    i = 0
    for variant in result:
        i += 1
        pp(variant)

    LOG.info("Number of variants found in database: %s", i)
1219693
def get_multisample_vcf(fnames, name, caller, data):
    
    unique_fnames = []
    for f in fnames:
        if f not in unique_fnames:
            unique_fnames.append(f)
    out_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "gemini"))
    if len(unique_fnames) > 1:
        gemini_vcf = os.path.join(out_dir, "%s-%s.vcf.gz" % (name, caller))
        vrn_file_batch = None
        for variant in data.get("variants", []):
            if variant["variantcaller"] == caller and variant.get("vrn_file_batch"):
                vrn_file_batch = variant["vrn_file_batch"]
        if vrn_file_batch:
            utils.symlink_plus(vrn_file_batch, gemini_vcf)
            return gemini_vcf
        else:
            return vcfutils.merge_variant_files(unique_fnames, gemini_vcf, dd.get_ref_file(data),
                                                data["config"])
    else:
        gemini_vcf = os.path.join(out_dir, "%s-%s%s" % (name, caller, utils.splitext_plus(unique_fnames[0])[1]))
        utils.symlink_plus(unique_fnames[0], gemini_vcf)
        return gemini_vcf
767030
def _back_compatible_gemini(conf_files, data):
    
    if vcfanno.is_human(data, builds=["37"]):
        for f in conf_files:
            if f and os.path.basename(f) == "gemini.conf" and os.path.exists(f):
                with open(f) as in_handle:
                    for line in in_handle:
                        if line.startswith("file"):
                            fname = line.strip().split("=")[-1].replace('"', '').strip()
                            if fname.find(".tidy.") > 0:
                                return install.get_gemini_dir(data)
    return None
874641
def _run_snpeff(snp_in, out_format, data):
    
    snpeff_db, datadir = get_db(data)
    if not snpeff_db:
        return None, None

    assert os.path.exists(os.path.join(datadir, snpeff_db)), \
        "Did not find %s snpEff genome data in %s" % (snpeff_db, datadir)
    ext = utils.splitext_plus(snp_in)[1] if out_format == "vcf" else ".tsv"
    out_file = "%s-effects%s" % (utils.splitext_plus(snp_in)[0], ext)
    stats_file = "%s-stats.html" % utils.splitext_plus(out_file)[0]
    csv_file = "%s-stats.csv" % utils.splitext_plus(out_file)[0]
    if not utils.file_exists(out_file):
        config_args = " ".join(_snpeff_args_from_config(data))
        if ext.endswith(".gz"):
            bgzip_cmd = "| %s -c" % tools.get_bgzip_cmd(data["config"])
        else:
            bgzip_cmd = ""
        with file_transaction(data, out_file) as tx_out_file:
            snpeff_cmd = _get_snpeff_cmd("eff", datadir, data, tx_out_file)
            cmd = ("{snpeff_cmd} {config_args} -noLog -i vcf -o {out_format} "
                   "-csvStats {csv_file} -s {stats_file} {snpeff_db} {snp_in} {bgzip_cmd} > {tx_out_file}")
            do.run(cmd.format(**locals()), "snpEff effects", data)
    if ext.endswith(".gz"):
        out_file = vcfutils.bgzip_and_index(out_file, data["config"])
    return out_file, [stats_file, csv_file]
732476
def run_vcfanno(vcf_file, data, decomposed=False):
    
    conf_files = dd.get_vcfanno(data)
    if conf_files:
        with_basepaths = collections.defaultdict(list)
        gemini_basepath = _back_compatible_gemini(conf_files, data)
        for f in conf_files:
            name = os.path.splitext(os.path.basename(f))[0]
            if f.endswith(".lua"):
                conf_file = None
                lua_file = f
            else:
                conf_file = f
                lua_file = "%s.lua" % utils.splitext_plus(conf_file)[0]
            if lua_file and not os.path.exists(lua_file):
                lua_file = None
            data_basepath = gemini_basepath if name == "gemini" else None
            if conf_file and os.path.exists(conf_file):
                with_basepaths[(data_basepath, name)].append(conf_file)
            if lua_file and os.path.exists(lua_file):
                with_basepaths[(data_basepath, name)].append(lua_file)
        conf_files = with_basepaths.items()
    out_file = None
    if conf_files:
        VcfannoIn = collections.namedtuple("VcfannoIn", ["conf", "lua"])
        bp_files = collections.defaultdict(list)
        for (data_basepath, name), anno_files in conf_files:
            anno_files = list(set(anno_files))
            if len(anno_files) == 1:
                cur = VcfannoIn(anno_files[0], None)
            elif len(anno_files) == 2:
                lua_files = [x for x in anno_files if x.endswith(".lua")]
                assert len(lua_files) == 1, anno_files
                lua_file = lua_files[0]
                anno_files.remove(lua_file)
                cur = VcfannoIn(anno_files[0], lua_file)
            else:
                raise ValueError("Unexpected annotation group %s" % anno_files)
            bp_files[data_basepath].append(cur)
        for data_basepath, anno_files in bp_files.items():
            ann_file = vcfanno.run(vcf_file, [x.conf for x in anno_files],
                                   [x.lua for x in anno_files], data,
                                   basepath=data_basepath,
                                   decomposed=decomposed)
            if ann_file:
                out_file = ann_file
                vcf_file = ann_file
    return out_file