notebooks/variant_effect_prediction.ipynb

{
 "cells": [
  {
   "source": [
    "# Variant effect prediction\n",
    "Variant effect prediction offers a simple way to predict effects of SNVs using any model that uses DNA sequence as an input. Many different scoring methods can be chosen, but the principle relies on in-silico mutagenesis. The default input is a VCF and the default output again is a TSV annotated with predictions of variant effects. \n",
    "\n",
    "This iPython notebook goes through the basic programmatic steps that are needed to preform variant effect prediction."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variant centered effect prediction\n",
    "Models that uses `kipoiseq.dataloaders.SeqIntervalDl` as a default dataloader can make use of variant-centered effect prediction. This procedure starts out from the query VCF and generates genomic regions of the length of the model input, centered on the individual variant in the VCF which are then mutated according to the alleles in the VCF. The model batch prediction function is then triggered for all mutated sequence sets and finally the scoring method is applied.\n",
    "\n",
    "The selected scoring methods compare model predicitons for sequences carrying the reference or alternative allele. Those scoring methods can be `Diff` for simple subtraction of prediction, `Logit` for substraction of logit-transformed model predictions, or `DeepSEA_effect` which is a combination of `Diff` and `Logit`, which was published in the Troyanskaya et al. (2015) publication.\n",
    "\n",
    "This ipython notebook assumes that it is executed in an environment in which kipoi-veff2 is installed. For more information check https://github.com/kipoi/kipoi-veff2#install-the-conda-environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Using downloaded and verified file: /Users/b260/.kipoi/models/DeepSEA/variantEffects/downloaded/model_files/weights/35956ab9c28960b5a3693f470fe980c1\n"
     ]
    }
   ],
   "source": [
    "from kipoi_veff2 import variant_centered\n",
    "\n",
    "vcf_file = \"example_data/clinvar_donor_acceptor_chr22.vcf\"\n",
    "fasta_file = \"example_data/hg19_chr22.fa\"\n",
    "output_file = \"output.tsv\"\n",
    "model_name = \"DeepSEA/variantEffects\"\n",
    "\n",
    "model_group = model_name.split(\"/\")[0]\n",
    "model_group_config_dict = (\n",
    "    variant_centered.VARIANT_CENTERED_MODEL_GROUP_CONFIGS.get(\n",
    "        model_group, {}\n",
    "    )\n",
    ")\n",
    "\n",
    "model_config = variant_centered.get_model_config(\n",
    "    model_name, **model_group_config_dict\n",
    ")\n",
    "\n",
    "variant_centered.score_variants(\n",
    "    model_config=model_config,\n",
    "    vcf_file=vcf_file,\n",
    "    fasta_file=fasta_file,\n",
    "    output_file=output_file,\n",
    ")\n"
   ]
  },
  {
   "source": [
    "Lets have a look at the output annotated tsv:"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "  #CHROM       POS   ID REF ALT  DeepSEA/variantEffects/8988T_DNase_None/diff  \\\n0  chr22  41320486    4   G   T                                     -0.001468   \n1  chr22  31009031    9   T   G                                     -0.038191   \n2  chr22  43024150   15   C   G                                      0.013784   \n3  chr22  43027392   16   A   G                                     -0.060475   \n4  chr22  37469571  122   C   T                                     -0.015216   \n\n   DeepSEA/variantEffects/AoSMC_DNase_None/diff  \\\n0                                      0.001205   \n1                                     -0.019323   \n2                                      0.001041   \n3                                     -0.186859   \n4                                      0.012377   \n\n   DeepSEA/variantEffects/Chorion_DNase_None/diff  \\\n0                                       -0.001497   \n1                                       -0.009417   \n2                                        0.007227   \n3                                       -0.087627   \n4                                       -0.003255   \n\n   DeepSEA/variantEffects/CLL_DNase_None/diff  \\\n0                                   -0.000233   \n1                                   -0.019266   \n2                                    0.000264   \n3                                   -0.007143   \n4                                   -0.006013   \n\n   DeepSEA/variantEffects/Fibrobl_DNase_None/diff  \n0                                        0.001492  \n1                                       -0.026329  \n2                                        0.005147  \n3                                       -0.088509  \n4                                        0.005796  \n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "output_dataframe = pd.read_csv(\"output.tsv\", sep='\\t')\n",
    "print(output_dataframe.iloc[: 5, : 10])"
   ]
  },
  {
   "source": [
    "## Interval based effect prediction\n",
    "\n",
    "It is also possible to extend this to models that already perform variant effect prediction using the same api. An example is - models under group MMSplice. We make use of the kipoi dataloader directly as it already takes a vcf, a fasta and a gtf file. After that, we extract variant information directly from the metadata of the dataloader output to annotate the output predictions with.\n",
    "This ipython notebook assumes that it is executed in an environment in which kipoi-veff2 is installed. For more information check https://github.com/kipoi/kipoi-veff2#install-the-conda-environment"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "from kipoi_veff2 import interval_based\n",
    "\n",
    "vcf_file = \"example_data/test.vcf\"\n",
    "fasta_file = \"example_data/test.fa\"\n",
    "gtf_file = \"example_data/test.gtf\"\n",
    "output_file = \"output_MMSplice_mtsplice.tsv\"\n",
    "model_name = \"MMSplice/mtsplice\"\n",
    "\n",
    "model_config = interval_based.INTERVAL_BASED_MODEL_CONFIGS[model_name]\n",
    "interval_based.score_variants(\n",
    "        model_config=model_config,\n",
    "        vcf_file=vcf_file,\n",
    "        fasta_file=fasta_file,\n",
    "        gtf_file=gtf_file,\n",
    "        output_file=output_file,\n",
    "    )\n"
   ]
  },
  {
   "source": [
    "Lets have a look at the annotated tsv"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "   #CHROM       POS                                       ID        REF ALT  \\\n0      17  41197805  17:41197805:ACATCTGCC>A:ENSE00001814242  ACATCTGCC   A   \n1      17  41197805  17:41197805:ACATCTGCC>A:ENSE00001312675  ACATCTGCC   A   \n2      17  41197805  17:41197805:ACATCTGCC>A:ENSE00001831829  ACATCTGCC   A   \n3      17  41197805  17:41197805:ACATCTGCC>A:ENSE00002914501  ACATCTGCC   A   \n4      17  41197805  17:41197805:ACATCTGCC>A:ENSE00001937547  ACATCTGCC   A   \n\n   MMSplice/mtsplice/Retina_Eye  MMSplice/mtsplice/RPE_Choroid_Sclera_Eye  \\\n0                      0.245053                                 -0.105601   \n1                      0.244987                                 -0.106446   \n2                      0.269683                                 -0.102571   \n3                      0.398744                                 -0.087164   \n4                      0.106025                                 -0.007285   \n\n   MMSplice/mtsplice/Subcutaneous_Adipose  \\\n0                               -0.070435   \n1                               -0.070364   \n2                               -0.061823   \n3                               -0.080671   \n4                                0.008817   \n\n   MMSplice/mtsplice/Visceral_Omentum_Adipose  MMSplice/mtsplice/Adrenal_Gland  \n0                                   -0.082142                        -0.066790  \n1                                   -0.081946                        -0.066908  \n2                                   -0.077245                        -0.071069  \n3                                   -0.085462                        -0.057516  \n4                                   -0.004503                         0.055896  \n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "output_dataframe = pd.read_csv(\"output_MMSplice_mtsplice.tsv\", sep='\\t')\n",
    "print(output_dataframe.iloc[: 5, : 10])"
   ]
  },
  {
   "source": [
    "## Variant effect prediction at scale\n",
    "\n",
    "A typical usecase of variant effect prediction pipeline is to predict across many models and many vcf-s or vcf/fasta pairs. To make use of high performance clusters and score variants at scale we use Snakemake. Below is an example Snakemake for running kipoi-veff2 across 24 vcf/fasta pairs from 1000genome project. VCF files are available here - http://hgdownload.cse.ucsc.edu/gbdb/hg19/1000Genomes/phase3/. We assume snakemake is run from a conda environment where kipoi-veff2 is installed. In this specific case, the total number of submitted jobs will be 24*(600+102) = 16848 - 600 models under pwm_HOCOMOCO model group and 102 under DeepBind/Homo_sapiens/RBP for each of 24 vcf/fasta pairs."
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import glob\n",
    "\n",
    "data_dir = \"/data/1000genomevcfs/\"\n",
    "\n",
    "def get_vcf_fasta_pair():\n",
    "    vcf_files = [f'{p.stem.replace(\".vcf\", \"\")}' for p in Path(data_dir).iterdir() if p.is_file() and all(ext in p.suffixes for ext in ['.vcf', '.gz'])]\n",
    "    vcf_fasta_pair = {vf: f'chr{vf.split(\".\")[1][3:]}_standardized.fa' for vf in vcf_files}\n",
    "    return vcf_fasta_pair\n",
    "\n",
    "def get_args(params):\n",
    "    \"\"\"Function returning appropriate parameters with the flag\n",
    "    for the corresponding model\n",
    "    \"\"\"\n",
    "    if \"MMSplice\" in params.model:\n",
    "        return \"-g input/test.gtf\"\n",
    "    else:\n",
    "        return \"-s diff\"\n",
    "\n",
    "groups = [ \"pwm_HOCOMOCO\", \"DeepBind/Homo_sapiens/RBP\"]\n",
    "\n",
    "def get_list_of_models():\n",
    "    \"\"\"Function returning list of models\n",
    "    that belongs to a list of model groups\"\"\"\n",
    "    from kipoi import list_models\n",
    "    all_models = list_models().model\n",
    "    group_to_models = {group : sorted(list(all_models[all_models.str.contains(group)])) for group in groups}\n",
    "    return group_to_models\n",
    "\n",
    "group_to_models = get_list_of_models()\n",
    "vcf_fasta_pair = get_vcf_fasta_pair()\n",
    "\n",
    "rule all:\n",
    "    input: \n",
    "        expand(\"merged__{group}__{id}.tsv\", group=groups, id=vcf_fasta_pair.keys())\n",
    "\n",
    "rule run_vep:\n",
    "    input: \n",
    "        vcf = data_dir+\"{id}.vcf.gz\",\n",
    "        fasta = lambda wildcards: data_dir+vcf_fasta_pair[wildcards.id], \n",
    "    output: \n",
    "        \"output__{model}__{id}.tsv\"\n",
    "    params: \n",
    "        model_args = get_args\n",
    "    shell: \n",
    "        \"kipoi_veff2_predict {input.vcf} {input.fasta} {params.model_args} {output} -m {wildcards.model}\"\n",
    "\n",
    "rule merge_per_group:\n",
    "    input:\n",
    "        lambda wildcards: expand(\"output__{model}__{{id}}.tsv\", model=group_to_models[wildcards.group])\n",
    "    output:\n",
    "        \"merged__{group}__{id}.tsv\"\n",
    "    resources:\n",
    "        mem_mb=28000\n",
    "    shell:\n",
    "        \"kipoi_veff2_merge {input} {output}\""
   ]
  },
  {
   "source": [
    "I highly recommend making a profile for your cluster type using https://github.com/Snakemake-Profiles/. Run the above Snakefile with \n",
    "```snakemake --profile <name_of_your_profile>.```"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "For more information and examples please check https://github.com/kipoi/kipoi-veff2"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "kipoi-veff2",
   "display_name": "kipoi-veff2",
   "language": "python"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}