<commit_before>"""Calculate potential effects of variations using external programs.

Supported:
  snpEff: http://sourceforge.net/projects/snpeff/
  VEP: http://www.ensembl.org/info/docs/tools/vep/index.html
"""
import os
import glob
import shutil
import subprocess

import toolz as tz
import yaml

from bcbio import utils
from bcbio.distributed.transaction import file_transaction
from bcbio.pipeline import config_utils, tools
from bcbio.provenance import do
from bcbio.variation import vcfutils

# ## Ensembl VEP

def vep_version(config):
    try:
        vep = config_utils.get_program("variant_effect_predictor.pl", config)
        help_str = subprocess.check_output([vep, "--help"])
        for line in help_str.split("\n"):
            if line.startswith("version"):
                return line.split()[-1].strip()
        return None
    except config_utils.CmdNotFound:
        return None
        return False

def _special_dbkey_maps(dbkey, ref_file):
    """Avoid duplicate VEP information for databases with chromosome differences like hg19/GRCh37.
    """
    remaps = {"hg19": "GRCh37"}
    if dbkey in remaps:
        base_dir = os.path.normpath(os.path.join(os.path.dirname(ref_file), os.pardir))
        vep_dir = os.path.normpath(os.path.join(base_dir, "vep"))
        other_dir = os.path.relpath(os.path.normpath(os.path.join(base_dir, os.pardir, remaps[dbkey], "vep")),
                                    base_dir)
        if not os.path.lexists(vep_dir):
            os.symlink(other_dir, vep_dir)
        return vep_dir
    else:
        return None

def prep_vep_cache(dbkey, ref_file, config=None):
    """Ensure correct installation of VEP cache file.
    """
    if config is None: config = {}
    resource_file = os.path.join(os.path.dirname(ref_file), "%s-resources.yaml" % dbkey)
    vepv = vep_version(config)
    if os.path.exists(resource_file) and vepv:
        with open(resource_file) as in_handle:
            resources = yaml.load(in_handle)
        ensembl_name = tz.get_in(["aliases", "ensembl"], resources)
        symlink_dir = _special_dbkey_maps(dbkey, ref_file)
        if symlink_dir:
            return symlink_dir, ensembl_name
        elif ensembl_name:
            vep_dir = utils.safe_makedir(os.path.normpath(os.path.join(
                os.path.dirname(os.path.dirname(ref_file)), "vep")))
            out_dir = os.path.join(vep_dir, ensembl_name, vepv)
            if not os.path.exists(out_dir):
                cmd = ["vep_install.pl", "-a", "c", "-s", ensembl_name,
                       "-c", vep_dir]
                do.run(cmd, "Prepare VEP directory for %s" % ensembl_name)
                cmd = ["vep_convert_cache.pl", "-species", ensembl_name, "-version", vepv,
                       "-d", vep_dir]
                do.run(cmd, "Convert VEP cache to tabix %s" % ensembl_name)
            tmp_dir = os.path.join(vep_dir, "tmp")
            if os.path.exists(tmp_dir):
                shutil.rmtree(tmp_dir)
            return vep_dir, ensembl_name
    return None, None

def run_vep(data):
    """Annotate input VCF file with Ensembl variant effect predictor.
    """
    out_file = utils.append_stem(data["vrn_file"], "-vepeffects")
    assert data["vrn_file"].endswith(".gz") and out_file.endswith(".gz")
    if not utils.file_exists(out_file):
        with file_transaction(out_file) as tx_out_file:
            vep_dir, ensembl_name = prep_vep_cache(data["genome_build"],
                                                   tz.get_in(["reference", "fasta", "base"], data))
            if vep_dir:
                vep = config_utils.get_program("variant_effect_predictor.pl", data["config"])
                cmd = [vep, "--vcf", "-o", "stdout", "--fork", "4",
                       "--species", ensembl_name,
                       "--no_stats",
                       "--cache", "--offline", "--dir", vep_dir,
                       "--sift", "b", "--polyphen", "b", "--symbol", "--numbers", "--biotype", "--total_length",
                       "--fields",
                       "Consequence,Codons,Amino_acids,Gene,SYMBOL,Feature,EXON,PolyPhen," +
                       "SIFT,Protein_position,BIOTYPE"] + _get_dbnsfp(data)
                cmd = "gunzip -c %s | %s | bgzip -c > %s" % (data["vrn_file"], " ".join(cmd), tx_out_file)
                do.run(cmd, "Ensembl variant effect predictor", data)
    if utils.file_exists(out_file):
        vcfutils.bgzip_and_index(out_file, data["config"])
        return out_file

