{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Variant effect prediction - Part II" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This tutorial illustrate an alternative for performing variant effect prediction using a native keras model (as opposed to using the Janggu model wrapper).\n", "\n", "This tutorial requires janggu>=0.10.0." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "\n", "import numpy as np\n", "import h5py\n", "from keras import backend as K\n", "from keras.layers import Conv2D\n", "from keras.layers import GlobalAveragePooling2D\n", "from pkg_resources import resource_filename\n", "\n", "from janggu import create_model\n", "from janggu import predict_variant_effect\n", "\n", "from janggu import inputlayer\n", "from janggu import outputdense\n", "from janggu.data import Bioseq\n", "from janggu.data import Cover\n", "from janggu.data import GenomicIndexer\n", "from janggu.data import ReduceDim\n", "from janggu.data import plotGenomeTrack\n", "from janggu.layers import DnaConv2D\n", "\n", "np.random.seed(1234)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define and fit a keras model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we need to specify the output directory in which the results are stored and load the datasets. We also specify the number of epochs to train the model and the sequence feature order." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "order = 3\n", "epochs = 100" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "os.environ['JANGGU_OUTPUT'] = '/home/wkopp/janggu_examples'\n", "\n", "# load the dataset\n", "# The pseudo genome represents just a concatenation of all sequences\n", "# in sample.fa and sample2.fa. Therefore, the results should be almost\n", "# identically to the models obtained from classify_fasta.py.\n", "REFGENOME = resource_filename('janggu', 'resources/pseudo_genome.fa')\n", "VCFFILE = resource_filename('janggu', 'resources/pseudo_snps.vcf')\n", "# ROI contains regions spanning positive and negative examples\n", "ROI_TRAIN_FILE = resource_filename('janggu', 'resources/roi_train.bed')\n", "ROI_TEST_FILE = resource_filename('janggu', 'resources/roi_test.bed')\n", "\n", "# PEAK_FILE only contains positive examples\n", "PEAK_FILE = resource_filename('janggu', 'resources/scores.bed')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Training input and labels are purely defined genomic coordinates\n", "DNA = Bioseq.create_from_refgenome('dna', refgenome=REFGENOME,\n", " roi=ROI_TRAIN_FILE,\n", " binsize=200,\n", " order=order)\n", "\n", "LABELS = Cover.create_from_bed('peaks', roi=ROI_TRAIN_FILE,\n", " bedfiles=PEAK_FILE,\n", " binsize=200,\n", " resolution=200)\n", "\n", "\n", "DNA_TEST = Bioseq.create_from_refgenome('dna', refgenome=REFGENOME,\n", " roi=ROI_TEST_FILE,\n", " binsize=200,\n", " order=order)\n", "\n", "LABELS_TEST = Cover.create_from_bed('peaks',\n", " roi=ROI_TEST_FILE,\n", " bedfiles=PEAK_FILE,\n", " binsize=200,\n", " resolution=200)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define and fit a new model" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "WARNING:tensorflow:From /home/wkopp/anaconda3/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:541: The name tf.placeholder is deprecated. Please use tf.compat.v1.placeholder instead.\n", "\n", "WARNING:tensorflow:From /home/wkopp/anaconda3/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:66: The name tf.get_default_graph is deprecated. Please use tf.compat.v1.get_default_graph instead.\n", "\n", "WARNING:tensorflow:From /home/wkopp/anaconda3/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:4432: The name tf.random_uniform is deprecated. Please use tf.random.uniform instead.\n", "\n", "WARNING:tensorflow:From /home/wkopp/anaconda3/lib/python3.7/site-packages/keras/optimizers.py:793: The name tf.train.Optimizer is deprecated. Please use tf.compat.v1.train.Optimizer instead.\n", "\n", "WARNING:tensorflow:From /home/wkopp/anaconda3/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:3657: The name tf.log is deprecated. Please use tf.math.log instead.\n", "\n", "WARNING:tensorflow:From /home/wkopp/anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/nn_impl.py:180: add_dispatch_support..wrapper (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n", "Instructions for updating:\n", "Use tf.where in 2.0, which has the same broadcast rule as np.where\n", "WARNING:tensorflow:From /home/wkopp/anaconda3/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:1033: The name tf.assign_add is deprecated. Please use tf.compat.v1.assign_add instead.\n", "\n", "Epoch 1/100\n", "7797/7797 [==============================] - 12s 1ms/step - loss: 0.5612 - acc: 0.7353\n", "Epoch 2/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.4607 - acc: 0.7981\n", "Epoch 3/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.4115 - acc: 0.8222\n", "Epoch 4/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.3785 - acc: 0.8424\n", "Epoch 5/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.3513 - acc: 0.8552\n", "Epoch 6/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.3255 - acc: 0.8656\n", "Epoch 7/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.3037 - acc: 0.8821\n", "Epoch 8/100\n", "7797/7797 [==============================] - 12s 1ms/step - loss: 0.2801 - acc: 0.8916\n", "Epoch 9/100\n", "7797/7797 [==============================] - 12s 1ms/step - loss: 0.2601 - acc: 0.9006\n", "Epoch 10/100\n", "7797/7797 [==============================] - 12s 1ms/step - loss: 0.2412 - acc: 0.9125\n", "Epoch 11/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.2240 - acc: 0.9207\n", "Epoch 12/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.2083 - acc: 0.9277\n", "Epoch 13/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1951 - acc: 0.9314\n", "Epoch 14/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1822 - acc: 0.9395\n", "Epoch 15/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1712 - acc: 0.9427\n", "Epoch 16/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1621 - acc: 0.9450\n", "Epoch 17/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1535 - acc: 0.9484\n", "Epoch 18/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1448 - acc: 0.9514\n", "Epoch 19/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1387 - acc: 0.9532\n", "Epoch 20/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1311 - acc: 0.9597\n", "Epoch 21/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1249 - acc: 0.9592\n", "Epoch 22/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1196 - acc: 0.9618\n", "Epoch 23/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1136 - acc: 0.9647\n", "Epoch 24/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1091 - acc: 0.9663\n", "Epoch 25/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.1039 - acc: 0.9668\n", "Epoch 26/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0991 - acc: 0.9708\n", "Epoch 27/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0965 - acc: 0.9699\n", "Epoch 28/100\n", "7797/7797 [==============================] - 12s 1ms/step - loss: 0.0921 - acc: 0.9724\n", "Epoch 29/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0887 - acc: 0.9750\n", "Epoch 30/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0860 - acc: 0.9761\n", "Epoch 31/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0826 - acc: 0.9763\n", "Epoch 32/100\n", "7797/7797 [==============================] - 12s 2ms/step - loss: 0.0795 - acc: 0.9774\n", "Epoch 33/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0769 - acc: 0.9810\n", "Epoch 34/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0743 - acc: 0.9794\n", "Epoch 35/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0714 - acc: 0.9811\n", "Epoch 36/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0696 - acc: 0.9819\n", "Epoch 37/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0669 - acc: 0.9837\n", "Epoch 38/100\n", "7797/7797 [==============================] - 14s 2ms/step - loss: 0.0642 - acc: 0.9845\n", "Epoch 39/100\n", "7797/7797 [==============================] - 12s 2ms/step - loss: 0.0620 - acc: 0.9845\n", "Epoch 40/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0605 - acc: 0.9858\n", "Epoch 41/100\n", "7797/7797 [==============================] - 12s 1ms/step - loss: 0.0584 - acc: 0.9858\n", "Epoch 42/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0565 - acc: 0.9867\n", "Epoch 43/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0545 - acc: 0.9876\n", "Epoch 44/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0538 - acc: 0.9882\n", "Epoch 45/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0517 - acc: 0.9895\n", "Epoch 46/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0503 - acc: 0.9887\n", "Epoch 47/100\n", "7797/7797 [==============================] - 12s 1ms/step - loss: 0.0489 - acc: 0.9887\n", "Epoch 48/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0472 - acc: 0.9905\n", "Epoch 49/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0459 - acc: 0.9905\n", "Epoch 50/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0450 - acc: 0.9895\n", "Epoch 51/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0432 - acc: 0.9900\n", "Epoch 52/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0420 - acc: 0.9917\n", "Epoch 53/100\n", "7797/7797 [==============================] - 12s 2ms/step - loss: 0.0403 - acc: 0.9927\n", "Epoch 54/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0394 - acc: 0.9922\n", "Epoch 55/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0383 - acc: 0.9924\n", "Epoch 56/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0368 - acc: 0.9933\n", "Epoch 57/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0361 - acc: 0.9935\n", "Epoch 58/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0352 - acc: 0.9932\n", "Epoch 59/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.0342 - acc: 0.9937\n", "Epoch 60/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0332 - acc: 0.9938\n", "Epoch 61/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0322 - acc: 0.9940\n", "Epoch 62/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.0314 - acc: 0.9944\n", "Epoch 63/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0305 - acc: 0.9944\n", "Epoch 64/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0296 - acc: 0.9958\n", "Epoch 65/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0291 - acc: 0.9959\n", "Epoch 66/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0284 - acc: 0.9955\n", "Epoch 67/100\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0277 - acc: 0.9949\n", "Epoch 68/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.0265 - acc: 0.9958\n", "Epoch 69/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.0260 - acc: 0.9964\n", "Epoch 70/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.0253 - acc: 0.9964\n", "Epoch 71/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.0245 - acc: 0.9967\n", "Epoch 72/100\n", "7797/7797 [==============================] - 9s 1ms/step - loss: 0.0239 - acc: 0.9963\n", "Epoch 73/100\n", "7797/7797 [==============================] - 10s 1ms/step - loss: 0.0230 - acc: 0.9976\n", "Epoch 74/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0228 - acc: 0.9967\n", "Epoch 75/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0218 - acc: 0.9976\n", "Epoch 76/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0216 - acc: 0.9974\n", "Epoch 77/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0208 - acc: 0.9973\n", "Epoch 78/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0201 - acc: 0.9982\n", "Epoch 79/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0197 - acc: 0.9977\n", "Epoch 80/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0192 - acc: 0.9982\n", "Epoch 81/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0187 - acc: 0.9979\n", "Epoch 82/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0183 - acc: 0.9985\n", "Epoch 83/100\n", "7797/7797 [==============================] - 12s 1ms/step - loss: 0.0177 - acc: 0.9986\n", "Epoch 84/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0172 - acc: 0.9987\n", "Epoch 85/100\n", "7797/7797 [==============================] - 13s 2ms/step - loss: 0.0167 - acc: 0.9987\n", "Epoch 86/100\n", "7797/7797 [==============================] - 13s 2ms/step - loss: 0.0165 - acc: 0.9988\n", "Epoch 87/100\n", "7797/7797 [==============================] - 13s 2ms/step - loss: 0.0159 - acc: 0.9992\n", "Epoch 88/100\n", "7797/7797 [==============================] - 13s 2ms/step - loss: 0.0156 - acc: 0.9990\n", "Epoch 89/100\n", "7797/7797 [==============================] - 13s 2ms/step - loss: 0.0152 - acc: 0.9990\n", "Epoch 90/100\n", "7797/7797 [==============================] - 12s 2ms/step - loss: 0.0148 - acc: 0.9988\n", "Epoch 91/100\n", "7797/7797 [==============================] - 12s 2ms/step - loss: 0.0145 - acc: 0.9991\n", "Epoch 92/100\n", "7797/7797 [==============================] - 12s 2ms/step - loss: 0.0144 - acc: 0.9991\n", "Epoch 93/100\n", "7797/7797 [==============================] - 12s 2ms/step - loss: 0.0137 - acc: 0.9991\n", "Epoch 94/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0132 - acc: 0.9995\n", "Epoch 95/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0130 - acc: 0.9991\n", "Epoch 96/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0124 - acc: 0.9996\n", "Epoch 97/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0125 - acc: 0.9988\n", "Epoch 98/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0123 - acc: 0.9994\n", "Epoch 99/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0118 - acc: 0.9991\n", "Epoch 100/100\n", "7797/7797 [==============================] - 11s 1ms/step - loss: 0.0116 - acc: 0.9995\n", "########################################\n", "loss: 0.01164846237482735, acc: 0.9994869821726305\n", "########################################\n" ] } ], "source": [ "@inputlayer\n", "@outputdense('sigmoid')\n", "def double_stranded_model_dnaconv(inputs, inp, oup, params):\n", " \"\"\" keras model for scanning both DNA strands.\n", "\n", " A more elegant way of scanning both strands for motif occurrences\n", " is achieved by the DnaConv2D layer wrapper, which internally\n", " performs the convolution operation with the normal kernel weights\n", " and the reverse complemented weights.\n", " \"\"\"\n", " with inputs.use('dna') as layer:\n", " # the name in inputs.use() should be the same as the dataset name.\n", " layer = DnaConv2D(Conv2D(params[0], (params[1], 1),\n", " activation=params[2]))(layer)\n", " output = GlobalAveragePooling2D(name='motif')(layer)\n", " return inputs, output\n", "\n", "\n", "# create a new model object\n", "model = create_model(template=double_stranded_model_dnaconv,\n", " modelparams=(30, 21, 'relu'),\n", " inputs=DNA,\n", " outputs=ReduceDim(LABELS))\n", "\n", "model.compile(optimizer='adadelta', loss='binary_crossentropy',\n", " metrics=['acc'])\n", "\n", "hist = model.fit(DNA, ReduceDim(LABELS), epochs=epochs)\n", "\n", "print('#' * 40)\n", "print('loss: {}, acc: {}'.format(hist.history['loss'][-1],\n", " hist.history['acc'][-1]))\n", "print('#' * 40)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Perform variant effect prediction" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Model: \"model_1\"\n", "_________________________________________________________________\n", "Layer (type) Output Shape Param # \n", "=================================================================\n", "dna (InputLayer) (None, 198, 1, 64) 0 \n", "_________________________________________________________________\n", "dna_conv2d_1 (DnaConv2D) (None, 178, 1, 30) 40350 \n", "_________________________________________________________________\n", "motif (GlobalAveragePooling2 (None, 30) 0 \n", "_________________________________________________________________\n", "peaks (Dense) (None, 1) 31 \n", "=================================================================\n", "Total params: 40,381\n", "Trainable params: 40,381\n", "Non-trainable params: 0\n", "_________________________________________________________________\n" ] } ], "source": [ "model.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the cell below, we use predict_variant_effect (rather than Janggu.predict_variant_effect).\n", "The function performs variant effect prediction equivalently compared to the method Janggu.predict_variant_effect.\n", "The only difference is that it requires a keras model as the first argument.\n", "\n", "Note also, that as of janggu>=0.10.0, predict_variant_effect and Janggu.predict_variant_effect accept the reference genome\n", "as in fasta format directly (in addition to the Bioseq object, which was supported before only)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# output directory for the variant effect prediction\n", "vcfoutput = os.path.join(os.environ['JANGGU_OUTPUT'], 'vcfoutput')\n", "os.makedirs(vcfoutput, exist_ok=True)\n", "\n", "# perform variant effect prediction using Bioseq object and\n", "# a VCF file\n", "scoresfile, variantsfile = predict_variant_effect(model,\n", " REFGENOME,\n", " VCFFILE,\n", " conditions=['feature'],\n", " output_folder=vcfoutput,\n", " order=order)\n", "\n", "scoresfile = os.path.join(vcfoutput, 'scores.hdf5')\n", "variantsfile = os.path.join(vcfoutput, 'snps.bed.gz')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "scores.hdf5 contains a variety of scores for each variant. The most important ones are refscore and altscore which are used to derive the score difference and the logoddsscore." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "altscore\n", "diffscore\n", "labels\n", "logoddsscore\n", "refscore\n" ] } ], "source": [ "# parse the variant effect predictions (difference between\n", "# reference and alternative variant) into a Cover object\n", "# for the purpose of visualization\n", "f = h5py.File(scoresfile, 'r')\n", "\n", "for name in f:\n", " print(name)\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0. ],\n", " [-0.000977],\n", " [-0.000977],\n", " [-0.000977],\n", " [-0.003906],\n", " [ 0. ]], dtype=float16)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f['diffscore'][:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Use the VariantStreamer for a custom model setup" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far, we have illustrated variant effect predictions by using: 1. Janggu.predict_variant_effect (method from the Janggu model wrapper),\n", " and 2. predict_variant_effect (using a pure keras model).\n", " \n", "In case you want to use a different setup (e.g. an sklearn model) you can use the VariantStreamer to obtain\n", "variants and supply them to the custom model in a loop directly.\n", "This also enables the possibility to adjust the setup if e.g. transformations to the one-hot encoding are desired,\n", "if the output should be reported in a different way (predict_variant_effect outputs a hdf5 file), or if\n", "you want to experiment with different variant effect prediction scores." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "from janggu.data import VariantStreamer" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "variantstreamer = VariantStreamer(REFGENOME,\n", " VCFFILE,\n", " binsize=200,\n", " batch_size=128,\n", " order=order)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We shall just use our existing keras model in the following, but any other model could be used.\n", "Furthermore, transformations on the one-hot encoding could be done as well if necessary." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The VariantStreamer produces mini-batches of pairs of reference and alternative allele variants embedded in the sequence context\n", "around the genomic locus." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0. ]\n", " [-0.000977]\n", " [-0.000977]\n", " [-0.000977]\n", " [-0.003906]\n", " [ 0. ]]\n" ] } ], "source": [ "for names, chroms, positions, ref_alleles, alt_alleles, references, alternatives in variantstreamer.flow():\n", " # names, chroms, positions, ref_alleles, alt_alleles are reported as meta information.\n", " \n", " # references, alternatives represent the one-hot encoded sequences\n", " \n", " # here you could employ any model and if necessary transform the one-hot\n", " # encoded sequences: references and alternatives\n", " ref_score = model.predict_on_batch(references)\n", " alt_score = model.predict_on_batch(alternatives)\n", " \n", " # high score difference indicates potentially high variant effect\n", " \n", " diff_score = alt_score.astype('float16') - ref_score.astype('float16')\n", " # float16 is used to reduce memory requirements, if many variants are tested.\n", " # Also, small variant effects are usually not of interest,\n", " # but this might be adapted if necessary.\n", " break\n", "\n", "print(diff_score)" ] } ], "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }