tutorials/sat_mut.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Saturation mutagenesis is a powerful tool both for dissecting a specific sequence of interest and understanding what the model learned. [basenji_sat.py](https://github.com/calico/basenji/blob/master/bin/basenji_sat.py) enables this analysis from a test set of data. [basenji_sat_vcf.py](https://github.com/calico/basenji/blob/master/bin/basenji_sat_vcf.py) lets you provide a VCF file for variant-centered mutagenesis.\n",
    "\n",
    "To do this, you'll need\n",
    " * Trained model\n",
    " * Input file (FASTA or HDF5 with test_in/test_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, you can either train your own model in the [Train/test tutorial](https://github.com/calico/basenji/blob/master/tutorials/train_test.ipynb) or use one that I pre-trained from the models subdirectory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll bash the PIM1 promoter to see what motifs drive its expression. I placed a 262 kb FASTA file surrounding the PIM1 TSS in data/pim1.fa, so we'll use [basenji_sat.py](https://github.com/calico/basenji/blob/master/bin/basenji_sat.py).\n",
    "\n",
    "The most relevant options are:\n",
    "\n",
    "| Option/Argument | Value | Note |\n",
    "|:---|:---|:---|\n",
    "| -g | | Plot the nucleotides proportional to the gain score, too. |\n",
    "| -f | 20 | Figure width, that I usually scale to 10x the saturation mutageneis region |\n",
    "| -l | 200 | Saturation mutagenesis region in the center of the given sequence(s) |\n",
    "| -o | pim1_sat | Outplot plot directory. |\n",
    "| --rc | | Predict forward and reverse complement versions and average the results. |\n",
    "| -t | 0,38 | Target indexes. 0 is a DNase and 38 is CAGE, as you can see in data/gm12878_wigs.txt. |\n",
    "| params_file | models/params_small_sat.txt | Table of parameters to setup the model architecture and optimization parameters. |\n",
    "| model_file | models/gm12878.tf | Trained saved model prefix. |\n",
    "| input_file | data/pim1.fa | Either FASTA or HDF5 with test_in/test_out keys. |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cnn_pool': [1, 2, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1], 'batch_size': 1, 'adam_beta1': 0.97, 'cnn_dilation': [1, 1, 1, 1, 1, 1, 2, 4, 8, 16, 32, 64, 128, 1], 'adam_beta2': 0.98, 'loss': 'poisson', 'num_targets': 39, 'cnn_filters': [196, 196, 235, 282, 338, 384, 64, 64, 64, 64, 64, 64, 64, 512], 'cnn_dropout': [0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.1], 'link': 'softplus', 'target_pool': 128, 'cnn_dense': [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0], 'batch_buffer': 16384, 'batch_renorm': 1, 'cnn_filter_sizes': [22, 1, 6, 6, 6, 3, 3, 3, 3, 3, 3, 3, 3, 3], 'learning_rate': 0.002}\n",
      "Targets pooled by 128 to length 2048\n",
      "Convolution w/ 196 4x22 filters strided 1, dilated 1\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 196 196x1 filters strided 1, dilated 1\n",
      "Batch normalization\n",
      "ReLU\n",
      "Max pool 2\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 235 196x6 filters strided 1, dilated 1\n",
      "Batch normalization\n",
      "ReLU\n",
      "Max pool 4\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 282 235x6 filters strided 1, dilated 1\n",
      "Batch normalization\n",
      "ReLU\n",
      "Max pool 4\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 338 282x6 filters strided 1, dilated 1\n",
      "Batch normalization\n",
      "ReLU\n",
      "Max pool 4\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 384 338x3 filters strided 1, dilated 1\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 64 384x3 filters strided 1, dilated 2\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 64 448x3 filters strided 1, dilated 4\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 64 512x3 filters strided 1, dilated 8\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 64 576x3 filters strided 1, dilated 16\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 64 640x3 filters strided 1, dilated 32\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 64 704x3 filters strided 1, dilated 64\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 64 768x3 filters strided 1, dilated 128\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.050\n",
      "Convolution w/ 512 832x3 filters strided 1, dilated 1\n",
      "Batch normalization\n",
      "ReLU\n",
      "Dropout w/ probability 0.100\n",
      "Convolution w/ 39 512x1 filters to final targets\n",
      "Model building time 9.124612\n",
      "2017-09-16 09:51:15.501533: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.2 instructions, but these are available on your machine and could speed up CPU computations.\n",
      "2017-09-16 09:51:15.501556: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use AVX instructions, but these are available on your machine and could speed up CPU computations.\n",
      "2017-09-16 09:51:15.501572: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use AVX2 instructions, but these are available on your machine and could speed up CPU computations.\n",
      "2017-09-16 09:51:15.501576: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use FMA instructions, but these are available on your machine and could speed up CPU computations.\n",
      "Mutating sequence 1 / 1\n",
      "/Users/davidkelley/anaconda3/lib/python3.5/site-packages/matplotlib/tight_layout.py:226: UserWarning: tight_layout : falling back to Agg renderer\n",
      "  warnings.warn(\"tight_layout : falling back to Agg renderer\")\n"
     ]
    }
   ],
   "source": [
    "! basenji_sat.py -g -f 20 -l 200 -o pim1_sat --rc -t 0,38 models/params_med.txt models/gm12878.tf data/pim1.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The saturated mutagenesis heatmaps go into pim1_sat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First the DNASE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"1200\"\n",
       "            height=\"400\"\n",
       "            src=\"pim1_sat/seq0_t0.pdf\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x106ba39b0>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import IFrame\n",
    "IFrame('pim1_sat/seq0_t0.pdf', width=1200, height=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second the CAGE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"1200\"\n",
       "            height=\"400\"\n",
       "            src=\"pim1_sat/seq0_t1.pdf\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x106ba31d0>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "IFrame('pim1_sat/seq0_t1.pdf', width=1200, height=400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}