def _get_dbnsfp(data):
    """Retrieve dbNSFP file options for VEP if downloaded and available.

    Uses high level combined annotations from this GEMINI discussion as a
    starting point:
    https://groups.google.com/d/msg/gemini-variation/WeZ6C2YvfUA/mII9uum_pGoJ
    """
    dbnsfp_file = tz.get_in(("genome_resources", "variation", "dbnsfp"), data)
    if dbnsfp_file and os.path.exists(dbnsfp_file):
        annotations = ["RadialSVM_score", "RadialSVM_pred", "LR_score", "LR_pred",
                       "CADD_raw", "CADD_phred", "Reliability_index"]
        return ["--plugin", "dbNSFP,%s,%s" % (dbnsfp_file, ",".join(annotations))]
    else:
        return []

# ## snpEff variant effects

def snpeff_effects(data):
    """Annotate input VCF file with effects calculated by snpEff.
    """
    vcf_in = data["vrn_file"]
    if vcfutils.vcf_has_variants(vcf_in):
        vcf_file = _run_snpeff(vcf_in, "vcf", data)
        return vcf_file

def _snpeff_args_from_config(data):
    """Retrieve snpEff arguments supplied through input configuration.
    """
    config = data["config"]
    args = []
    # General supplied arguments
    resources = config_utils.get_resources("snpeff", config)<commit_msg>VEP: Add canonical and CCDS annotations for prioritizing transcripts. Ensure dbNSFP annotations added to CSQ tag features<commit_after>"""Calculate potential effects of variations using external programs.

Supported:
  snpEff: http://sourceforge.net/projects/snpeff/
  VEP: http://www.ensembl.org/info/docs/tools/vep/index.html
"""
import os
import glob
import shutil
import subprocess

import toolz as tz
import yaml

from bcbio import utils
from bcbio.distributed.transaction import file_transaction
from bcbio.pipeline import config_utils, tools
from bcbio.provenance import do
from bcbio.variation import vcfutils

# ## Ensembl VEP

def vep_version(config):
    try:
        vep = config_utils.get_program("variant_effect_predictor.pl", config)
        help_str = subprocess.check_output([vep, "--help"])
        for line in help_str.split("\n"):
            if line.startswith("version"):
                return line.split()[-1].strip()
        return None
    except config_utils.CmdNotFound:
        return None
        return False

def _special_dbkey_maps(dbkey, ref_file):
    """Avoid duplicate VEP information for databases with chromosome differences like hg19/GRCh37.
    """
    remaps = {"hg19": "GRCh37"}
    if dbkey in remaps:
        base_dir = os.path.normpath(os.path.join(os.path.dirname(ref_file), os.pardir))
        vep_dir = os.path.normpath(os.path.join(base_dir, "vep"))
        other_dir = os.path.relpath(os.path.normpath(os.path.join(base_dir, os.pardir, remaps[dbkey], "vep")),
                                    base_dir)
        if not os.path.lexists(vep_dir):
            os.symlink(other_dir, vep_dir)
        return vep_dir
    else:
        return None

