tutorials/sad.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import os\n",
    "import subprocess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Precursors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isfile('data/hg19.ml.fa'):\n",
    "    subprocess.call('curl -o data/hg19.ml.fa https://storage.googleapis.com/basenji_tutorial_data/hg19.ml.fa', shell=True)\n",
    "    subprocess.call('curl -o data/hg19.ml.fa.fai https://storage.googleapis.com/basenji_tutorial_data/hg19.ml.fa.fai', shell=True)                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isdir('models/heart'):\n",
    "    os.mkdir('models/heart')\n",
    "if not os.path.isfile('models/heart/model_best.h5'):\n",
    "    subprocess.call('curl -o models/heart/model_best.h5 https://storage.googleapis.com/basenji_tutorial_data/model_best.h5', shell=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "lines = [['index','identifier','file','clip','sum_stat','description']]\n",
    "lines.append(['0', 'CNhs11760', 'data/CNhs11760.bw', '384', 'sum', 'aorta'])\n",
    "lines.append(['1', 'CNhs12843', 'data/CNhs12843.bw', '384', 'sum', 'artery'])\n",
    "lines.append(['2', 'CNhs12856', 'data/CNhs12856.bw', '384', 'sum', 'pulmonic_valve'])\n",
    "\n",
    "samples_out = open('data/heart_wigs.txt', 'w')\n",
    "for line in lines:\n",
    "    print('\\t'.join(line), file=samples_out)\n",
    "samples_out.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SNP activity difference compute"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analyzing noncoding variation associated with disease is a major application of Basenji. I now offer several tools to enable that analysis. If you have a small set of variants and know what datasets are most relevant, [basenji_sat_vcf.py](https://github.com/calico/basenji/blob/master/bin/basenji_sat_vcf.py) lets you perform a saturation mutagenesis of the variant and surrounding region to see the relevant nearby motifs.\n",
    "\n",
    "If you want scores measuring the influence of those variants on all datasets,\n",
    " * [basenji_sad.py](https://github.com/calico/basenji/blob/master/bin/basenji_sad.py) computes my SNP activity difference (SAD) score--the predicted change in aligned fragments to the region.\n",
    " * [basenji_sed.py](https://github.com/calico/basenji/blob/master/bin/basenji_sed.py) computes my SNP expression difference (SED) score--the predicted change in aligned fragments to gene TSS's.\n",
    "\n",
    "Here, I'll demonstrate those two programs. 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": [
    "As an example, we'll study a prostate cancer susceptibility allele of rs339331 that increases RFX6 expression by modulating HOXB13 chromatin binding (http://www.nature.com/ng/journal/v46/n2/full/ng.2862.html).\n",
    "\n",
    "First, we'll use [basenji_sad.py](https://github.com/calico/basenji/blob/master/bin/basenji_sad.py) to predict across the region for each allele and compute stats about the mean and max differences.\n",
    "\n",
    "The most relevant options are:\n",
    "\n",
    "| Option/Argument | Value | Note |\n",
    "|:---|:---|:---|\n",
    "| -f | data/hg19.ml.fa | Genome fasta. |\n",
    "| -g | data/human.hg19.genome | Genome assembly chromosome length to bound gene sequences. |\n",
    "| -o | rfx6_sad | Outplot plot directory. |\n",
    "| --rc | True | Ensemble predictions for forward and reverse complement sequences. |\n",
    "| --shift | 1,0,-1 | Ensemble predictions for sequences shifted by 1, 0, and -1 bp. |\n",
    "| -t | data/heart_wigs.txt | Target labels. |\n",
    "| params_file | models/params_small.json | JSON specified parameters to setup the model architecture and optimization. |\n",
    "| model_file | models/heart/model_best.h5 | Trained saved model parameters. |\n",
    "| vcf_file | data/rs339331.vcf | VCF file specifying variants to score. |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2021-02-12 15:49:10.346186: I tensorflow/compiler/jit/xla_cpu_device.cc:41] Not creating XLA devices, tf_xla_enable_xla_devices not set\n",
      "2021-02-12 15:49:10.346327: I tensorflow/core/platform/cpu_feature_guard.cc:142] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n",
      "Model: \"model_1\"\n",
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "sequence (InputLayer)           [(None, 131072, 4)]  0                                            \n",
      "__________________________________________________________________________________________________\n",
      "stochastic_reverse_complement ( ((None, 131072, 4),  0           sequence[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "stochastic_shift (StochasticShi (None, 131072, 4)    0           stochastic_reverse_complement[0][\n",
      "__________________________________________________________________________________________________\n",
      "gelu (GELU)                     (None, 131072, 4)    0           stochastic_shift[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv1d (Conv1D)                 (None, 131072, 64)   3840        gelu[0][0]                       \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization (BatchNorma (None, 131072, 64)   256         conv1d[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "max_pooling1d (MaxPooling1D)    (None, 16384, 64)    0           batch_normalization[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "gelu_1 (GELU)                   (None, 16384, 64)    0           max_pooling1d[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_1 (Conv1D)               (None, 16384, 64)    20480       gelu_1[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_1 (BatchNor (None, 16384, 64)    256         conv1d_1[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "max_pooling1d_1 (MaxPooling1D)  (None, 4096, 64)     0           batch_normalization_1[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "gelu_2 (GELU)                   (None, 4096, 64)     0           max_pooling1d_1[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_2 (Conv1D)               (None, 4096, 72)     23040       gelu_2[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_2 (BatchNor (None, 4096, 72)     288         conv1d_2[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "max_pooling1d_2 (MaxPooling1D)  (None, 1024, 72)     0           batch_normalization_2[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "gelu_3 (GELU)                   (None, 1024, 72)     0           max_pooling1d_2[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_3 (Conv1D)               (None, 1024, 32)     6912        gelu_3[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_3 (BatchNor (None, 1024, 32)     128         conv1d_3[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "gelu_4 (GELU)                   (None, 1024, 32)     0           batch_normalization_3[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_4 (Conv1D)               (None, 1024, 72)     2304        gelu_4[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_4 (BatchNor (None, 1024, 72)     288         conv1d_4[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "dropout (Dropout)               (None, 1024, 72)     0           batch_normalization_4[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "add (Add)                       (None, 1024, 72)     0           max_pooling1d_2[0][0]            \n",
      "                                                                 dropout[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "gelu_5 (GELU)                   (None, 1024, 72)     0           add[0][0]                        \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_5 (Conv1D)               (None, 1024, 32)     6912        gelu_5[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_5 (BatchNor (None, 1024, 32)     128         conv1d_5[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "gelu_6 (GELU)                   (None, 1024, 32)     0           batch_normalization_5[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_6 (Conv1D)               (None, 1024, 72)     2304        gelu_6[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_6 (BatchNor (None, 1024, 72)     288         conv1d_6[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "dropout_1 (Dropout)             (None, 1024, 72)     0           batch_normalization_6[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "add_1 (Add)                     (None, 1024, 72)     0           add[0][0]                        \n",
      "                                                                 dropout_1[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "gelu_7 (GELU)                   (None, 1024, 72)     0           add_1[0][0]                      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_7 (Conv1D)               (None, 1024, 32)     6912        gelu_7[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_7 (BatchNor (None, 1024, 32)     128         conv1d_7[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "gelu_8 (GELU)                   (None, 1024, 32)     0           batch_normalization_7[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_8 (Conv1D)               (None, 1024, 72)     2304        gelu_8[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_8 (BatchNor (None, 1024, 72)     288         conv1d_8[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "dropout_2 (Dropout)             (None, 1024, 72)     0           batch_normalization_8[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "add_2 (Add)                     (None, 1024, 72)     0           add_1[0][0]                      \n",
      "                                                                 dropout_2[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "gelu_9 (GELU)                   (None, 1024, 72)     0           add_2[0][0]                      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_9 (Conv1D)               (None, 1024, 32)     6912        gelu_9[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_9 (BatchNor (None, 1024, 32)     128         conv1d_9[0][0]                   \n",
      "__________________________________________________________________________________________________\n",
      "gelu_10 (GELU)                  (None, 1024, 32)     0           batch_normalization_9[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_10 (Conv1D)              (None, 1024, 72)     2304        gelu_10[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_10 (BatchNo (None, 1024, 72)     288         conv1d_10[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "dropout_3 (Dropout)             (None, 1024, 72)     0           batch_normalization_10[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "add_3 (Add)                     (None, 1024, 72)     0           add_2[0][0]                      \n",
      "                                                                 dropout_3[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "gelu_11 (GELU)                  (None, 1024, 72)     0           add_3[0][0]                      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_11 (Conv1D)              (None, 1024, 32)     6912        gelu_11[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_11 (BatchNo (None, 1024, 32)     128         conv1d_11[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "gelu_12 (GELU)                  (None, 1024, 32)     0           batch_normalization_11[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_12 (Conv1D)              (None, 1024, 72)     2304        gelu_12[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_12 (BatchNo (None, 1024, 72)     288         conv1d_12[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "dropout_4 (Dropout)             (None, 1024, 72)     0           batch_normalization_12[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "add_4 (Add)                     (None, 1024, 72)     0           add_3[0][0]                      \n",
      "                                                                 dropout_4[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "gelu_13 (GELU)                  (None, 1024, 72)     0           add_4[0][0]                      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_13 (Conv1D)              (None, 1024, 32)     6912        gelu_13[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_13 (BatchNo (None, 1024, 32)     128         conv1d_13[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "gelu_14 (GELU)                  (None, 1024, 32)     0           batch_normalization_13[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_14 (Conv1D)              (None, 1024, 72)     2304        gelu_14[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_14 (BatchNo (None, 1024, 72)     288         conv1d_14[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "dropout_5 (Dropout)             (None, 1024, 72)     0           batch_normalization_14[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "add_5 (Add)                     (None, 1024, 72)     0           add_4[0][0]                      \n",
      "                                                                 dropout_5[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "gelu_15 (GELU)                  (None, 1024, 72)     0           add_5[0][0]                      \n",
      "__________________________________________________________________________________________________\n",
      "conv1d_15 (Conv1D)              (None, 1024, 64)     4608        gelu_15[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_15 (BatchNo (None, 1024, 64)     256         conv1d_15[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "dropout_6 (Dropout)             (None, 1024, 64)     0           batch_normalization_15[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "gelu_16 (GELU)                  (None, 1024, 64)     0           dropout_6[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "dense (Dense)                   (None, 1024, 3)      195         gelu_16[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "switch_reverse (SwitchReverse)  (None, 1024, 3)      0           dense[0][0]                      \n",
      "                                                                 stochastic_reverse_complement[0][\n",
      "==================================================================================================\n",
      "Total params: 111,011\n",
      "Trainable params: 109,235\n",
      "Non-trainable params: 1,776\n",
      "__________________________________________________________________________________________________\n",
      "None\n",
      "model_strides [128]\n",
      "target_lengths [1024]\n",
      "target_crops [0]\n",
      "2021-02-12 15:49:14.227988: I tensorflow/compiler/mlir/mlir_graph_optimization_pass.cc:116] None of the MLIR optimization passes are enabled (registered 2)\n"
     ]
    }
   ],
   "source": [
    "! basenji_sad.py -f data/hg19.ml.fa -o output/rfx6_sad --rc --shift \"1,0,-1\" -t data/heart_wigs.txt models/params_small.json models/heart/model_best.h5 data/rs339331.vcf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SNP activity difference output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output HDF5 stores the SNP and target information and predicted scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['SAD',\n",
       " 'SAD_pct',\n",
       " 'alt_allele',\n",
       " 'chr',\n",
       " 'percentiles',\n",
       " 'pos',\n",
       " 'ref_allele',\n",
       " 'snp',\n",
       " 'target_ids',\n",
       " 'target_labels']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sad_h5 = h5py.File('output/rfx6_sad/sad.h5', 'r')\n",
    "list(sad_h5.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "snp [b'rs339331']\n",
      "chr [b'chr6']\n",
      "pos [117210052]\n",
      "ref_allele [b'T']\n"
     ]
    }
   ],
   "source": [
    "for snp_key in ['snp', 'chr', 'pos', 'ref_allele']:\n",
    "    print(snp_key, sad_h5[snp_key][:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 0   0.0018  b'CNhs11760'  b'aorta'\n",
      " 1  -0.0014  b'CNhs12843'  b'artery'\n",
      " 2  -0.0028  b'CNhs12856'  b'pulmonic_valve'\n"
     ]
    }
   ],
   "source": [
    "for ti in range(3):\n",
    "    cols = (ti, sad_h5['SAD'][0,ti], sad_h5['target_ids'][ti], sad_h5['target_labels'][ti])\n",
    "    print('%2d  %7.4f  %12s  %s' % cols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are inconclusive small effect sizes, not surprising given that we're only studying heart CAGE. The proper cell types and experiments would shed more light."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SNP expression difference compute"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can directly query the predictions at gene TSS's. Since moving the codebase to Tensorflow2, these scripts need to be rewritten. In my experience thus far, the more general SNP scores above work well for nearly all applications. The model isn't terribly good at predicting that one gene versus another will be more affected right now."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}