Document/Supplementary_Data_File_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DeepBind integration in Motif-Raptor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "In this notebook we show how to integrate the binding scores of DeepBind (Alipanahi et. Al Nature Biotech 2015: https://www.nature.com/articles/nbt.3300) in Motif-Raptor. We will focus on the step 2 of the pipeline i.e. the ranking of potential causal TFs for each of the prioritized cell types obtained in step 1 based on a set of SNPs in cell type specific DHS.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Requirements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A working installation of Motif-Raptor (see instructions here: https://github.com/pinellolab/MotifRaptor#installation)\n",
    "\n",
    "A working installation of DeepBind (see instructions here: http://tools.genes.toronto.edu/deepbind/)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we will use one DeepBind model for NFKB1 (https://github.com/pinellolab/MotifRaptor#installation) using cell type specific chromatin accessible sites in CD8+ cells and summary statistics for rheumatoid arthritis (Okada et al., 2014). However this preocedure is general and can be extended to all the precomputed models available in DeepBind and other cell type and GWAS studies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The input to this integration is the output from the step 2 from Motif-Raptor and in particular the following files (see the tutorial for more information on how to run step 1 and step 2: https://github.com/pinellolab/MotifRaptor#tutorial ):\n",
    "\n",
    "(1)\"hit_snps_seq_pickle_celltypeid.df.txt\":\n",
    "\n",
    "    This file is from Motif-Raptor step 2 (step2_out). It contains SNP hits defined by the summary statistics and overlapping the DHS of the target cell type. From this file, you can easily get the flanking sequences around each SNP for reference and alternative alleles.\n",
    "    \n",
    "    An example of this file will be used in the step 1 of this notebook.\n",
    "\n",
    "(2) \"motifid.scores\":\n",
    "\n",
    "    This file contains background SNP IDs from Motif-Raptor step 2 (step2_out/motif_result/background_files). It's handy to use this file to generate the background data set. However, you can also define your own background data set as long as it follows the required format as input to Deepbind.\n",
    "    \n",
    "    An example of this file will be used in the step 3 of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Get list of SNPs from Motif-Raptor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This file is for Rheumatoid Arthritis in CD8+ cells, which can be obtained from Motif-Raptor step2 output. It's described in the above introduction section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ID</th>\n",
       "      <th>CHR</th>\n",
       "      <th>POS</th>\n",
       "      <th>REF</th>\n",
       "      <th>ALT</th>\n",
       "      <th>seq_ref</th>\n",
       "      <th>seq_alt</th>\n",
       "      <th>pos_start_ref</th>\n",
       "      <th>pos_end_ref</th>\n",
       "      <th>pos_start_alt</th>\n",
       "      <th>pos_end_alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>rs11585048</td>\n",
       "      <td>1</td>\n",
       "      <td>2534087</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>ACCAGGTCACGTCCTGTCATACCCCCACTCTTCTCGGGCAGAGTCT...</td>\n",
       "      <td>ACCAGGTCACGTCCTGTCATACCCCCACTCCTCTCGGGCAGAGTCT...</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>rs10797435</td>\n",
       "      <td>1</td>\n",
       "      <td>2534801</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>GCAGGGTGAGAGTCGGCTGAGGGAGTGGGGGACTTGCCCCTTCCCA...</td>\n",
       "      <td>GCAGGGTGAGAGTCGGCTGAGGGAGTGGGGCACTTGCCCCTTCCCA...</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>rs4648652</td>\n",
       "      <td>1</td>\n",
       "      <td>2535758</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>AAAGGACGCCGGGCAGAGCTGGGTGGCCTCAGCCACCATGGCCCCC...</td>\n",
       "      <td>AAAGGACGCCGGGCAGAGCTGGGTGGCCTCGGCCACCATGGCCCCC...</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           ID  CHR      POS REF ALT  \\\n",
       "0  rs11585048    1  2534087   T   C   \n",
       "1  rs10797435    1  2534801   G   C   \n",
       "2   rs4648652    1  2535758   A   G   \n",
       "\n",
       "                                             seq_ref  \\\n",
       "0  ACCAGGTCACGTCCTGTCATACCCCCACTCTTCTCGGGCAGAGTCT...   \n",
       "1  GCAGGGTGAGAGTCGGCTGAGGGAGTGGGGGACTTGCCCCTTCCCA...   \n",
       "2  AAAGGACGCCGGGCAGAGCTGGGTGGCCTCAGCCACCATGGCCCCC...   \n",
       "\n",
       "                                             seq_alt  pos_start_ref  \\\n",
       "0  ACCAGGTCACGTCCTGTCATACCCCCACTCCTCTCGGGCAGAGTCT...             30   \n",
       "1  GCAGGGTGAGAGTCGGCTGAGGGAGTGGGGCACTTGCCCCTTCCCA...             30   \n",
       "2  AAAGGACGCCGGGCAGAGCTGGGTGGCCTCGGCCACCATGGCCCCC...             30   \n",
       "\n",
       "   pos_end_ref  pos_start_alt  pos_end_alt  \n",
       "0           30             30           30  \n",
       "1           30             30           30  \n",
       "2           30             30           30  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filename=\"RA/step2_out/hit_snps_seq_pickle_ENCFF512IML.df.txt\"\n",
    "hit_snp_pd=pd.read_csv(filename,sep=\"\\t\",header=0)\n",
    "hit_snp_pd.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We collect flanking sequences for reference and alternative allels around each SNP, and then we write these sequences to a sequence file *.seq, which is required for deepbind program.\n",
    "In the \"demo_RA_CD8.seq\" generated below, we use consecutive two lines to store the sequences (reference and alternative) for each SNP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seq_list=hit_snp_pd['seq_ref'].tolist()+hit_snp_pd['seq_alt'].tolist()\n",
    "outfile=open(\"demo_RA_CD8.seq\",\"w\")\n",
    "for seq in seq_list:\n",
    "    outfile.write(seq+\"\\n\")\n",
    "outfile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Run a deepbind model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "deepbind_id=\"D00540.002\" # This is NFKB1 model from deepbind"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to write deepbind_id (each ID in one line) into a file \"demo_input.ids\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "!touch demo_input.ids\n",
    "!echo \"D00540.002\" > demo_input.ids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the obtained file (\"demo_RA_CD8.seq\"), you can run deepbind using the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!./deepbind demo_input.ids < demo_RA_CD8.seq > demo_deepbind_out_RA_CD8.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Build the null model for a deepbind model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is similar to the previous step:\n",
    "You need to organize the background SNP using the same format as in \"demo_RA_CD8.seq\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also define your own background SNP set based on the similar format as following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Genome</th>\n",
       "      <th>pos_SNP</th>\n",
       "      <th>SNP_in_ref</th>\n",
       "      <th>SNP_in_alt</th>\n",
       "      <th>ID</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>13118</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>rs200579949</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>540426</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>rs540135598</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>568709</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>rs148329687</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Genome  pos_SNP SNP_in_ref SNP_in_alt           ID\n",
       "0   chr1    13118          A          G  rs200579949\n",
       "1   chr1   540426          T          C  rs540135598\n",
       "2   chr1   568709          A          G  rs148329687"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "backgroundfilename=\"RA/step2_out/motif_result/background_files/MA0105.1__NFKB1.scores\"\n",
    "background_df=pd.read_csv(backgroundfilename,sep=\"\\t\",header=0)\n",
    "background_df[['Genome','pos_SNP','SNP_in_ref','SNP_in_alt','ID']].head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the flanking sequences using \"twobitreader\". If you have installed Motif-Raptor, this tool will be already available.\n",
    "\n",
    "You may also use the command line version from UCSC: https://genome.ucsc.edu/goldenPath/help/twoBit.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import twobitreader\n",
    "#change the Database path depending on where you unpack the Database.zip from Motif-Raptor\n",
    "genome = twobitreader.TwoBitFile(\"Database/hg19/hg19.2bit\")\n",
    "\n",
    "def getseq_from_genome(from_genome,chromosomenum, start_pos, end_pos):\n",
    "    chromosome = from_genome[chromosomenum]\n",
    "    seq = chromosome.get_slice(start_pos, end_pos)\n",
    "    return seq.upper()\n",
    "\n",
    "#we only demo a few SNPs from the background, you may change this number\n",
    "testing_num=500\n",
    "testing_df=background_df[['Genome','pos_SNP','SNP_in_ref','SNP_in_alt','ID']].head(testing_num)\n",
    "\n",
    "window_size=30\n",
    "outfile=open(\"demo_RA_CD8_background.seq\",\"w\")\n",
    "for index in testing_df.index:\n",
    "    chrom=testing_df.loc[index]['Genome']\n",
    "    pos=int(testing_df.loc[index]['pos_SNP'])\n",
    "    start_pos=pos-window_size\n",
    "    end_pos=pos+window_size\n",
    "    ref=str(testing_df.loc[index]['SNP_in_ref']).upper()\n",
    "    alt=str(testing_df.loc[index]['SNP_in_alt']).upper()\n",
    "    seq=getseq_from_genome(genome,chrom, start_pos, end_pos)\n",
    "    ref_seq=seq[0:window_size]+ref+seq[window_size+1:2*window_size+1]\n",
    "    alt_seq=seq[0:window_size]+alt+seq[window_size+1:2*window_size+1]\n",
    "    outfile.write(ref_seq+\"\\n\")\n",
    "    outfile.write(ref_seq+\"\\n\")\n",
    "outfile.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the obtained file (\"demo_RA_CD8_background.seq\"), you can run deepbind using the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!./deepbind demo_input.ids < demo_RA_CD8_background.seq > demo_deepbind_out_RA_CD8_background.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Create motif modulation score for deepbind"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function is used to calculate binding scores, disruption scores, and modulation scores for SNPs in \"filename\"\n",
    "\n",
    "[1] filename: deepbind output\n",
    "\n",
    "[2] deepbind_id: deepbind model ID\n",
    "\n",
    "[3] bind_scaler,dirsupt_scaler_pos,disrupt_scaler_neg: these are scalers used to normalize the score. If not set, it will calcualte these scalers on the current data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_SNP_scores(filename, deepbind_id, bind_scaler=0,dirsupt_scaler_pos=0,dirsupt_scaler_neg=0):\n",
    "    bind_scores=list()\n",
    "    disrupt_scores=list()\n",
    "    deepbind_df=pd.read_csv(filename,sep=\"\\t\",header=0)\n",
    "    deepbind_column=deepbind_df[deepbind_id].tolist()\n",
    "    n=int(len(deepbind_column)/2)\n",
    "    \n",
    "    for i in range(0,n):\n",
    "        a=deepbind_column[i]\n",
    "        b=deepbind_column[i+1]\n",
    "        bind_score = max(a,b)\n",
    "        if bind_score>=0:\n",
    "            bind_scores.append(bind_score)\n",
    "            disrupt_score=b-a\n",
    "            disrupt_scores.append(disrupt_score)\n",
    "    if bind_scaler==0:\n",
    "        bind_scaler=max(bind_scores)\n",
    "        dirsupt_scaler_pos=max([i for i in disrupt_scores if i >= 0] or None)\n",
    "        dirsupt_scaler_neg=min([i for i in disrupt_scores if i < 0] or None)\n",
    "    scaled_bind_scores=list(np.divide(np.array(bind_scores),bind_scaler))\n",
    "    scaled_disrupt_scores=list()\n",
    "    for score in disrupt_scores:\n",
    "        if score>=0:\n",
    "            newscore=score/dirsupt_scaler_pos\n",
    "        else:\n",
    "            newscore=-score/dirsupt_scaler_neg\n",
    "        scaled_disrupt_scores.append(newscore)\n",
    "    scaled_mod_scores=list(np.multiply(scaled_bind_scores,scaled_disrupt_scores))\n",
    "    return bind_scaler,dirsupt_scaler_pos,dirsupt_scaler_neg,scaled_bind_scores,scaled_disrupt_scores,scaled_mod_scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate scores for target and background SNPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_filename=\"demo_deepbind_out_RA_CD8.txt\"\n",
    "background_filename=\"demo_deepbind_out_RA_CD8_background.txt\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "bind_scaler,dirsupt_scaler_pos,dirsupt_scaler_neg,scaled_bind_scores_bg,scaled_disrupt_scores_bg,scaled_mod_scores_bg \\\n",
    "    = get_SNP_scores(background_filename,deepbind_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "bind_scaler,dirsupt_scaler_pos,dirsupt_scaler_neg,scaled_bind_scores_target,scaled_disrupt_scores_target,scaled_mod_scores_target \\\n",
    "    = get_SNP_scores(target_filename,deepbind_id,bind_scaler,dirsupt_scaler_pos,dirsupt_scaler_neg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Assess the significance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the target and background modulation scores (i.e. \"scaled_mod_scores_target\" and \"scaled_mod_scores_bg\"), you can calculate p-values for potential disruption, enhancement and dual events for TF binding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Disruption p-value: 0.26782865836842923\n",
      "Enhancement p-value: 0.7321713416315707\n",
      "Dual p-value: 0.03529465255820752\n"
     ]
    }
   ],
   "source": [
    "target_data=scaled_mod_scores_target\n",
    "bg_data=scaled_mod_scores_bg\n",
    "\n",
    "pop_mean=np.mean(bg_data)\n",
    "pop_var=np.var(bg_data)\n",
    "test_mean=pop_mean\n",
    "\n",
    "test_var=pop_var/(len(target_data))\n",
    "\n",
    "\n",
    "t=(np.mean(target_data)-test_mean)/((test_var)**(0.5))\n",
    "report_pvalue_left=norm.cdf(t)\n",
    "report_pvalue_right=1-report_pvalue_left\n",
    "#print(str(test_mean))\n",
    "#print(str(test_var)) \n",
    "abs_target_data=np.abs(target_data)\n",
    "abs_bg_data=np.abs(bg_data)\n",
    "abs_pop_mean=np.mean(abs_bg_data)\n",
    "abs_pop_var=np.var(abs_bg_data)\n",
    "abs_test_mean=abs_pop_mean\n",
    "abs_test_var=abs_pop_var/(len(abs_target_data))\n",
    "if abs_test_var==0:\n",
    "    abs_test_var=1E-20\n",
    "abs_t=(np.mean(abs_target_data)-abs_test_mean)/((abs_test_var)**(0.5))\n",
    "report_pvalue_abs=1-norm.cdf(abs_t)\n",
    "print(\"Disruption p-value: \"+str(report_pvalue_left))\n",
    "print(\"Enhancement p-value: \"+str(report_pvalue_right))\n",
    "print(\"Dual p-value: \"+str(report_pvalue_abs))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "yao_py36",
   "language": "python",
   "name": "yao_py36"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}