src/SpliceAI_Predict_Skipping.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/gl/g6/wlouie/anaconda3/envs/skipguide_data_processing/lib/python3.7/site-packages/ipykernel_launcher.py:28: MatplotlibDeprecationWarning: \n",
      "The examples.directory rcparam was deprecated in Matplotlib 3.0 and will be removed in 3.2. In the future, examples will be found relative to the 'datapath' directory.\n",
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "from config import *\n",
    "from utils import *\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import regex\n",
    "import copy\n",
    "import numpy as np\n",
    "import collections\n",
    "import multiprocessing\n",
    "import pickle\n",
    "\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "# Suppress pandas future warning, which messes tqdm\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "import pandas as pd\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "inline_rc = dict(mpl.rcParams)\n",
    "import subprocess\n",
    "\n",
    "import time\n",
    "\n",
    "from keras.models import load_model\n",
    "from pkg_resources import resource_filename\n",
    "from spliceai.utils import one_hot_encode\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SpliceAI Predict\n",
    "SpliceAI can predict the probability of a base being part of an acceptor given its surrounding sequence context."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SpliceAI Predict dat-B\n",
    "For every repair product observed in dat-B, use SpliceAI to detect any acceptor regions, by taking max(SpliceAI(base) -> P(acceptor) for each base). Save the scores for later analysis in MetaSplice_SkipGuide_Evaluation.ipynb."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "indel_splice_precas_count_map = load_bc_seq(INDEL_SPLICE_PRECAS_COUNT_MAP)\n",
    "indel_splice_postcas_count_map = load_bc_seq(INDEL_SPLICE_POSTCAS_COUNT_MAP)\n",
    "gt_splice_count_map = load_bc_seq(GT_SPLICE_COUNT_MAP)\n",
    "gt_precas_splice_count_map = load_bc_seq(GT_PRECAS_SPLICE_COUNT_MAP)\n",
    "predicted_gt_indel_dist_map = load_var(PREDICTED_GT_INDEL_DIST_MAP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "EXON_A = 'CAAGATCCGCCACAACATCGAG'\n",
    "INTRON_TARGETSTART = 'GTAAGTTATCACCTTCGTGGCTACAGAGTTTCCTTATTTGTCTCTGTTGCCGGCTTATATGGACAAGCATATCACAGCCATTTATCGGAGCGCCTCCGTACACGCTATTATCGGACGCCTCGCGAGATCAATACGATTACCAGCTGCCCTCGTCGAC'\n",
    "TARGETEND_EXON_B = 'TGATTACACATATAGACACGCGAGCAGCCATCTTTTATAGAATGGGTAGAACCCGTCCTAAGGACTCAGATTGAGCATCGTTTGCTTCTCGAGTACTACCTGGTACAGATGTCTCTTCAAACAG'\n",
    "EXON_C_BCSTART = 'GACGGCAGCGTGCAGCTCGCC'\n",
    "BCEND_EXON_C = 'GACCACTACCAGCAGAACACCCC'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /data/gl/g6/wlouie/anaconda3/envs/skipguide_data_processing/lib/python3.7/site-packages/tensorflow/python/framework/op_def_library.py:263: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Colocations handled automatically by placer.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/gl/g6/wlouie/anaconda3/envs/skipguide_data_processing/lib/python3.7/site-packages/keras/engine/saving.py:292: UserWarning: No training configuration found in save file: the model was *not* compiled. Compile it manually.\n",
      "  warnings.warn('No training configuration found in save file: '\n"
     ]
    }
   ],
   "source": [
    "# Load spliceai models\n",
    "context = 10000\n",
    "paths = ('models/spliceai{}.h5'.format(x) for x in range(1, 6))\n",
    "models = [load_model(resource_filename('spliceai', x)) for x in paths]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spliceai_predict(payload):\n",
    "    i, model, x = payload\n",
    "    with tf.device('/device:GPU:' + str(i)):\n",
    "        return model.predict(x, verbose=verbose)\n",
    "    \n",
    "\n",
    "def predict_SpliceAIScores(indel_splice_count_map, verbose=1):\n",
    "    # target len 61, plus 1 possible insertion = 62 max\n",
    "    max_seq_len = context + len(INTRON_TARGETSTART) + 62 + len(TARGETEND_EXON_B)\n",
    "    x = np.zeros(\n",
    "        (sum(indel[1] == 'N' or indel[1] in DELETION_SIGNATURES or (indel[1] in INSERTION_SIGNATURES and indel[2] == 1) for i, (indel, splice_count_map) in enumerate(indel_iterator(indel_splice_count_map))),\n",
    "        max_seq_len,\n",
    "        4))\n",
    "    \n",
    "    for i, (indel, splice_count_map) in enumerate(indel_iterator(indel_splice_count_map)):\n",
    "        target = exp_tid_target_map[exp_gid_tid_map[splice_count_map['gid']]]\n",
    "        if indel[1] == 'N':\n",
    "            s = target\n",
    "        elif indel[1] in DELETION_SIGNATURES:\n",
    "            s = get_simulated_product(indel, target)\n",
    "        elif indel[1] in INSERTION_SIGNATURES and indel[2] == 1:\n",
    "            s = get_simulated_product(indel, target)\n",
    "        else:\n",
    "            continue\n",
    "        \n",
    "        input_sequence = 'N'*(context//2) + INTRON_TARGETSTART + s + TARGETEND_EXON_B\n",
    "        input_sequence += 'N'*(max_seq_len - len(input_sequence))\n",
    "        x[i, :, :] = one_hot_encode(input_sequence)[None, :]\n",
    "\n",
    "    y = np.mean([models[m].predict(x, verbose=verbose) for m in range(5)], axis=0)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_acceptor_score(preds):\n",
    "    return np.max(preds[:, :, 1], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(SPLICEAI_PRECAS_PREDS_PATH):\n",
    "    spliceai_precas_preds = predict_SpliceAIScores(indel_splice_precas_count_map)\n",
    "    with open(SPLICEAI_PRECAS_PREDS_PATH, 'wb') as f:\n",
    "        np.save(f, spliceai_precas_preds)\n",
    "else:\n",
    "    with open(SPLICEAI_PRECAS_PREDS_PATH, 'rb') as f:\n",
    "        spliceai_precas_preds = np.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(SPLICEAI_POSTCAS_PREDS_PATH):\n",
    "    spliceai_postcas_preds = predict_SpliceAIScores(indel_splice_postcas_count_map)\n",
    "    with open(SPLICEAI_POSTCAS_PREDS_PATH, 'wb') as f:\n",
    "        np.save(f, spliceai_postcas_preds)\n",
    "else:\n",
    "    with open(SPLICEAI_POSTCAS_PREDS_PATH, 'rb') as f:\n",
    "        spliceai_postcas_preds = np.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SpliceAI Predict on inDelphi Predictions\n",
    "For every designed guide/target pair, we used inDelphi to predict repair genotypes and their frequencies. For each predicted repair product, use SpliceAI to detect any acceptor regions, by taking max(SpliceAI(base) -> P(acceptor) for each base). Save the scores for later analysis in MetaSplice_SkipGuide_Evaluation.ipynb."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_spliceai_predict_gt_indel_dist(pairs):\n",
    "    gid_spliceai_repair_preds = {}\n",
    "    for i, p in enumerate(tqdm(pairs)):\n",
    "        gid = exp_grna_gid_map[p[0]][0]\n",
    "        cutsite = get_cutsite(*p)\n",
    "        indel_splice_count_map = collections.defaultdict(lambda: collections.defaultdict(int))\n",
    "        distribution = predicted_gt_indel_dist_map[p]\n",
    "        for deletion_size in range(1, MAX_INDEL_LEN + 1):\n",
    "            for genotype_pos in distribution[-deletion_size]:\n",
    "                if distribution[-deletion_size][genotype_pos] > 0:\n",
    "                    indel = ('', 'DS', deletion_size, genotype_pos, cutsite)\n",
    "                    indel_splice_count_map[indel]['gid'] = gid\n",
    "                    indel_splice_count_map[indel]['indelphifreq'] = distribution[-deletion_size][genotype_pos]\n",
    "        for base in 'AGTC':\n",
    "            if distribution[1][base] > 0:\n",
    "                indel = ('', 'IS', 1, base, cutsite)\n",
    "                indel_splice_count_map[indel]['gid'] = gid\n",
    "                indel_splice_count_map[indel]['indelphifreq'] = distribution[1][base]\n",
    "                \n",
    "        preds = predict_SpliceAIScores(indel_splice_count_map, verbose=0)\n",
    "        y_preds = get_acceptor_score(preds)\n",
    "        indelphifreq = [indel_splice_count_map[indel]['indelphifreq'] for indel, _ in indel_iterator(indel_splice_count_map)]\n",
    "        \n",
    "        gid_spliceai_repair_preds[gid] = [y_preds.tolist(), indelphifreq]\n",
    "        \n",
    "    save_var(gid_spliceai_repair_preds, SPLICEAI_GT_REPAIR_PREDS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not pickle_exists(SPLICEAI_GT_REPAIR_PREDS):\n",
    "    lib_pairs = [gid_to_gt(gid) for gid in exp_gid_tid_map]\n",
    "    save_spliceai_predict_gt_indel_dist(lib_pairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "gid_spliceai_repair_preds = load_var(SPLICEAI_GT_REPAIR_PREDS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:skipguide_data_processing] *",
   "language": "python",
   "name": "conda-env-skipguide_data_processing-py"
  },
  "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}