def prep_vep_cache(dbkey, ref_file, config=None):
    """Ensure correct installation of VEP cache file.
    """
    if config is None: config = {}
    resource_file = os.path.join(os.path.dirname(ref_file), "%s-resources.yaml" % dbkey)
    vepv = vep_version(config)
    if os.path.exists(resource_file) and vepv:
        with open(resource_file) as in_handle:
            resources = yaml.load(in_handle)
        ensembl_name = tz.get_in(["aliases", "ensembl"], resources)
        symlink_dir = _special_dbkey_maps(dbkey, ref_file)
        if symlink_dir:
            return symlink_dir, ensembl_name
        elif ensembl_name:
            vep_dir = utils.safe_makedir(os.path.normpath(os.path.join(
                os.path.dirname(os.path.dirname(ref_file)), "vep")))
            out_dir = os.path.join(vep_dir, ensembl_name, vepv)
            if not os.path.exists(out_dir):
                cmd = ["vep_install.pl", "-a", "c", "-s", ensembl_name,
                       "-c", vep_dir]
                do.run(cmd, "Prepare VEP directory for %s" % ensembl_name)
                cmd = ["vep_convert_cache.pl", "-species", ensembl_name, "-version", vepv,
                       "-d", vep_dir]
                do.run(cmd, "Convert VEP cache to tabix %s" % ensembl_name)
            tmp_dir = os.path.join(vep_dir, "tmp")
            if os.path.exists(tmp_dir):
                shutil.rmtree(tmp_dir)
            return vep_dir, ensembl_name
    return None, None

def run_vep(data):
    """Annotate input VCF file with Ensembl variant effect predictor.
    """
    out_file = utils.append_stem(data["vrn_file"], "-vepeffects")
    assert data["vrn_file"].endswith(".gz") and out_file.endswith(".gz")
    if not utils.file_exists(out_file):
        with file_transaction(out_file) as tx_out_file:
            vep_dir, ensembl_name = prep_vep_cache(data["genome_build"],
                                                   tz.get_in(["reference", "fasta", "base"], data))
            if vep_dir:
                vep = config_utils.get_program("variant_effect_predictor.pl", data["config"])
                dbnsfp_args, dbnsfp_fields = _get_dbnsfp(data)
                std_fields = ["Consequence", "Codons", "Amino_acids", "Gene", "SYMBOL", "Feature",
                              "EXON", "PolyPhen", "SIFT", "Protein_position", "BIOTYPE", "CANONICAL", "CCDS"]
                cmd = [vep, "--vcf", "-o", "stdout", "--fork", "4",
                       "--species", ensembl_name,
                       "--no_stats",
                       "--cache", "--offline", "--dir", vep_dir,
                       "--sift", "b", "--polyphen", "b", "--symbol", "--numbers", "--biotype", "--total_length",
                       "--canonical", "--ccds",
                       "--fields", ",".join(std_fields + dbnsfp_fields)] + dbnsfp_args
                cmd = "gunzip -c %s | %s | bgzip -c > %s" % (data["vrn_file"], " ".join(cmd), tx_out_file)
                do.run(cmd, "Ensembl variant effect predictor", data)
    if utils.file_exists(out_file):
        vcfutils.bgzip_and_index(out_file, data["config"])
        return out_file

def _get_dbnsfp(data):
    """Retrieve dbNSFP file options for VEP if downloaded and available.

    Uses high level combined annotations from this GEMINI discussion as a
    starting point:
    https://groups.google.com/d/msg/gemini-variation/WeZ6C2YvfUA/mII9uum_pGoJ
    """
    dbnsfp_file = tz.get_in(("genome_resources", "variation", "dbnsfp"), data)
    if dbnsfp_file and os.path.exists(dbnsfp_file):
        annotations = ["RadialSVM_score", "RadialSVM_pred", "LR_score", "LR_pred",
                       "CADD_raw", "CADD_phred", "Reliability_index"]
        return ["--plugin", "dbNSFP,%s,%s" % (dbnsfp_file, ",".join(annotations))], annotations
    else:
        return [], []

# ## snpEff variant effects

def snpeff_effects(data):
    """Annotate input VCF file with effects calculated by snpEff.
    """
    vcf_in = data["vrn_file"]
    if vcfutils.vcf_has_variants(vcf_in):
        vcf_file = _run_snpeff(vcf_in, "vcf", data)
        return vcf_file

def _snpeff_args_from_config(data):
    """Retrieve snpEff arguments supplied through input configuration.
    """
    config = data["config"]
    args = []
    # General supplied arguments
    resources = config_utils.get_resources("snpeff", config)
{'fraction_non_alphanumeric': 0.1283427014137285, 'fraction_numerical': 0.0022142735479475386, 'mean_word_length': 2.5338549503460728, 'pattern_counts': {'":': 2, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4861463', 'n_tokens_mistral': 3914, 'n_tokens_neox': 3802, 'n_words': 827}