Datasets
Models
Downloads: 1
718 lines (717 with data), 75.0 kB
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Input attribution and Variant effect prediction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this tutorial we illustrate input feature importance attribution and variant effect prediction."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/wkopp/anaconda3/envs/jdev/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
" from ._conv import register_converters as _register_converters\n",
"Using TensorFlow backend.\n"
]
}
],
"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 Janggu\n",
"from janggu import Scorer\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.data import LineTrack\n",
"from janggu.data import SeqTrack\n",
"from janggu.layers import DnaConv2D\n",
"from janggu import input_attribution\n",
"\n",
"np.random.seed(1234)"
]
},
{
"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",
" store_whole_genome=True,\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",
" store_whole_genome=True,\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": [
"Generated model-id: '6acf30da11d48d4d96ed0669bf3a52f4'\n",
"Epoch 1/100\n",
"244/244 [==============================] - 4s 14ms/step - loss: 0.6253 - acc: 0.6480\n",
"Epoch 2/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.5209 - acc: 0.7661\n",
"Epoch 3/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.4645 - acc: 0.7958\n",
"Epoch 4/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.4248 - acc: 0.8201\n",
"Epoch 5/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.3963 - acc: 0.8330\n",
"Epoch 6/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.3720 - acc: 0.8441\n",
"Epoch 7/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.3488 - acc: 0.8535\n",
"Epoch 8/100\n",
"244/244 [==============================] - 2s 9ms/step - loss: 0.3287 - acc: 0.8642\n",
"Epoch 9/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.3064 - acc: 0.8775\n",
"Epoch 10/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.2868 - acc: 0.8881\n",
"Epoch 11/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.2666 - acc: 0.8989\n",
"Epoch 12/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.2482 - acc: 0.9084\n",
"Epoch 13/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.2305 - acc: 0.9148\n",
"Epoch 14/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.2159 - acc: 0.9233\n",
"Epoch 15/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.2011 - acc: 0.9297\n",
"Epoch 16/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1889 - acc: 0.9346\n",
"Epoch 17/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1777 - acc: 0.9396\n",
"Epoch 18/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1668 - acc: 0.9434\n",
"Epoch 19/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1571 - acc: 0.9465\n",
"Epoch 20/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1495 - acc: 0.9500\n",
"Epoch 21/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1412 - acc: 0.9549\n",
"Epoch 22/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1343 - acc: 0.9576\n",
"Epoch 23/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1277 - acc: 0.9589\n",
"Epoch 24/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1208 - acc: 0.9616\n",
"Epoch 25/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1161 - acc: 0.9632\n",
"Epoch 26/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1106 - acc: 0.9668\n",
"Epoch 27/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1066 - acc: 0.9661\n",
"Epoch 28/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.1021 - acc: 0.9689\n",
"Epoch 29/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0989 - acc: 0.9694\n",
"Epoch 30/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0948 - acc: 0.9721\n",
"Epoch 31/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0911 - acc: 0.9709\n",
"Epoch 32/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0873 - acc: 0.9752\n",
"Epoch 33/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0841 - acc: 0.9750\n",
"Epoch 34/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0811 - acc: 0.9766\n",
"Epoch 35/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0783 - acc: 0.9781\n",
"Epoch 36/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0759 - acc: 0.9789\n",
"Epoch 37/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0732 - acc: 0.9804\n",
"Epoch 38/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0704 - acc: 0.9809\n",
"Epoch 39/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0682 - acc: 0.9834\n",
"Epoch 40/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0655 - acc: 0.9839\n",
"Epoch 41/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0637 - acc: 0.9842\n",
"Epoch 42/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0622 - acc: 0.9841\n",
"Epoch 43/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0600 - acc: 0.9851\n",
"Epoch 44/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0579 - acc: 0.9869\n",
"Epoch 45/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0567 - acc: 0.9858\n",
"Epoch 46/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0547 - acc: 0.9862\n",
"Epoch 47/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0535 - acc: 0.9866\n",
"Epoch 48/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0520 - acc: 0.9874\n",
"Epoch 49/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0502 - acc: 0.9889\n",
"Epoch 50/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0490 - acc: 0.9883\n",
"Epoch 51/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0472 - acc: 0.9889\n",
"Epoch 52/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0462 - acc: 0.9892\n",
"Epoch 53/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0449 - acc: 0.9895\n",
"Epoch 54/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0436 - acc: 0.9914\n",
"Epoch 55/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0423 - acc: 0.9900\n",
"Epoch 56/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0407 - acc: 0.9908\n",
"Epoch 57/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0400 - acc: 0.9913\n",
"Epoch 58/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0388 - acc: 0.9912\n",
"Epoch 59/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0372 - acc: 0.9922\n",
"Epoch 60/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0368 - acc: 0.9926\n",
"Epoch 61/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0355 - acc: 0.9930\n",
"Epoch 62/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0347 - acc: 0.9937\n",
"Epoch 63/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0341 - acc: 0.9932\n",
"Epoch 64/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0329 - acc: 0.9940\n",
"Epoch 65/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0319 - acc: 0.9939\n",
"Epoch 66/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0311 - acc: 0.9950\n",
"Epoch 67/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0310 - acc: 0.9946\n",
"Epoch 68/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0293 - acc: 0.9951\n",
"Epoch 69/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0289 - acc: 0.9946\n",
"Epoch 70/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0281 - acc: 0.9954\n",
"Epoch 71/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0272 - acc: 0.9956\n",
"Epoch 72/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0267 - acc: 0.9959\n",
"Epoch 73/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0259 - acc: 0.9963\n",
"Epoch 74/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0250 - acc: 0.9965\n",
"Epoch 75/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0244 - acc: 0.9965\n",
"Epoch 76/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0245 - acc: 0.9967\n",
"Epoch 77/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0235 - acc: 0.9967\n",
"Epoch 78/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0226 - acc: 0.9974\n",
"Epoch 79/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0221 - acc: 0.9976\n",
"Epoch 80/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0214 - acc: 0.9980\n",
"Epoch 81/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0212 - acc: 0.9978\n",
"Epoch 82/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0203 - acc: 0.9974\n",
"Epoch 83/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0201 - acc: 0.9976\n",
"Epoch 84/100\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"244/244 [==============================] - 2s 10ms/step - loss: 0.0194 - acc: 0.9981\n",
"Epoch 85/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0190 - acc: 0.9982\n",
"Epoch 86/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0184 - acc: 0.9986\n",
"Epoch 87/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0180 - acc: 0.9987\n",
"Epoch 88/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0177 - acc: 0.9983\n",
"Epoch 89/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0171 - acc: 0.9990\n",
"Epoch 90/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0165 - acc: 0.9990\n",
"Epoch 91/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0161 - acc: 0.9991\n",
"Epoch 92/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0161 - acc: 0.9987\n",
"Epoch 93/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0156 - acc: 0.9991\n",
"Epoch 94/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0151 - acc: 0.9995\n",
"Epoch 95/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0145 - acc: 0.9991\n",
"Epoch 96/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0140 - acc: 0.9991\n",
"Epoch 97/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0140 - acc: 0.9996\n",
"Epoch 98/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0136 - acc: 0.9995\n",
"Epoch 99/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0132 - acc: 0.9992\n",
"Epoch 100/100\n",
"244/244 [==============================] - 2s 10ms/step - loss: 0.0129 - acc: 0.9995\n",
"########################################\n",
"loss: 0.012950192026199592, 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 = Janggu.create(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": [
"The toy example illustrates a binary classification example composed of Oct4 and Mafk binding sites. \n",
"The true Oct4 binding sites have been labeled with ones while the Mafk binding sites are labeled zeros."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a sanity check we inspect the predicted values for a few data points."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Oct4 predictions scores should be greater than Mafk scores:\n",
"Prediction score examples for Oct4\n",
"0.: [[[[0.99882954]]]]\n",
"1.: [[[[0.9025411]]]]\n",
"2.: [[[[0.99821955]]]]\n",
"3.: [[[[0.98289168]]]]\n",
"Prediction score examples for Mafk\n",
"1.: [[[[0.59950411]]]]\n",
"2.: [[[[0.00201734]]]]\n",
"3.: [[[[0.00082711]]]]\n",
"4.: [[[[2.10625944e-06]]]]\n"
]
}
],
"source": [
"pred = model.predict(DNA_TEST)\n",
"cov_pred = Cover.create_from_array('BindingProba', pred, LABELS_TEST.gindexer)\n",
"\n",
"print('Oct4 predictions scores should be greater than Mafk scores:')\n",
"print('Prediction score examples for Oct4')\n",
"for i in range(4):\n",
" print('{}.: {}'.format(i, cov_pred[i]))\n",
"print('Prediction score examples for Mafk')\n",
"for i in range(1, 5):\n",
" print('{}.: {}'.format(i, cov_pred[-i]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to perform input feature attribution, we utilize the 'input_attribution' method for a genomic region of interest.\n",
"\n",
"Underneath, the result is illustrated visually. It highlights an Oct4 binding sites occuring at the left peak."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAExCAYAAADx4e+wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFTVJREFUeJzt3X2QpVWdH/Dvb5kZGtGIjLOuOMKMDsQKVBwVLbXUEDC1rJUtBTfZTjaLGDcWFeMU1FIuuqZCsmqUTdW6lC5VvmQVK8xgUIjJVqjd5cWoCYozS8wgqdp2hGVUcBhfVmUHkZz80beHBpuZpvt239vnfj5VU31fn/Prnnvv833OOc+51VoLAABr3y+MugAAAIZDsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AB2oqk9U1XtGXQcwWoIdwISpqt+rqv9TVT+rqstHXQ8wPIIdwOSZSfKOJH8y6kKA4RLsAIagqu6uqndW1der6vtV9cdVNVVVz6yq/1ZVP6iq71XVF6rqFwbPOamqPlNVB6rqm1W1Y972HjO0WlVnVdX+eddfVFV7qupHVXVtkqnH1fMvqmpm0Obnquqkuftaa59srf33JD9ayb8JsPoEO4Dh+Y0kv5zk+UlOS/LuJL+dZH+STUmeleRdSdog3P3XJP87yXOSnJPk4qr65aM1UlUbktyQ5FNJTkzyn5O8cd79Zyf590n+cZJnJ7knya6h/IbAWBPsAIbnQ621e1tr30vy3iT/JMnDmQ1Xp7TWHm6tfaHNfkn3S5Nsaq39u9baT1tr+5J8NMn0Itp5eZL1ST442OZ1SW6fd/9vJPmPrbU9rbWHkrwzySuqasuQfk9gTAl2AMNz77zL9yQ5KcnvZ3ZO259W1b6qumxw/ylJThoM0f6gqn6Q2d68Zy2inZOSfGsQEOe3N//+w9dbaz9OcjCzPYNAx9aNugCAjjx33uWTk3y7tfajzA7H/nZVnZ7klqq6PbMh8JuttVOfYFs/SfKUedd/ad7l7yR5TlXVvHB3cpJvDC5/O7PBMUlSVccn2ZjkW0v7tYC1Qo8dwPC8rao2V9WJme19u7aq/mFVbauqSvLXSR4Z/PtKkr+uqt+pquOq6piqOqOqXjrY1h1JXldVJ1bVLyW5eF47/yvJz5LsqKp1VXV+kpfNu/+aJG+uqu1VdWyS9yX5cmvt7iSpqvVVNZXZfcC6wUkex6zQ3wRYRYIdwPBck+RPk+wb/HtPklOT/HmSH2c2kP1Ra+3W1tojSX41yfYk30zyQJKPJXn6YFufyuyJFXcPtnntXCOttZ8mOT/JhUm+n+TXk3x23v03JfnXST6T2d695+exc/c+muRvMjsH8HcHl39zGH8AYLTqsVM0AFiKqro7yW+11v581LUAk0uPHQBAJwQ7AIBOGIoFAOiEHjsAgE4IdgAAnbBAcT+MqQPA2lfLebIeOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnajW2qhrYAiqam+SQ6OuAwBYlqnW2hlLffK6YVbCSB1qrZ056iIAgKWrqq8u5/mGYgEAOiHYAQB0QrDrx0dGXQAAsGzL2p87eQIAoBN67AAAOiHYAQB0QrBbI6rqhKq6rqr+b1XdVVWvqKoTq+rPquovBz+fMXhsVdWVVTVTVV+rqhePun4AIKmqS6rqzqraW1U7q2qqqrZW1ZcH+/Nrq2rD4LHHDq7PDO7fcrTtC3Zrxx8mubG19oIkL0xyV5LLktzUWjs1yU2D60nyK0lOHfx7a5KrVr9cAGC+qnpOkh1JzhwsQnxMkukkH0jyB4P9+feTvGXwlLck+X5rbVuSPxg87ogEuzWgqv5Wktck+XiStNZ+2lr7QZLXJ/nk4GGfTPKGweXXJ7m6zbotyQlV9exVLhsA+HnrkhxXVeuSPCXJd5KcneS6wf2P35/P7eevS3JOVdWRNi7YrQ3PS3IgyR9X1V9U1ceq6vgkz2qtfSdJBj9/cfD45yS5d97z9w9uAwBGpLX2rST/IclfZTbQ/TDJ7iQ/aK39bPCw+fvsw/vzwf0/TLLxSG0IdmvDuiQvTnJVa+1FSX6SR4ddF7JQmreuDQCM0GAu/OuTbE1yUpLjMzt96vHm9tlPen8u2K0N+5Psb619eXD9uswGvfvnhlgHP7877/HPnff8zUm+vUq1AgALe22Sb7bWDrTWHk7y2SSvzOyUqXWDx8zfZx/enw/uf3qS7x2pAcFuDWit3Zfk3qr624Obzkny9SSfS/KmwW1vSvJfBpc/l+SCwdmxL0/yw7khWwBgZP4qycur6imDuXJz+/Nbkvza4DGP35/P7ed/LcnN7SjfLOGbJ9aIqtqe5GNJNiTZl+TNmQ3mn05ycmZfLP+otfa9wYvlQ0nOTfJgkje31r46ksIBgMOq6t8m+fUkP0vyF0l+K7Nz6XYlOXFw2z9rrT1UVVNJPpXkRZntqZture074vYFOwCAPhiKBQDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOrBt1ASvt3HPPbTfeeOOoyziaWu4GnvnMZ7YtW7YMoRSY9d2H9iVJfvHY5424EoDJsXv37gdaa5uW+vzug90DDzww6hJWxZYtW/LVr3511GXQkStnppMkO7btGnElAJOjqu5ZzvMNxQIAdEKwAwDohGAHANCJ7ufYAQCT4eGHH87+/ftz6NChUZdyVFNTU9m8eXPWr18/1O0KdgBAF/bv35+nPe1p2bJlS6qWveDEimmt5eDBg9m/f3+2bt061G0bigUAunDo0KFs3LhxrENdklRVNm7cuCI9i4IdANCNcQ91c1aqTsEOAKAT5tgBAF2auXB6qNvb9onFLdh+/fXX5/zzz89dd92VF7zgBUOt4Wj02AEADNHOnTvzqle9Krt2rf439wh2AABD8uMf/zhf+tKX8vGPf1ywAwBYy2644Yace+65Oe2003LiiSdmz549q9q+YAcAMCQ7d+7M9PTs3L7p6ens3LlzVdt38gQAwBAcPHgwN998c/bu3ZuqyiOPPJKqyhVXXLFqy7DosQMAGILrrrsuF1xwQe65557cfffduffee7N169Z88YtfXLUa9NgBAF1a7PIkw7Jz585cdtllj7ntjW98Y6655pq8+tWvXpUaBDsAgCG49dZbf+62HTt2rGoNhmIBADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJyx3AgB06cqZ6aFub8e2xa2Ld9999+Xiiy/O7bffnmOPPTZbtmzJBz/4wZx22mlDrWcheuwAAIaktZbzzjsvZ511Vr7xjW/k61//et73vvfl/vvvX5X29dgBAAzJLbfckvXr1+eiiy46fNv27dtXrX09dgAAQ7J379685CUvGVn7gh0AQCcEOwCAITn99NOze/fukbUv2AEADMnZZ5+dhx56KB/96EcP33b77bfn85///Kq07+QJAKBLi12eZJiqKtdff30uvvjivP/978/U1NTh5U5Wg2AHADBEJ510Uj796U+PpG1DsQAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATljsBALo0vXdmqNvbdca2RT3u/vvvzyWXXJLbbrstz3jGM7Jhw4a84x3vyHnnnTfUehaixw4AYEhaa3nDG96Q17zmNdm3b192796dXbt2Zf/+/avSvmAHHDa9d2boR7gwLDMXTmfmwulRlwFHdPPNN2fDhg256KKLDt92yimn5O1vf/uqtC/YAQAMyZ133pkXv/jFI2tfsAMAWCFve9vb8sIXvjAvfelLV6U9wQ4AYEhOP/307Nmz5/D1D3/4w7npppty4MCBVWlfsAMAGJKzzz47hw4dylVXXXX4tgcffHDV2rfcCQDQpcUuTzJMVZUbbrghl1xySa644ops2rQpxx9/fD7wgQ+sSvuCHQDAED372c/Orl27RtK2oVgAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcudAABdunp6Zqjbu2DXkdfFO3jwYM4555wkyX333ZdjjjkmmzZtSpJ85StfyYYNG4Zaz0IEOwCAIdi4cWPuuOOOJMnll1+epz71qbn00ktXtQZDsQCsKVdPzwy9JwZ6IdjBBJi5cDozF06PugwAVphgBwDQCcEOAKATgh0AQCecFQsAdOloy5P0SLADABiyyy+/fCTtGooFAOiEYAcA0AnBDgDoRmtt1CUsykrVKdgBQ2dBZGAUpqamcvDgwbEPd621HDx4MFNTU0PftpMnAIAubN68Ofv378+BAwdGXcpRTU1NZfPmzUPfrmAHAHRh/fr12bp166jLGClDsQAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCeWHOyq6l3zLp9QVf9yOYVU1eVVdelytgEAMMmW02P3rnmXT0iyYLCrqmOW0QYAAIu0qG+eqKobkjw3yVSSP0zyvCTHVdUdSe5MckyS5w+u/1mSP0nyb5J8J8n2JH/nCbb7u0kuSHJvkgNJdg9uvzXJl5P8/cyGxre01r5QVVuSfCrJ8YNN/KvW2v9cYLtvTfLWJDn55JMX8ysCAKx5i/1KsX/eWvteVR2X5PYkfy+zoWp7kgwC1xnzrp+V5GWD27650Aar6iVJppO8aFDHngyC3VxtrbWXVdXrMhsSX5vku0n+QWvtUFWdmmRnkjMfv+3W2keSfCRJzjzzzPH+JmAAgCFZbLDbUVXnDS4/N8mpi3jOV54o1A28Osn1rbUHk6SqPve4+z87+Lk7yZbB5fVJPlRV25M8kuS0RdQBADARjhrsBr1vr03yitbag4Nh0qlFbPsni3jMkXrTHhr8fCSP1nlJkvuTvDCz8wMPLaINADo1vXcmSbLrjG0jrgTGw2JOnnh6ku8PQt0Lkrx8cPvDVbV+cPlHSZ72JNv+H0nOq6rjquppSX51kbV8p7X2/5L8Zmbn9gEAkMUFuxuTrKuqryX5vSS3DW7/SJKvVdV/aq0dTPKlqtpbVb+/mIZba3uSXJvkjiSfSfKFRTztj5K8qapuy+ww7GJ6BQEAJsJRh2Jbaw8l+ZUF7ro1ye/Me9w/XeD+o237vUneu8DtZ827/EAGc+xaa3+Z5O/Oe+g7j9YGAMCk8M0TAACdWOxZsUtWVRuT3LTAXecMhnABABiCFQ92g/C2faXbAQCYdIZiAQA6IdgBAHRixYdigfF39fTsIq9592jrAGB59NgBAHRCsAOOaObC6cxcOD3qMgBYBMEOAKATgh0AQCcEOwCATgh2AACdEOxgjFw5M50rZ5yoAMDSCHYAAJ0Q7AAAOiHYAcAKmd47k+m9MyNr3zqUk0ewAwDohGAHANAJwQ4ARuDq6ZlcPT26YVr6JNgBAHRCsAMA6MS6URcAjMbcENAFu7aNuBJgOeYWNd+xbdfP3fa6kVTEKOmxA4AhMGeOcSDYAQB0QrADuqPnhNWy0ALAXn+MkmAHANAJwQ4AoBOCHYyI4RoAhk2wA1aMLyDn8RzQwMoS7FgTrpyZPrwuEwCwMMEOOqOXjEm0kgd/ehlZSwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADYNmctAPjYd2oCwCGY+6svVdOjbgQAEZGsIM1aC7EXbBr24grGS+P9hi9Z8nb8LcF1jLBDiaM9bjo3VzA3/aJXSOuBFafOXYAAJ0Q7ACewPTemUzvXdkeTicdAMMk2AEAdMIcOwB4kvSyMq702MEY8qXjACyFYAdDcOXMdK6ccQS/VKsxlw1/Z9YOn6lLJ9gBrAC9ro9ayt/CSSWwNObYAbAi5npcTljGgtGT7PG9q7vOGM9FsxdaN3Du/37HNmsJrjbBDkbMNx0szI4BRmMuqL370tlAPq6BkoUZimXVTNrQyrB+33GaF2V4cTIs57U7Tq9XmESCHSMxaSEPWPsc2LAWCHbAmrYaBwk9HoiM41mH41jTuNIzyhMxx45uzX3omR8CMHxzn7FOjRkveuwAYEz02DvM6tJjx4qbm5PyyqkRF7IEC53GvxzOgH2swzuwSx3zAwyDYAewgoZ9cNCL+Qc5h09ImOB8r5du8Y72npr0aTiGYnlCJufyRBY6O9DrZbzM//8wvAeTQ48da46Fa2H19Dx9QG/q6M0/C9pn+nAIdvyc5fS6CF30bqGg03P4YbwMo1d8XD+nVyNoT0KYF+xgYBLe8GudALWyJnX5ikmek+U91R9z7BgpK7nD4ox68V7vVVgb9NixbMPq1l+NI8dJPjIftXHsEZ3/mpt7HZ8wRv1V4zrsOw41QOIzfSGCHU+KN9F4WYs72KO9hg6fvfnu1apoPCzlAGkS3o+T8DsyWr29xgzFsmZZwoFRGvXQKE/Ok/28sHzP2mKqwKP02LFk49SzMo7DfONgUifDj7txPSuRR3vBV/JzbRyH/XmstdyLJ9hNuEnYwQh9AJNtkg5yBTtYBY8OAS39Y2UtH0Eya7G93Hp0GFfLmdc7vyPB18itHMEOABiKcZqisxQ9zKsU7CbIJAxJDqtXy0kZ/ZmE1/9STMJ0DCbTpI5yCHYTaDV2cJOyszCcMC8EXzqef4S5/6NXTq18W5Pyup+zFneca7HmSTb/PTXuB9zjcvAo2HXOhxg9MOfsyMZlhzJuJu3vMmkHFixMsOvUuPeiHMlC3wYwbPM/8OfaeN0SnguwGA6yx89y9pPjvB8Q7Dqx728eWtHt9zChdLU82aDI0qzmEOta4r3KE5l7bXjPPDlzf7fXDfLf0c4IXijEr2YQrNbaijcySlV1Y2vt3FHXsdKq6kCSe0ZdBwCwLKe01jYt9cndBzsAgEnhu2IBADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOvH/AW7gn/jbjS1ZAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 720x360 with 3 Axes>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAExCAYAAADx4e+wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFTVJREFUeJzt3X2QpVWdH/Dvb5kZGtGIjLOuOMKMDsQKVBwVLbXUEDC1rJUtBTfZTjaLGDcWFeMU1FIuuqZCsmqUTdW6lC5VvmQVK8xgUIjJVqjd5cWoCYozS8wgqdp2hGVUcBhfVmUHkZz80beHBpuZpvt239vnfj5VU31fn/Prnnvv833OOc+51VoLAABr3y+MugAAAIZDsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AB2oqk9U1XtGXQcwWoIdwISpqt+rqv9TVT+rqstHXQ8wPIIdwOSZSfKOJH8y6kKA4RLsAIagqu6uqndW1der6vtV9cdVNVVVz6yq/1ZVP6iq71XVF6rqFwbPOamqPlNVB6rqm1W1Y972HjO0WlVnVdX+eddfVFV7qupHVXVtkqnH1fMvqmpm0Obnquqkuftaa59srf33JD9ayb8JsPoEO4Dh+Y0kv5zk+UlOS/LuJL+dZH+STUmeleRdSdog3P3XJP87yXOSnJPk4qr65aM1UlUbktyQ5FNJTkzyn5O8cd79Zyf590n+cZJnJ7knya6h/IbAWBPsAIbnQ621e1tr30vy3iT/JMnDmQ1Xp7TWHm6tfaHNfkn3S5Nsaq39u9baT1tr+5J8NMn0Itp5eZL1ST442OZ1SW6fd/9vJPmPrbU9rbWHkrwzySuqasuQfk9gTAl2AMNz77zL9yQ5KcnvZ3ZO259W1b6qumxw/ylJThoM0f6gqn6Q2d68Zy2inZOSfGsQEOe3N//+w9dbaz9OcjCzPYNAx9aNugCAjjx33uWTk3y7tfajzA7H/nZVnZ7klqq6PbMh8JuttVOfYFs/SfKUedd/ad7l7yR5TlXVvHB3cpJvDC5/O7PBMUlSVccn2ZjkW0v7tYC1Qo8dwPC8rao2V9WJme19u7aq/mFVbauqSvLXSR4Z/PtKkr+uqt+pquOq6piqOqOqXjrY1h1JXldVJ1bVLyW5eF47/yvJz5LsqKp1VXV+kpfNu/+aJG+uqu1VdWyS9yX5cmvt7iSpqvVVNZXZfcC6wUkex6zQ3wRYRYIdwPBck+RPk+wb/HtPklOT/HmSH2c2kP1Ra+3W1tojSX41yfYk30zyQJKPJXn6YFufyuyJFXcPtnntXCOttZ8mOT/JhUm+n+TXk3x23v03JfnXST6T2d695+exc/c+muRvMjsH8HcHl39zGH8AYLTqsVM0AFiKqro7yW+11v581LUAk0uPHQBAJwQ7AIBOGIoFAOiEHjsAgE4IdgAAnbBAcT+MqQPA2lfLebIeOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnajW2qhrYAiqam+SQ6OuAwBYlqnW2hlLffK6YVbCSB1qrZ056iIAgKWrqq8u5/mGYgEAOiHYAQB0QrDrx0dGXQAAsGzL2p87eQIAoBN67AAAOiHYAQB0QrBbI6rqhKq6rqr+b1XdVVWvqKoTq+rPquovBz+fMXhsVdWVVTVTVV+rqhePun4AIKmqS6rqzqraW1U7q2qqqrZW1ZcH+/Nrq2rD4LHHDq7PDO7fcrTtC3Zrxx8mubG19oIkL0xyV5LLktzUWjs1yU2D60nyK0lOHfx7a5KrVr9cAGC+qnpOkh1JzhwsQnxMkukkH0jyB4P9+feTvGXwlLck+X5rbVuSPxg87ogEuzWgqv5Wktck+XiStNZ+2lr7QZLXJ/nk4GGfTPKGweXXJ7m6zbotyQlV9exVLhsA+HnrkhxXVeuSPCXJd5KcneS6wf2P35/P7eevS3JOVdWRNi7YrQ3PS3IgyR9X1V9U1ceq6vgkz2qtfSdJBj9/cfD45yS5d97z9w9uAwBGpLX2rST/IclfZTbQ/TDJ7iQ/aK39bPCw+fvsw/vzwf0/TLLxSG0IdmvDuiQvTnJVa+1FSX6SR4ddF7JQmreuDQCM0GAu/OuTbE1yUpLjMzt96vHm9tlPen8u2K0N+5Psb619eXD9uswGvfvnhlgHP7877/HPnff8zUm+vUq1AgALe22Sb7bWDrTWHk7y2SSvzOyUqXWDx8zfZx/enw/uf3qS7x2pAcFuDWit3Zfk3qr624Obzkny9SSfS/KmwW1vSvJfBpc/l+SCwdmxL0/yw7khWwBgZP4qycur6imDuXJz+/Nbkvza4DGP35/P7ed/LcnN7SjfLOGbJ9aIqtqe5GNJNiTZl+TNmQ3mn05ycmZfLP+otfa9wYvlQ0nOTfJgkje31r46ksIBgMOq6t8m+fUkP0vyF0l+K7Nz6XYlOXFw2z9rrT1UVVNJPpXkRZntqZture074vYFOwCAPhiKBQDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOrBt1ASvt3HPPbTfeeOOoyziaWu4GnvnMZ7YtW7YMoRSY9d2H9iVJfvHY5424EoDJsXv37gdaa5uW+vzug90DDzww6hJWxZYtW/LVr3511GXQkStnppMkO7btGnElAJOjqu5ZzvMNxQIAdEKwAwDohGAHANCJ7ufYAQCT4eGHH87+/ftz6NChUZdyVFNTU9m8eXPWr18/1O0KdgBAF/bv35+nPe1p2bJlS6qWveDEimmt5eDBg9m/f3+2bt061G0bigUAunDo0KFs3LhxrENdklRVNm7cuCI9i4IdANCNcQ91c1aqTsEOAKAT5tgBAF2auXB6qNvb9onFLdh+/fXX5/zzz89dd92VF7zgBUOt4Wj02AEADNHOnTvzqle9Krt2rf439wh2AABD8uMf/zhf+tKX8vGPf1ywAwBYy2644Yace+65Oe2003LiiSdmz549q9q+YAcAMCQ7d+7M9PTs3L7p6ens3LlzVdt38gQAwBAcPHgwN998c/bu3ZuqyiOPPJKqyhVXXLFqy7DosQMAGILrrrsuF1xwQe65557cfffduffee7N169Z88YtfXLUa9NgBAF1a7PIkw7Jz585cdtllj7ntjW98Y6655pq8+tWvXpUaBDsAgCG49dZbf+62HTt2rGoNhmIBADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJyx3AgB06cqZ6aFub8e2xa2Ld9999+Xiiy/O7bffnmOPPTZbtmzJBz/4wZx22mlDrWcheuwAAIaktZbzzjsvZ511Vr7xjW/k61//et73vvfl/vvvX5X29dgBAAzJLbfckvXr1+eiiy46fNv27dtXrX09dgAAQ7J379685CUvGVn7gh0AQCcEOwCAITn99NOze/fukbUv2AEADMnZZ5+dhx56KB/96EcP33b77bfn85///Kq07+QJAKBLi12eZJiqKtdff30uvvjivP/978/U1NTh5U5Wg2AHADBEJ510Uj796U+PpG1DsQAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATljsBALo0vXdmqNvbdca2RT3u/vvvzyWXXJLbbrstz3jGM7Jhw4a84x3vyHnnnTfUehaixw4AYEhaa3nDG96Q17zmNdm3b192796dXbt2Zf/+/avSvmAHHDa9d2boR7gwLDMXTmfmwulRlwFHdPPNN2fDhg256KKLDt92yimn5O1vf/uqtC/YAQAMyZ133pkXv/jFI2tfsAMAWCFve9vb8sIXvjAvfelLV6U9wQ4AYEhOP/307Nmz5/D1D3/4w7npppty4MCBVWlfsAMAGJKzzz47hw4dylVXXXX4tgcffHDV2rfcCQDQpcUuTzJMVZUbbrghl1xySa644ops2rQpxx9/fD7wgQ+sSvuCHQDAED372c/Orl27RtK2oVgAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcudAABdunp6Zqjbu2DXkdfFO3jwYM4555wkyX333ZdjjjkmmzZtSpJ85StfyYYNG4Zaz0IEOwCAIdi4cWPuuOOOJMnll1+epz71qbn00ktXtQZDsQCsKVdPzwy9JwZ6IdjBBJi5cDozF06PugwAVphgBwDQCcEOAKATgh0AQCecFQsAdOloy5P0SLADABiyyy+/fCTtGooFAOiEYAcA0AnBDgDoRmtt1CUsykrVKdgBQ2dBZGAUpqamcvDgwbEPd621HDx4MFNTU0PftpMnAIAubN68Ofv378+BAwdGXcpRTU1NZfPmzUPfrmAHAHRh/fr12bp166jLGClDsQAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCeWHOyq6l3zLp9QVf9yOYVU1eVVdelytgEAMMmW02P3rnmXT0iyYLCrqmOW0QYAAIu0qG+eqKobkjw3yVSSP0zyvCTHVdUdSe5MckyS5w+u/1mSP0nyb5J8J8n2JH/nCbb7u0kuSHJvkgNJdg9uvzXJl5P8/cyGxre01r5QVVuSfCrJ8YNN/KvW2v9cYLtvTfLWJDn55JMX8ysCAKx5i/1KsX/eWvteVR2X5PYkfy+zoWp7kgwC1xnzrp+V5GWD27650Aar6iVJppO8aFDHngyC3VxtrbWXVdXrMhsSX5vku0n+QWvtUFWdmmRnkjMfv+3W2keSfCRJzjzzzPH+JmAAgCFZbLDbUVXnDS4/N8mpi3jOV54o1A28Osn1rbUHk6SqPve4+z87+Lk7yZbB5fVJPlRV25M8kuS0RdQBADARjhrsBr1vr03yitbag4Nh0qlFbPsni3jMkXrTHhr8fCSP1nlJkvuTvDCz8wMPLaINADo1vXcmSbLrjG0jrgTGw2JOnnh6ku8PQt0Lkrx8cPvDVbV+cPlHSZ72JNv+H0nOq6rjquppSX51kbV8p7X2/5L8Zmbn9gEAkMUFuxuTrKuqryX5vSS3DW7/SJKvVdV/aq0dTPKlqtpbVb+/mIZba3uSXJvkjiSfSfKFRTztj5K8qapuy+ww7GJ6BQEAJsJRh2Jbaw8l+ZUF7ro1ye/Me9w/XeD+o237vUneu8DtZ827/EAGc+xaa3+Z5O/Oe+g7j9YGAMCk8M0TAACdWOxZsUtWVRuT3LTAXecMhnABABiCFQ92g/C2faXbAQCYdIZiAQA6IdgBAHRixYdigfF39fTsIq9592jrAGB59NgBAHRCsAOOaObC6cxcOD3qMgBYBMEOAKATgh0AQCcEOwCATgh2AACdEOxgjFw5M50rZ5yoAMDSCHYAAJ0Q7AAAOiHYAcAKmd47k+m9MyNr3zqUk0ewAwDohGAHANAJwQ4ARuDq6ZlcPT26YVr6JNgBAHRCsAMA6MS6URcAjMbcENAFu7aNuBJgOeYWNd+xbdfP3fa6kVTEKOmxA4AhMGeOcSDYAQB0QrADuqPnhNWy0ALAXn+MkmAHANAJwQ4AoBOCHYyI4RoAhk2wA1aMLyDn8RzQwMoS7FgTrpyZPrwuEwCwMMEOOqOXjEm0kgd/ehlZSwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADYNmctAPjYd2oCwCGY+6svVdOjbgQAEZGsIM1aC7EXbBr24grGS+P9hi9Z8nb8LcF1jLBDiaM9bjo3VzA3/aJXSOuBFafOXYAAJ0Q7ACewPTemUzvXdkeTicdAMMk2AEAdMIcOwB4kvSyMq702MEY8qXjACyFYAdDcOXMdK6ccQS/VKsxlw1/Z9YOn6lLJ9gBrAC9ro9ayt/CSSWwNObYAbAi5npcTljGgtGT7PG9q7vOGM9FsxdaN3Du/37HNmsJrjbBDkbMNx0szI4BRmMuqL370tlAPq6BkoUZimXVTNrQyrB+33GaF2V4cTIs57U7Tq9XmESCHSMxaSEPWPsc2LAWCHbAmrYaBwk9HoiM41mH41jTuNIzyhMxx45uzX3omR8CMHxzn7FOjRkveuwAYEz02DvM6tJjx4qbm5PyyqkRF7IEC53GvxzOgH2swzuwSx3zAwyDYAewgoZ9cNCL+Qc5h09ImOB8r5du8Y72npr0aTiGYnlCJufyRBY6O9DrZbzM//8wvAeTQ48da46Fa2H19Dx9QG/q6M0/C9pn+nAIdvyc5fS6CF30bqGg03P4YbwMo1d8XD+nVyNoT0KYF+xgYBLe8GudALWyJnX5ikmek+U91R9z7BgpK7nD4ox68V7vVVgb9NixbMPq1l+NI8dJPjIftXHsEZ3/mpt7HZ8wRv1V4zrsOw41QOIzfSGCHU+KN9F4WYs72KO9hg6fvfnu1apoPCzlAGkS3o+T8DsyWr29xgzFsmZZwoFRGvXQKE/Ok/28sHzP2mKqwKP02LFk49SzMo7DfONgUifDj7txPSuRR3vBV/JzbRyH/XmstdyLJ9hNuEnYwQh9AJNtkg5yBTtYBY8OAS39Y2UtH0Eya7G93Hp0GFfLmdc7vyPB18itHMEOABiKcZqisxQ9zKsU7CbIJAxJDqtXy0kZ/ZmE1/9STMJ0DCbTpI5yCHYTaDV2cJOyszCcMC8EXzqef4S5/6NXTq18W5Pyup+zFneca7HmSTb/PTXuB9zjcvAo2HXOhxg9MOfsyMZlhzJuJu3vMmkHFixMsOvUuPeiHMlC3wYwbPM/8OfaeN0SnguwGA6yx89y9pPjvB8Q7Dqx728eWtHt9zChdLU82aDI0qzmEOta4r3KE5l7bXjPPDlzf7fXDfLf0c4IXijEr2YQrNbaijcySlV1Y2vt3FHXsdKq6kCSe0ZdBwCwLKe01jYt9cndBzsAgEnhu2IBADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOvH/AW7gn/jbjS1ZAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 720x360 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Extract the 4th interval to perform input feature importance attribution\n",
"# which represents an Oct4 bound region\n",
"gi = DNA.gindexer[3]\n",
"chrom = gi.chrom\n",
"start = gi.start\n",
"end = gi.end\n",
"attr_oct = input_attribution(model, DNA, chrom=chrom, start=start, end=end)\n",
"\n",
"# visualize the important sequence features\n",
"plotGenomeTrack(SeqTrack(attr_oct[0]), chrom, start, end)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By comparison, the input attribution for a Mafk binding sites highlights a Mafk motif in the center."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAoMAAAExCAYAAAAHhKhtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFjxJREFUeJzt3X+U5WV9H/D3R9h1rJIIy9aIK+wqcDzCiaui0Ry1FGwknsaK9jSTpkHbtJZE5UD0GH/1yPFXFdsTw4nSSrQWT9nVaKAmaUwUxKgNgkuoWfAkDghlFRAWbTRkEcnTP+bOMuDs7MzOnbl37vN6nTNnv/e5d77fz52de7/veZ7n+9xqrQUAgD49YtQFAAAwOsIgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGASYAFX10ap656jrANYfYRCgM1X1jqr6y6r6UVWdP+p6gNESBgH6M5PkDUn+aNSFAKMnDAIMQVXdUlVvqqobq+q7VfXfqmqqqo6uqj+squ9V1T1V9cWqesTge46pqk9V1V1V9c2qOmfe/h4y7FtVp1bVnnm3n15V11XV96vq40mmHlbPv6uqmcExP11Vx8zd11r77621P07y/dX8mQDrgzAIMDy/nORFSZ6c5MQkb03yuiR7kmxO8rgkb07SBoHwD5L8nyRPSHJ6knOr6kUHO0hVbUxyeZKPJTkqye8lefm8+09L8h+T/Iskj09ya5KdQ3mGwMQRBgGG53daa7e11u5J8q4kv5Tk/swGsuNaa/e31r7YZj8U/llJNrfW3t5a+2Fr7eYkFyeZXsJxnpNkQ5L3D/b5ySTXzrv/l5N8pLV2XWvtviRvSvLcqto6pOcJTBBhEGB4bpu3fWuSY5K8L7Nz9P60qm6uqjcO7j8uyTGD4ePvVdX3Mttr+LglHOeYJN8ahMr5x5t///7brbUfJNmb2R5IgIc4fNQFAEyQJ87bPjbJt1tr38/sUPHrquqkJJ+vqmszGxy/2Vo74QD7+tsk/2De7Z+at317kidUVc0LhMcmuWmw/e3Mhs0kSVU9OsmmJN86tKcFTDI9gwDD8+qq2lJVR2W2l+/jVfVPq+r4qqokf5PkgcHXNUn+pqp+s6oeVVWHVdXJVfWswb6uT/Liqjqqqn4qybnzjvPnSX6U5JyqOryqXpbk2fPuvzTJv66q7VX1yCTvTvKV1totSVJVG6pqKrPngMMHF7octko/E2DMCYMAw3Npkj9NcvPg651JTkjyuSQ/yGyI+2Br7arW2gNJfiHJ9iTfTHJ3kt9N8pODfX0ssxeX3DLY58fnDtJa+2GSlyV5ZZLvJvnFJL8/7/4rkvyHJJ/KbC/ik/PQuYgXJ/m7zM5pfMtg+1eG8QMA1p966JQTAA5FVd2S5N+21j436loAlkPPIABAx4RBAICOGSYGAOiYnkEAgI4JgwAAHbPo9OQw3g8A61+t9QH1DAIAdEwYBADomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB0TBgEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JgwCADQMWEQAKBjwiAAQMeEQQCAjgmDAAAdEwYBADomDAIAdEwYBADomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB0TBgEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6Fi11kZdA0NQVbuT7Bt1HQDAiky11k5eywMevpYHY1Xta62dMuoiAIBDV1VfXetjGiYGAOiYMAgA0DFhcHJ8aNQFAAArtubncxeQAAB0TM8gAEDHhEEAgI4Jg2usqj5SVd8ZrAs413Z+VX2rqq4ffL140L61qv5uXvt/mfc9v1hVX6uqG6rqgnntv1FVNw7uu6Kqjpt33yuq6huDr1fMa39mVf1lVc1U1YVVVav/kwCA9W21z+nz7v/nVdWq6pR5bW8anLf/qqpeNK/9jEHbTFW9cSnPQxhcex9NcsYC7b/VWts++Ppf89pvmtd+dpJU1aYk70tyemvtpCSPq6rTB4//iySntNZ+Osknk1ww+J6jkrwtyc8keXaSt1XVkYPvuSjJq5KcMPhaqD4A4KE+mtU9p6eqjkhyTpKvzGt7apLpJCcNjv/Bqjqsqg5L8oEkP5/kqUl+afDYRQmDa6y19mdJ7lnhbp6U5K9ba3cNbn8uycsH+/98a+3eQfvVSbYMtl+U5LOttXtaa99N8tkkZ1TV45P8RGvtz9vs1USXJHnpCusDgIm32uf0gXdktmNn/qeM/bMkO1tr97XWvplkJrMdPc9OMtNau7m19sMkOwePXZQwOD5eM+gi/si8Hrsk2VZVf1FVX6iq5w/aZpI8ZdDlfHhmw9sTF9jnryb548H2E5LcNu++PYO2Jwy2H94OAByaoZzTq+rpSZ7YWvvDh+1/sXP6Qu2LEgbHw0VJnpxke5Lbk/znQfvtSY5trT09yW8kubSqfmLQs/drST6e5ItJbknyo/k7rKp/leSUzHY9J8lC8wDbIu0AwPIN5ZxeVY9I8ltJXrfAMYZ6ThcGx0Br7c7W2gOttb9PcnFmu3kz6P7dO9jeleSmJCcObv9Ba+1nWmvPTfJXSb4xt7+qemGStyR5SWvtvkHznjy093BLkm8P2rcs0A4ALNMQz+lHJDk5yVVVdUuS5yT59OAiksXO6Qu1L0oYHAODeXtzzkyye9C+eTAZNFX1pMxe3HHz4PY/HPx7ZJJfT/K7g9tPT/JfMxsEvzNvv3+S5Oeq6sjB9/xckj9prd2e5PtV9ZzBVcRnJfmfq/ZkAWCCDeuc3lr7f621o1trW1trWzN7HcBLWmtfTfLpJNNV9ciq2jbY1zVJrk1yQlVtq6qNmb3I5NMHq/nwITxvlqGqdiQ5NcnRVbUns1f4nlpV2zPblXtLkn8/ePgLkry9qn6U5IEkZ7fW5iaq/nZVPW2w/fbW2l8Ptt+X5DFJfm+wQsz/ba29pLV2T1W9I7O/KHPfM7evX8vsFVGPyuwcw7l5hgDAAazBOX1BrbUbquoTSW7M7DSxV7fWHhjU9JrMdgAdluQjrbUbDvo8fBwdAEC/DBMDAHRMGAQA6JgwCADQMWEQAKBjwiAAQMeEQQCAjgmDAAAdEwYBADomDAIAdEwYBADomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB0TBgEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JgwCADQMWEQAKBjwiAAQMeEQQCAjgmDAAAdEwYBADomDAIAdEwYBADomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB0TBgEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JgwCADQMWEQAKBjwiAAQMcOH3UBq+2MM85on/nMZ0ZdxsHUSndw9NFHt61btw6hFABgVHbt2nV3a23zWh5z4sPg3XffPeoS1sTWrVvz1a9+ddRlAAArUFW3rvUxDRMDAHRMGAQA6JgwCADQsYmfMwgA9OH+++/Pnj17sm/fvlGXclBTU1PZsmVLNmzYMOpShEEAYDLs2bMnRxxxRLZu3ZqqFS/UsWpaa9m7d2/27NmTbdu2jbocw8QAwGTYt29fNm3aNNZBMEmqKps2bRqbHkxhEGDMXDI9k0umZ0ZdBqxL4x4E54xTncIgAEDHzBkEACbSzCunh7q/4z+6c0mPu+yyy/Kyl70sX//61/OUpzxlqDWsBj2DAABDtGPHjjzvec/Lzp1LC4+jJgwCAAzJD37wg3z5y1/Ohz/8YWEQAKA3l19+ec4444yceOKJOeqoo3LdddeNuqSDEgYBAIZkx44dmZ6enas4PT2dHTt2jLiig3MBCQDAEOzduzdXXnlldu/enarKAw88kKrKBRdcMFZLyTycnkEAgCH45Cc/mbPOOiu33nprbrnlltx2223Ztm1bvvSlL426tEXpGQQAJtJSl4IZlh07duSNb3zjQ9pe/vKX59JLL83zn//8Na1lOYRBAIAhuOqqq36s7Zxzzln7QpbJMDEAQMeEQQCAjgmDAGNqevdMpnfPjLoMYMIJgwAAHRMGAQA6JgwCAHTM0jIAwES6cGZ6qPs75/ilrVt4xx135Nxzz821116bRz7ykdm6dWve//7358QTTxxqPcOiZxAAYEhaaznzzDNz6qmn5qabbsqNN96Yd7/73bnzzjtHXdoB6RkEABiSz3/+89mwYUPOPvvs/W3bt28fYUUHp2cQAGBIdu/enWc+85mjLmNZhEGAEbKWIDBqwiAAwJCcdNJJ2bVr16jLWBZhEABgSE477bTcd999ufjii/e3XXvttfnCF74wwqoW5wISAGAiLXUpmGGqqlx22WU599xz8573vCdTU1P7l5YZV8IgAMAQHXPMMfnEJz4x6jKWzDAxAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JilZQCAiTTsj3rcefLxS3rcnXfemfPOOy9XX311jjzyyGzcuDFveMMbcuaZZw61nmHRMwgAMCSttbz0pS/NC17wgtx8883ZtWtXdu7cmT179oy6tAMSBgHW2CXTM7lkerg9FsB4uPLKK7Nx48acffbZ+9uOO+64vPa1rx1hVYsTBgEAhuSGG27IM57xjFGXsSzCIADAKnn1q1+dpz3taXnWs5416lIOSBgEABiSk046Kdddd93+2x/4wAdyxRVX5K677hphVYsTBgEAhuS0007Lvn37ctFFF+1vu/fee0dY0cFZWgYAmEhLXQpmmKoql19+ec4777xccMEF2bx5cx796Efnve9975rXslTCIADAED3+8Y/Pzp07R13GkhkmBgDomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB2ztAwAMJEumZ4Z6v7O2rn4uoV79+7N6aefniS54447cthhh2Xz5s1JkmuuuSYbN24caj3DIgwCAAzBpk2bcv311ydJzj///DzmMY/J61//+hFXdXCGiQEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY64mBgAm0sGWgmGWMAgwBi6cmU6SnHP8zhFXAgzD+eefP+oSlswwMQBAx4RBAICOCYMAwMRorY26hCUZpzqFQQBgIkxNTWXv3r1jFbQW0lrL3r17MzU1NepSkriABACYEFu2bMmePXty1113jbqUg5qamsqWLVtGXUYSYRAAmBAbNmzItm3bRl3GumOYGACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6Nghh8GqevO87cdW1a+vpJCqOr+qXr+SfQAAsDwr6Rl887ztxyZZMAxW1WErOAYAAKtoSZ9AUlWXJ3likqkkv53kSUkeVVXXJ7khyWFJnjy4/dkkf5TkbUluT7I9yVMPsN+3JDkryW1J7kqya9B+VZKvJPnHmQ2av9pa+2JVbU3ysSSPHuziNa21/73Afl+V5FVJcuyxxy7lKQIAdGmpH0f3b1pr91TVo5Jcm+QfZTaIbU+SQUg7ed7tU5M8e9D2zYV2WFXPTDKd5OmDOq7LIAzO1dZae3ZVvTizwfKFSb6T5J+01vZV1QlJdiQ55eH7bq19KMmHkuSUU04Z70+rBgAYoaWGwXOq6szB9hOTnLCE77nmQEFw4PlJLmut3ZskVfXph93/+4N/dyXZOtjekOR3qmp7kgeSnLiEOgAAOICDhsFBL98Lkzy3tXbvYAh3agn7/tslPGaxXrv7Bv8+kAfrPC/JnUmeltn5jvuWcAwAAA5gKReQ/GSS7w6C4FOSPGfQfn9VbRhsfz/JEcs89p8lObOqHlVVRyT5hSXWcntr7e+T/Epm5yoCAHCIlhIGP5Pk8Kr6WpJ3JLl60P6hJF+rqv/RWtub5MtVtbuq3reUA7fWrkvy8STXJ/lUki8u4ds+mOQVVXV1ZoeIl9L7CADAARx0mLi1dl+Sn1/grquS/Oa8x/3LBe4/2L7fleRdC7SfOm/77gzmDLbWvpHkp+c99E0HOwYAAAfmE0gAADq21KuJD1lVbUpyxQJ3nT4YXgaAdWl690ySZOfJx+eS6dnts3YeP8qSYNlWPQwOAt/21T4OAADLZ5gYAKBjwiAADFw4M50LZ6ZHXQasKWEQAKBjwiAAQMeEQQCAjgmDABNoevfM/mVPABYjDAIAdEwYBADomDAIANAxYRAAoGPCIMA6N/PK6cy80kLJwKERBgEAOiYMAgB0TBgEAOiYMAgA0DFhEKATLjQBFiIMAgB07PBRFwDQiwd75d455P0BHDo9gwAAHRMGAQA6JgwCAHRMGAQA6JgwCLBClmxhvgtnpnPhjN8H1g9hEACgY8IgAEDHhEEAuja9eybTu2dGXQaMjDAIANAxYRDowiXTM7lkev33/rg4ARg2YRAAJsCk/MHD2hMGAQA6JgwCAHRMGAQA6JgwCACrxKfTsB4IgwCwTEIek+TwURcAwPDsv5r0raOtA1g/9AwCwAhYM5JxIQwCAD9GWO2HMAhAEvPgoFfCIABAx4RBuucjnADomTAIsAyGUsef/6O++QN/+SwtQ5fm3ijO2nn8iCthtT0YCt450jpYXV7TcOj0DAIwkfQQwdIIgwBDYngSWI+EQQAm3nJ7CUcV7K3txygIg8DITe+eyfTuyRzOM1TJOJuU8Ol1tjLCILBkhkGBQ+X9Y3wJgwAd0pPSr0nuiefQWFoGWNQ4LNkxd+LaefL4LBsyjjWxPPvD8FtHWweMmp5BAICO6RmECTE3Cfyc43eOuJL1yc9v9fjZ9s08wfGnZxCAA1rLSf+TcmUrrDd6Bpl4K+mVmDsJHv9RPRrA+tbzHEnv5YvTMwhLtNAVeGvZa2JZhodyNSyMr2Fdsay3eG3oGQTW1MF6aiflKt0Hg/s7R1oH42fud/xnp0ZcyAiNQ0/dYjWMwyoKa0kYhDG0Fm9E4/Bmtz8wvf7HA9M41Md4ckHKeHj4SIUh2PVLGGRd8Oa/uobxV/rBwtv+IV0dZWNvHHptYBTmem1fPOI61po5gzBhzC08OPMNAR6kZxDGiIDCapjrWV/N3o7lDuubBrC+6C2ebMIgXZk7KT52DcYqJ/lkN8nPbSGmKbBSa31hVG+vUVZGGIQxN3cSGdepdq6ahdUl2M3SO7l6hEGYYE4iMDxr/Xpa6A8tr2lWgzAIACM2Tgsr7w+hHX5SSa+EwQ6N61+W+4dD/9PsO1AvQwGGPpZuVPOujIA/1FJ/Z4cRKuYfS0iB1SEMcsjW4sS83Df/UQXd+RcYLDa087NTww26gmQ/VisIreR1PKzX2zj1ijEe1mPwX8/vx8IgK3awk8mahsZ11IWzkgtDrCM4HJPy0Xcs3aEEWEs+PWg9B55D1cOyScLghFrJSW7uxf7WwUeELWcfi71RDOtNZKn1rdYQn2VGxlePJ6px4/9gbfg5M0zC4AQ5lDeHpQablQTEhfQ6P3AcrOVai6MyDkOfk2Zcw4deu8lxKL9jy/34uIVe3yt5zS/1YzjH/f1EGCTJ8N/ox/XEsZj18qJdTQtO1j+E0LiWH2C/Hn/XerYep3SMu8UC0VL/MBrWiMfc8VaT1/zwVWtt1DWsqqr6TGvtjFHXsdqq6q4kt466DgBgRY5rrW1eywNOfBgEAODAHjHqAgAAGB1hEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JgwCADQsf8P5hZ9yHaPZ3AAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x360 with 3 Axes>"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAoMAAAExCAYAAAAHhKhtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFjxJREFUeJzt3X+U5WV9H/D3R9h1rJIIy9aIK+wqcDzCiaui0Ry1FGwknsaK9jSTpkHbtJZE5UD0GH/1yPFXFdsTw4nSSrQWT9nVaKAmaUwUxKgNgkuoWfAkDghlFRAWbTRkEcnTP+bOMuDs7MzOnbl37vN6nTNnv/e5d77fz52de7/veZ7n+9xqrQUAgD49YtQFAAAwOsIgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGASYAFX10ap656jrANYfYRCgM1X1jqr6y6r6UVWdP+p6gNESBgH6M5PkDUn+aNSFAKMnDAIMQVXdUlVvqqobq+q7VfXfqmqqqo6uqj+squ9V1T1V9cWqesTge46pqk9V1V1V9c2qOmfe/h4y7FtVp1bVnnm3n15V11XV96vq40mmHlbPv6uqmcExP11Vx8zd11r77621P07y/dX8mQDrgzAIMDy/nORFSZ6c5MQkb03yuiR7kmxO8rgkb07SBoHwD5L8nyRPSHJ6knOr6kUHO0hVbUxyeZKPJTkqye8lefm8+09L8h+T/Iskj09ya5KdQ3mGwMQRBgGG53daa7e11u5J8q4kv5Tk/swGsuNaa/e31r7YZj8U/llJNrfW3t5a+2Fr7eYkFyeZXsJxnpNkQ5L3D/b5ySTXzrv/l5N8pLV2XWvtviRvSvLcqto6pOcJTBBhEGB4bpu3fWuSY5K8L7Nz9P60qm6uqjcO7j8uyTGD4ePvVdX3Mttr+LglHOeYJN8ahMr5x5t///7brbUfJNmb2R5IgIc4fNQFAEyQJ87bPjbJt1tr38/sUPHrquqkJJ+vqmszGxy/2Vo74QD7+tsk/2De7Z+at317kidUVc0LhMcmuWmw/e3Mhs0kSVU9OsmmJN86tKcFTDI9gwDD8+qq2lJVR2W2l+/jVfVPq+r4qqokf5PkgcHXNUn+pqp+s6oeVVWHVdXJVfWswb6uT/Liqjqqqn4qybnzjvPnSX6U5JyqOryqXpbk2fPuvzTJv66q7VX1yCTvTvKV1totSVJVG6pqKrPngMMHF7octko/E2DMCYMAw3Npkj9NcvPg651JTkjyuSQ/yGyI+2Br7arW2gNJfiHJ9iTfTHJ3kt9N8pODfX0ssxeX3DLY58fnDtJa+2GSlyV5ZZLvJvnFJL8/7/4rkvyHJJ/KbC/ik/PQuYgXJ/m7zM5pfMtg+1eG8QMA1p966JQTAA5FVd2S5N+21j436loAlkPPIABAx4RBAICOGSYGAOiYnkEAgI4JgwAAHbPo9OQw3g8A61+t9QH1DAIAdEwYBADomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB0TBgEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JgwCADQMWEQAKBjwiAAQMeEQQCAjgmDAAAdEwYBADomDAIAdEwYBADomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB0TBgEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6Fi11kZdA0NQVbuT7Bt1HQDAiky11k5eywMevpYHY1Xta62dMuoiAIBDV1VfXetjGiYGAOiYMAgA0DFhcHJ8aNQFAAArtubncxeQAAB0TM8gAEDHhEEAgI4Jg2usqj5SVd8ZrAs413Z+VX2rqq4ffL140L61qv5uXvt/mfc9v1hVX6uqG6rqgnntv1FVNw7uu6Kqjpt33yuq6huDr1fMa39mVf1lVc1U1YVVVav/kwCA9W21z+nz7v/nVdWq6pR5bW8anLf/qqpeNK/9jEHbTFW9cSnPQxhcex9NcsYC7b/VWts++Ppf89pvmtd+dpJU1aYk70tyemvtpCSPq6rTB4//iySntNZ+Osknk1ww+J6jkrwtyc8keXaSt1XVkYPvuSjJq5KcMPhaqD4A4KE+mtU9p6eqjkhyTpKvzGt7apLpJCcNjv/Bqjqsqg5L8oEkP5/kqUl+afDYRQmDa6y19mdJ7lnhbp6U5K9ba3cNbn8uycsH+/98a+3eQfvVSbYMtl+U5LOttXtaa99N8tkkZ1TV45P8RGvtz9vs1USXJHnpCusDgIm32uf0gXdktmNn/qeM/bMkO1tr97XWvplkJrMdPc9OMtNau7m19sMkOwePXZQwOD5eM+gi/si8Hrsk2VZVf1FVX6iq5w/aZpI8ZdDlfHhmw9sTF9jnryb548H2E5LcNu++PYO2Jwy2H94OAByaoZzTq+rpSZ7YWvvDh+1/sXP6Qu2LEgbHw0VJnpxke5Lbk/znQfvtSY5trT09yW8kubSqfmLQs/drST6e5ItJbknyo/k7rKp/leSUzHY9J8lC8wDbIu0AwPIN5ZxeVY9I8ltJXrfAMYZ6ThcGx0Br7c7W2gOttb9PcnFmu3kz6P7dO9jeleSmJCcObv9Ba+1nWmvPTfJXSb4xt7+qemGStyR5SWvtvkHznjy093BLkm8P2rcs0A4ALNMQz+lHJDk5yVVVdUuS5yT59OAiksXO6Qu1L0oYHAODeXtzzkyye9C+eTAZNFX1pMxe3HHz4PY/HPx7ZJJfT/K7g9tPT/JfMxsEvzNvv3+S5Oeq6sjB9/xckj9prd2e5PtV9ZzBVcRnJfmfq/ZkAWCCDeuc3lr7f621o1trW1trWzN7HcBLWmtfTfLpJNNV9ciq2jbY1zVJrk1yQlVtq6qNmb3I5NMHq/nwITxvlqGqdiQ5NcnRVbUns1f4nlpV2zPblXtLkn8/ePgLkry9qn6U5IEkZ7fW5iaq/nZVPW2w/fbW2l8Ptt+X5DFJfm+wQsz/ba29pLV2T1W9I7O/KHPfM7evX8vsFVGPyuwcw7l5hgDAAazBOX1BrbUbquoTSW7M7DSxV7fWHhjU9JrMdgAdluQjrbUbDvo8fBwdAEC/DBMDAHRMGAQA6JgwCADQMWEQAKBjwiAAQMeEQQCAjgmDAAAdEwYBADomDAIAdEwYBADomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB0TBgEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JgwCADQMWEQAKBjwiAAQMeEQQCAjgmDAAAdEwYBADomDAIAdEwYBADomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB0TBgEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JgwCADQMWEQAKBjwiAAQMcOH3UBq+2MM85on/nMZ0ZdxsHUSndw9NFHt61btw6hFABgVHbt2nV3a23zWh5z4sPg3XffPeoS1sTWrVvz1a9+ddRlAAArUFW3rvUxDRMDAHRMGAQA6JgwCADQsYmfMwgA9OH+++/Pnj17sm/fvlGXclBTU1PZsmVLNmzYMOpShEEAYDLs2bMnRxxxRLZu3ZqqFS/UsWpaa9m7d2/27NmTbdu2jbocw8QAwGTYt29fNm3aNNZBMEmqKps2bRqbHkxhEGDMXDI9k0umZ0ZdBqxL4x4E54xTncIgAEDHzBkEACbSzCunh7q/4z+6c0mPu+yyy/Kyl70sX//61/OUpzxlqDWsBj2DAABDtGPHjjzvec/Lzp1LC4+jJgwCAAzJD37wg3z5y1/Ohz/8YWEQAKA3l19+ec4444yceOKJOeqoo3LdddeNuqSDEgYBAIZkx44dmZ6enas4PT2dHTt2jLiig3MBCQDAEOzduzdXXnlldu/enarKAw88kKrKBRdcMFZLyTycnkEAgCH45Cc/mbPOOiu33nprbrnlltx2223Ztm1bvvSlL426tEXpGQQAJtJSl4IZlh07duSNb3zjQ9pe/vKX59JLL83zn//8Na1lOYRBAIAhuOqqq36s7Zxzzln7QpbJMDEAQMeEQQCAjgmDAGNqevdMpnfPjLoMYMIJgwAAHRMGAQA6JgwCAHTM0jIAwES6cGZ6qPs75/ilrVt4xx135Nxzz821116bRz7ykdm6dWve//7358QTTxxqPcOiZxAAYEhaaznzzDNz6qmn5qabbsqNN96Yd7/73bnzzjtHXdoB6RkEABiSz3/+89mwYUPOPvvs/W3bt28fYUUHp2cQAGBIdu/enWc+85mjLmNZhEGAEbKWIDBqwiAAwJCcdNJJ2bVr16jLWBZhEABgSE477bTcd999ufjii/e3XXvttfnCF74wwqoW5wISAGAiLXUpmGGqqlx22WU599xz8573vCdTU1P7l5YZV8IgAMAQHXPMMfnEJz4x6jKWzDAxAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JilZQCAiTTsj3rcefLxS3rcnXfemfPOOy9XX311jjzyyGzcuDFveMMbcuaZZw61nmHRMwgAMCSttbz0pS/NC17wgtx8883ZtWtXdu7cmT179oy6tAMSBgHW2CXTM7lkerg9FsB4uPLKK7Nx48acffbZ+9uOO+64vPa1rx1hVYsTBgEAhuSGG27IM57xjFGXsSzCIADAKnn1q1+dpz3taXnWs5416lIOSBgEABiSk046Kdddd93+2x/4wAdyxRVX5K677hphVYsTBgEAhuS0007Lvn37ctFFF+1vu/fee0dY0cFZWgYAmEhLXQpmmKoql19+ec4777xccMEF2bx5cx796Efnve9975rXslTCIADAED3+8Y/Pzp07R13GkhkmBgDomDAIANAxYRAAoGPCIABAx4RBAICOCYMAAB2ztAwAMJEumZ4Z6v7O2rn4uoV79+7N6aefniS54447cthhh2Xz5s1JkmuuuSYbN24caj3DIgwCAAzBpk2bcv311ydJzj///DzmMY/J61//+hFXdXCGiQEAOiYMAgB0TBgEAOiYMAgA0DFhEACgY64mBgAm0sGWgmGWMAgwBi6cmU6SnHP8zhFXAgzD+eefP+oSlswwMQBAx4RBAICOCYMAwMRorY26hCUZpzqFQQBgIkxNTWXv3r1jFbQW0lrL3r17MzU1NepSkriABACYEFu2bMmePXty1113jbqUg5qamsqWLVtGXUYSYRAAmBAbNmzItm3bRl3GumOYGACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6Nghh8GqevO87cdW1a+vpJCqOr+qXr+SfQAAsDwr6Rl887ztxyZZMAxW1WErOAYAAKtoSZ9AUlWXJ3likqkkv53kSUkeVVXXJ7khyWFJnjy4/dkkf5TkbUluT7I9yVMPsN+3JDkryW1J7kqya9B+VZKvJPnHmQ2av9pa+2JVbU3ysSSPHuziNa21/73Afl+V5FVJcuyxxy7lKQIAdGmpH0f3b1pr91TVo5Jcm+QfZTaIbU+SQUg7ed7tU5M8e9D2zYV2WFXPTDKd5OmDOq7LIAzO1dZae3ZVvTizwfKFSb6T5J+01vZV1QlJdiQ55eH7bq19KMmHkuSUU04Z70+rBgAYoaWGwXOq6szB9hOTnLCE77nmQEFw4PlJLmut3ZskVfXph93/+4N/dyXZOtjekOR3qmp7kgeSnLiEOgAAOICDhsFBL98Lkzy3tXbvYAh3agn7/tslPGaxXrv7Bv8+kAfrPC/JnUmeltn5jvuWcAwAAA5gKReQ/GSS7w6C4FOSPGfQfn9VbRhsfz/JEcs89p8lObOqHlVVRyT5hSXWcntr7e+T/Epm5yoCAHCIlhIGP5Pk8Kr6WpJ3JLl60P6hJF+rqv/RWtub5MtVtbuq3reUA7fWrkvy8STXJ/lUki8u4ds+mOQVVXV1ZoeIl9L7CADAARx0mLi1dl+Sn1/grquS/Oa8x/3LBe4/2L7fleRdC7SfOm/77gzmDLbWvpHkp+c99E0HOwYAAAfmE0gAADq21KuJD1lVbUpyxQJ3nT4YXgaAdWl690ySZOfJx+eS6dnts3YeP8qSYNlWPQwOAt/21T4OAADLZ5gYAKBjwiAADFw4M50LZ6ZHXQasKWEQAKBjwiAAQMeEQQCAjgmDABNoevfM/mVPABYjDAIAdEwYBADomDAIANAxYRAAoGPCIMA6N/PK6cy80kLJwKERBgEAOiYMAgB0TBgEAOiYMAgA0DFhEKATLjQBFiIMAgB07PBRFwDQiwd75d455P0BHDo9gwAAHRMGAQA6JgwCAHRMGAQA6JgwCLBClmxhvgtnpnPhjN8H1g9hEACgY8IgAEDHhEEAuja9eybTu2dGXQaMjDAIANAxYRDowiXTM7lkev33/rg4ARg2YRAAJsCk/MHD2hMGAQA6JgwCAHRMGAQA6JgwCACrxKfTsB4IgwCwTEIek+TwURcAwPDsv5r0raOtA1g/9AwCwAhYM5JxIQwCAD9GWO2HMAhAEvPgoFfCIABAx4RBuucjnADomTAIsAyGUsef/6O++QN/+SwtQ5fm3ijO2nn8iCthtT0YCt450jpYXV7TcOj0DAIwkfQQwdIIgwBDYngSWI+EQQAm3nJ7CUcV7K3txygIg8DITe+eyfTuyRzOM1TJOJuU8Ol1tjLCILBkhkGBQ+X9Y3wJgwAd0pPSr0nuiefQWFoGWNQ4LNkxd+LaefL4LBsyjjWxPPvD8FtHWweMmp5BAICO6RmECTE3Cfyc43eOuJL1yc9v9fjZ9s08wfGnZxCAA1rLSf+TcmUrrDd6Bpl4K+mVmDsJHv9RPRrA+tbzHEnv5YvTMwhLtNAVeGvZa2JZhodyNSyMr2Fdsay3eG3oGQTW1MF6aiflKt0Hg/s7R1oH42fud/xnp0ZcyAiNQ0/dYjWMwyoKa0kYhDG0Fm9E4/Bmtz8wvf7HA9M41Md4ckHKeHj4SIUh2PVLGGRd8Oa/uobxV/rBwtv+IV0dZWNvHHptYBTmem1fPOI61po5gzBhzC08OPMNAR6kZxDGiIDCapjrWV/N3o7lDuubBrC+6C2ebMIgXZk7KT52DcYqJ/lkN8nPbSGmKbBSa31hVG+vUVZGGIQxN3cSGdepdq6ahdUl2M3SO7l6hEGYYE4iMDxr/Xpa6A8tr2lWgzAIACM2Tgsr7w+hHX5SSa+EwQ6N61+W+4dD/9PsO1AvQwGGPpZuVPOujIA/1FJ/Z4cRKuYfS0iB1SEMcsjW4sS83Df/UQXd+RcYLDa087NTww26gmQ/VisIreR1PKzX2zj1ijEe1mPwX8/vx8IgK3awk8mahsZ11IWzkgtDrCM4HJPy0Xcs3aEEWEs+PWg9B55D1cOyScLghFrJSW7uxf7WwUeELWcfi71RDOtNZKn1rdYQn2VGxlePJ6px4/9gbfg5M0zC4AQ5lDeHpQablQTEhfQ6P3AcrOVai6MyDkOfk2Zcw4deu8lxKL9jy/34uIVe3yt5zS/1YzjH/f1EGCTJ8N/ox/XEsZj18qJdTQtO1j+E0LiWH2C/Hn/XerYep3SMu8UC0VL/MBrWiMfc8VaT1/zwVWtt1DWsqqr6TGvtjFHXsdqq6q4kt466DgBgRY5rrW1eywNOfBgEAODAHjHqAgAAGB1hEACgY8IgAEDHhEEAgI4JgwAAHRMGAQA6JgwCAHRMGAQA6JgwCADQsf8P5hZ9yHaPZ3AAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x360 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# For the comparison, extract an interval\n",
"# representing a Mafk bound region and visualize the\n",
"# important features.\n",
"gi = DNA.gindexer[7796]\n",
"chrom = gi.chrom\n",
"start = gi.start\n",
"end = gi.end\n",
"attr_mafk = input_attribution(model, DNA, chrom=chrom, start=start, end=end)\n",
"\n",
"plotGenomeTrack(SeqTrack(attr_mafk[0]), chrom, start, end)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to perform variant effect prediction, we need the DNA sequence loaded for the whole genome into a Bioseq object and a VCF file containing single nucleotide variant.\n",
"\n",
"The result of this analysis is stored in two files: scores.hdf5 and snps.bed.gz."
]
},
{
"cell_type": "code",
"execution_count": 9,
"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 = model.predict_variant_effect(DNA,\n",
" VCFFILE,\n",
" conditions=['feature'],\n",
" output_folder=vcfoutput)\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": 10,
"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": "markdown",
"metadata": {},
"source": [
"Finally, we can convert the variant predictions (the score differences in this case) along with the genomic context with other genomic tracks."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"gindexer = GenomicIndexer.create_from_file(variantsfile, None, None)\n",
"\n",
"snpcov = Cover.create_from_array('snps', f['diffscore'],\n",
" gindexer,\n",
" store_whole_genome=True,\n",
" padding_value=np.nan)\n",
"#snpcov = Cover.create_from_array('snps', f['diffscore'],\n",
"# gindexer,\n",
"# store_whole_genome=False,\n",
"# padding_value=np.nan)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAExCAYAAADx4e+wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF/9JREFUeJzt3X2wnFd9H/Dvz7aEsKGAhUIwxpbBuExNgjEKAxmgrk0mCm2Gt6TcRq2BknjcUjz2xCW8NXWamICTmRAG4qmAhpcaCWqwQ8vEJMHmtQVkGQcEtEEYCwuwscRLAGFjnNM/diVfy1e617p7d+89+/nM3NndZ599zm/v3bvPd895nrPVWgsAACvfUZMuAACA0RDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh1AB6rqHVX1B5OuA5gswQ5gylTV71fVF6rqp1V1yaTrAUZHsAOYPjuTvCLJhyZdCDBagh3ACFTVzVX1qqr6UlV9t6r+vKrWVNXDq+p/VdX3quo7VfWJqjpq+JgTqur9VXV7VX2tqi6Ytb17Da1W1VlVtXvW7SdV1Q1V9YOqem+SNQfV81tVtXPY5ger6oT997XW3tla+8skP1jK3wkwfoIdwOhsSvLLSR6b5LQkr03y20l2J1mX5BFJXp2kDcPd/0zyt0keleScJBdW1S/P10hVrU5ydZJ3Jzk+yf9I8oJZ95+d5A+T/Mskj0yyK8nWkTxDYFkT7ABG582ttVtaa99JcmmSf5XkrgzC1cmttbtaa59ogy/p/oUk61pr/6W19pPW2k1J3ppkZgHtPDXJqiRvHG7zyiTbZt2/Kcl/a63d0Fq7M8mrkjytqtaP6HkCy5RgBzA6t8y6vivJCUn+KINj2v6qqm6qqlcO7z85yQnDIdrvVdX3MujNe8QC2jkhyTeGAXF2e7PvP3C7tfbDJHsz6BkEOnbMpAsA6MijZ10/Kck3W2s/yGA49rer6vQk11XVtgxC4Ndaa487xLZ+lOTYWbd/dtb1byV5VFXVrHB3UpKvDq9/M4PgmCSpquOSrE3yjSN7WsBKoccOYHReVlUnVtXxGfS+vbeq/kVVnVpVleTvk9w9/Plskr+vqt+pqgdW1dFV9YSq+oXhtm5M8uyqOr6qfjbJhbPa+T9Jfprkgqo6pqqen+Qps+5/T5KXVNUZVfWAJK9L8pnW2s1JUlWrqmpNBvuAY4YneRy9RL8TYIwEO4DReU+Sv0py0/DnD5I8LsnfJPlhBoHsz1prH22t3Z3kV5OckeRrSfYkeVuShwy39e4MTqy4ebjN9+5vpLX2kyTPT/LiJN9N8sIkH5h1/0eS/Kck78+gd++xufexe29N8uMMjgF8zfD6vxnFLwCYrLr3IRoAHImqujnJb7bW/mbStQDTS48dAEAnBDsAgE4YigUA6IQeOwCATgh2AACdMEFxP4ypA8DKV4t5sB47AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANCJaq1NugZGoKp2JLlj0nUAAIuyprX2hCN98DGjrISJuqO1tmHSRQAAR66qrl/M4w3FAgB0QrADAOiEYNePzZMuAABYtEXtz508AQDQCT12AACdEOwAADoh2K0QVfXQqrqyqv5vVX25qp5WVcdX1V9X1VeGlw8brltV9aaq2llVn6+qMyddPwCQVNVFVfXFqtpRVVuqak1VnVJVnxnuz99bVauH6z5geHvn8P71821fsFs5/jTJNa21xyd5YpIvJ3llko+01h6X5CPD20nyK0keN/w5L8nl4y8XAJitqh6V5IIkG4aTEB+dZCbJG5L8yXB//t0kLx0+5KVJvttaOzXJnwzXOyzBbgWoqn+U5JlJ3p4krbWftNa+l+Q5Sd45XO2dSZ47vP6cJO9qA59O8tCqeuSYywYA7uuYJA+sqmOSHJvkW0nOTnLl8P6D9+f79/NXJjmnqupwGxfsVobHJLk9yZ9X1eeq6m1VdVySR7TWvpUkw8ufGa7/qCS3zHr87uEyAGBCWmvfSPLHSb6eQaD7fpLtSb7XWvvpcLXZ++wD+/Ph/d9PsvZwbQh2K8MxSc5Mcnlr7UlJfpR7hl3nMleaN68NAEzQ8Fj45yQ5JckJSY7L4PCpg+3fZ9/v/blgtzLsTrK7tfaZ4e0rMwh6t+0fYh1efnvW+o+e9fgTk3xzTLUCAHN7VpKvtdZub63dleQDSX4xg0OmjhmuM3uffWB/Prz/IUm+c7gGBLsVoLV2a5JbquofDxedk+RLST6Y5EXDZS9K8hfD6x9Mcu7w7NinJvn+/iFbAGBivp7kqVV17PBYuf378+uS/NpwnYP35/v387+W5No2zzdL+OaJFaKqzkjytiSrk9yU5CUZBPP3JTkpgxfLr7fWvjN8sbw5ycYk+5K8pLV2/UQKBwAOqKrfS/LCJD9N8rkkv5nBsXRbkxw/XPavW2t3VtWaJO9O8qQMeupmWms3HXb7gh0AQB8MxQIAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCeOmXQBS23jxo3tmmuumXQZ86nFbmCFPE8A4PAWlQm677Hbs2fPpEsYi2l5nsDKcsUVyfr1yVFHDS6vuGLSFUHfuu+xA2AyrrgiOe+8ZN++we1duwa3k2TTpsnVBT3rvscOgMl4zWvuCXX77ds3WA4sDcEOgCXx9a/fv+XA4gl2ACyJk066f8uBxRPsAFgSl16aHHvsvZcde+xgObA0BDsAlsSmTcnmzcnJJydVg8vNm504AUvJWbEALJlNmwQ5GCc9dgAAnRDsoHMmiAWYHoZioWMmiAWYLnrsoGMmiAWYLoIddMwEsQDTRbCDjpkgFmC6CHbQMRPEAkwXwQ46ZoJYgOnirFjonAliAaaHHjsAgE4IdjCF5pu02KTGACuToViYMvNNWmxSY4CVS48dTJn5Ji02qTHAyiXYwZSZb9JikxoDrFyCHUyZ+SYtNqkxwMol2MGUmW/SYpMaA6xcgh1MmfkmLTapMcDKVa21SdewpDZs2NCuv/76SZcxn1rsBlbI8wQADm9RmUCPHZDE3HUAPTCPHWDuOoBO6LEDzF0H0AnBDjB3HUAnBDvA3HUAnRDsAHPXAXRCsAPMXQfQCWfFAkkGIU6QA1jZ9NgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdGGuwq6rjqupDVfW3VbWjql5YVTdX1e9V1Q1V9YWqevxw3Uuq6t1VdW1VfaWqfmu4/JFV9fGqunG4jWeM8zkAACxX4+6x25jkm621J7bWnpDkmuHyPa21M5NcnuTiWev/fJJ/nuRpSX63qk5I8htJPtxaOyPJE5PceHAjVXVeVV1fVdfffvvtS/h0AACWj3EHuy8keVZVvaGqntFa+/5w+QeGl9uTrJ+1/l+01n7cWtuT5LokT0myLclLquqSJD/XWvvBwY201ja31ja01jasW7duqZ4LAMCyMtZg11r7uyRPziDg/WFV/e7wrjuHl3fn3l9z1u67ifbxJM9M8o0k766qc5ewZACAFWPcx9idkGRfa+2/J/njJGfO85DnVNWaqlqb5Kwk26rq5CTfbq29NcnbF7ANAICpcMz8q4zUzyX5o6r6hyR3Jfl3Sa48zPqfTfKhJCcl+f3W2jer6kVJ/mNV3ZXkh0n02AEAZMzBrrX24SQfPmjx+ln3X59Bz9x+f9daO++gbbwzyTuXqEQAgBXLPHYAAJ0Y91DsgrXWLpl0DQAAK4keOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAw7piiuS9euTo44aXF5xxaQrAuBwjpl0AcDydMUVyXnnJfv2DW7v2jW4nSSbNk2uLgAOTY8dMKfXvOaeULffvn2D5QAsT4IdMKevf/3+LQdg8gQ7YE4nnXT/lgMweYIdMKdLL02OPfbey449drAcgOVJsAPmtGlTsnlzcvLJSdXgcvNmJ04ALGfVWpt0DUuqqq5prW2cdB1LrapuT7Jr0nUAAIuyZzG5pftgBwAwLQzFAgB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATx0y6gKW2cePGds0110y6jPnUYjfw8Ic/vK1fv34EpcDAt++8KUnyMw94zIQrAZge27dv39NaW3ekj+8+2O3Zs2fSJYzF+vXrc/3110+6DDrypp0zSZILTt064UoApkdV7VrM4w3FAgB0QrADAOiEYAcA0Inuj7EDAKbDXXfdld27d+eOO+6YdCnzWrNmTU488cSsWrVqpNsV7ACALuzevTsPfvCDs379+lQtesKJJdNay969e7N79+6ccsopI922oVgAoAt33HFH1q5du6xDXZJUVdauXbskPYuCHQDQjeUe6vZbqjoFOwCATjjGDgDo0s4Xz4x0e6e+Y2ETtl911VV5/vOfny9/+ct5/OMfP9Ia5qPHDgBghLZs2ZKnP/3p2bp1/N/cI9gBAIzID3/4w3zqU5/K29/+dsEOAGAlu/rqq7Nx48acdtppOf7443PDDTeMtX3BDgBgRLZs2ZKZmcGxfTMzM9myZctY23fyBADACOzduzfXXnttduzYkarK3XffnarKZZddNrZpWPTYAQCMwJVXXplzzz03u3btys0335xbbrklp5xySj75yU+OrQY9dgBAlxY6PcmobNmyJa985SvvtewFL3hB3vOe9+QZz3jGWGoQ7AAARuCjH/3ofZZdcMEFY63BUCwAQCcEOwCATgh2AACdEOwAADoh2AEHzOzYmZkdOyddBsxp54tnRv6l7tAbwQ4AoBOmOwEAuvSmnaPt4b3g1IXNi3frrbfmwgsvzLZt2/KABzwg69evzxvf+MacdtppI61nLnrsAABGpLWW5z3veTnrrLPy1a9+NV/60pfyute9LrfddttY2tdjBwAwItddd11WrVqV888//8CyM844Y2zt67EDABiRHTt25MlPfvLE2hfsAAA6IdgBAIzI6aefnu3bt0+sfcEOgBXlXTM7864Z8y2yPJ199tm5884789a3vvXAsm3btuVjH/vYWNp38gRMgf2Tup76joWdqg/Qg4VOTzJKVZWrrroqF154YV7/+tdnzZo1B6Y7GQfBDgBghE444YS8733vm0jbhmIBADoh2AEAdEKwAwDohGAHjNzOF88cOGEDgPER7AAAOiHYAQB0wnQnAECXZnaMdiLrrU84dUHr3Xbbbbnooovy6U9/Og972MOyevXqvOIVr8jznve8kdYzlyPusauqV8+6/tCq+veLKaSqLqmqixezDQCASWqt5bnPfW6e+cxn5qabbsr27duzdevW7N69eyztL2Yo9tWzrj80yZzBrqqOXkQbAAArxrXXXpvVq1fn/PPPP7Ds5JNPzstf/vKxtL+godiqujrJo5OsSfKnSR6T5IFVdWOSLyY5Osljh7f/OsmHkvznJN9KckaSf3KI7b4myblJbklye5Ltw+UfTfKZJP8sg9D40tbaJ6pqfZJ3JzluuIn/0Fr733Ns97wk5yXJSSedtJCnCACwaF/84hdz5plnTqz9hR5j929ba9+pqgcm2Zbkn2YQqs5IkmHgesKs22clecpw2dfm2mBVPTnJTJInDeu4IcNgt7+21tpTqurZGYTEZyX5dpJfaq3dUVWPS7IlyYaDt91a25xkc5Js2LChLfA5AgCM1Mte9rJ88pOfzOrVq7Nt27Ylb2+hwe6Cqtp/xN+jkzxuAY/57KFC3dAzklzVWtuXJFX1wYPu/8DwcnuS9cPrq5K8uarOSHJ3ktMWUAcAndp/cPxCD2qHpXb66afn/e9//4Hbb3nLW7Jnz55s2HCffqglMe8xdsPet2cleVpr7YlJPpfBkOx8frSAdQ7Xm3bn8PLu3BNAL0pyW5InZtBTt3oBbQAAjMXZZ5+dO+64I5dffvmBZfv27Rtb+wvpsXtIku+21vZV1eOTPHW4/K6qWtVauyvJD5I8+H62/fEk76iq1w/r+NUk/3UBtexurf1DVb0og2P7AADuYxI9uVWVq6++OhdddFEuu+yyrFu3Lscdd1ze8IY3jKX9hQS7a5KcX1WfT/L/knx6uHxzks9X1Q2ttU1V9amq2pHkLzM4eeKwWms3VNV7k9yYZFeSTyyglj9L8v6q+vUk12VhvYIAAGPzyEc+Mlu3bp1I2/MGu9banUl+ZY67Pprkd2at9xtz3D/fti9Ncukcy8+adX1PhsfYtda+kuTnZ636qvnaAACYFr55Asi7Zoazs792snUAsDhLHuyqam2Sj8xx1zmttb1L3T6wODtfPJMkOfUdkxlWAGDhljzYDcPbGUvdDgDAtFvMV4oBI/amnTN5086ZSZcBwAol2AEAdMLJEwCwRCb9zRjTfozsgRPDRuTcrYf/O+7duzfnnHNOkuTWW2/N0UcfnXXr1iVJPvvZz2b16qX/XgXBDgBgBNauXZsbb7wxSXLJJZfkQQ96UC6++OKx1mAoFgAm4F0zO0feowSCHUwpOxWA/hiKBYAR2P9Bab7jsEZt/5n0F5y69T7Lnj3WSlgO9NgB3dEbybjsfPHMgRMU9vP6Y5IEOwCAThiKhQmZ1LDNOE37VAvc1zS87lk+pvF1JtixIsx1DAlzE6aYRkv5HiGMciQuueSSibRrKBYAoBOCHQBAJwQ7ABZtrrNDYRJaa5MuYUGWqk7H2MEKNNcxP/uX/eKaiZS0LNwTLP7giLfheKqVz3Gm02vNmjXZu3dv1q5dm6qadDmH1FrL3r17s2bN6N+wBTuYMubXAnp14oknZvfu3bn99tsnXcq81qxZkxNPPHHk2xXsAA5hZscgBG99wtL13uldgtFZtWpVTjnllEmXMVGCHQDcT44nZLly8gSMwJt2zhyYR2sUpu0riWZ27DzQO8bSmdTv+Uhez07GmG6jfk+dJnrsAJaAkzDumTT4oYs4mWWaHRzCl/KQgMO1P1+7cx1OYFL5ydFjx9jM/gQ+DZ/GF/oc5+vNWE69WePsSfSJfXIW8/+5nF6vHJlpeH/umWAHrGjj2An1uKNbjsF5OdY023I6REKA5lAMxdKtcZzRCDCt9r/HGmhfXgQ7ltxKnjh31FNROO7q3g70gl3c767BdCZzm/2/cKAXrN+XwYJ5vcxvvt/RtH+oF+w4pGn/5+DQ5gqoXi/Ly+y/x2ICdM8fRpZTiOptqH+hZg+9O9FiNAQ77mMxx22M40woZ1tNr+UQMg73dW4rPfwYWlv+RnFc3XJ9Dx1H0F5OYX6pCHZM1HLaIU7DPzz3tpxef/M5MHXIawexa9w1zzV0uhJ+bwsxzb3Nvf0tEey4n+Z6AxzVp79xvMH09ga+kt6U9wfn1w6HAw/1NzgwJPXasZR1wKTnXJvr/2i+v29vr+e5TOo5TsPvtgej+Dv19rcW7Dhik9oBH9z+qe/YqrftEAytLZ3FfKBZrkNhPVvoe8SBEzmW8H1t0h8iejTqD7krOewJdlNuGnYwyyH03XNg9JG/ka/kNxoGFvphyI6f5WoxAWr2/mbcZ0JP04dcwY6uCD8cynII+NC7SY/kzHYk+4MeJn0W7KbIwafTL+UObqX3BC506gHzby3/uejGOY/ipF73iwmti6l5uX+Qmuv3stxr5t5mvz6X+5Qwy+XDo2DXiZt+fOe9bi/0QPXlaHZX/1J9vdDsf8D9bTz7CB7LeBiaZKUbR6CcL6Tvr2ElTha/FBbzoXQ57weqtTbpGpZUVV3TWts46TqWWlXdnmTXpOsAABbl5NbauiN9cPfBDgBgWhw16QIAABgNwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0In/DxsmzW9pChKJAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 720x360 with 5 Axes>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAExCAYAAADx4e+wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF/9JREFUeJzt3X2wnFd9H/Dvz7aEsKGAhUIwxpbBuExNgjEKAxmgrk0mCm2Gt6TcRq2BknjcUjz2xCW8NXWamICTmRAG4qmAhpcaCWqwQ8vEJMHmtQVkGQcEtEEYCwuwscRLAGFjnNM/diVfy1e617p7d+89+/nM3NndZ599zm/v3bvPd895nrPVWgsAACvfUZMuAACA0RDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh1AB6rqHVX1B5OuA5gswQ5gylTV71fVF6rqp1V1yaTrAUZHsAOYPjuTvCLJhyZdCDBagh3ACFTVzVX1qqr6UlV9t6r+vKrWVNXDq+p/VdX3quo7VfWJqjpq+JgTqur9VXV7VX2tqi6Ytb17Da1W1VlVtXvW7SdV1Q1V9YOqem+SNQfV81tVtXPY5ger6oT997XW3tla+8skP1jK3wkwfoIdwOhsSvLLSR6b5LQkr03y20l2J1mX5BFJXp2kDcPd/0zyt0keleScJBdW1S/P10hVrU5ydZJ3Jzk+yf9I8oJZ95+d5A+T/Mskj0yyK8nWkTxDYFkT7ABG582ttVtaa99JcmmSf5XkrgzC1cmttbtaa59ogy/p/oUk61pr/6W19pPW2k1J3ppkZgHtPDXJqiRvHG7zyiTbZt2/Kcl/a63d0Fq7M8mrkjytqtaP6HkCy5RgBzA6t8y6vivJCUn+KINj2v6qqm6qqlcO7z85yQnDIdrvVdX3MujNe8QC2jkhyTeGAXF2e7PvP3C7tfbDJHsz6BkEOnbMpAsA6MijZ10/Kck3W2s/yGA49rer6vQk11XVtgxC4Ndaa487xLZ+lOTYWbd/dtb1byV5VFXVrHB3UpKvDq9/M4PgmCSpquOSrE3yjSN7WsBKoccOYHReVlUnVtXxGfS+vbeq/kVVnVpVleTvk9w9/Plskr+vqt+pqgdW1dFV9YSq+oXhtm5M8uyqOr6qfjbJhbPa+T9Jfprkgqo6pqqen+Qps+5/T5KXVNUZVfWAJK9L8pnW2s1JUlWrqmpNBvuAY4YneRy9RL8TYIwEO4DReU+Sv0py0/DnD5I8LsnfJPlhBoHsz1prH22t3Z3kV5OckeRrSfYkeVuShwy39e4MTqy4ebjN9+5vpLX2kyTPT/LiJN9N8sIkH5h1/0eS/Kck78+gd++xufexe29N8uMMjgF8zfD6vxnFLwCYrLr3IRoAHImqujnJb7bW/mbStQDTS48dAEAnBDsAgE4YigUA6IQeOwCATgh2AACdMEFxP4ypA8DKV4t5sB47AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANCJaq1NugZGoKp2JLlj0nUAAIuyprX2hCN98DGjrISJuqO1tmHSRQAAR66qrl/M4w3FAgB0QrADAOiEYNePzZMuAABYtEXtz508AQDQCT12AACdEOwAADoh2K0QVfXQqrqyqv5vVX25qp5WVcdX1V9X1VeGlw8brltV9aaq2llVn6+qMyddPwCQVNVFVfXFqtpRVVuqak1VnVJVnxnuz99bVauH6z5geHvn8P71821fsFs5/jTJNa21xyd5YpIvJ3llko+01h6X5CPD20nyK0keN/w5L8nl4y8XAJitqh6V5IIkG4aTEB+dZCbJG5L8yXB//t0kLx0+5KVJvttaOzXJnwzXOyzBbgWoqn+U5JlJ3p4krbWftNa+l+Q5Sd45XO2dSZ47vP6cJO9qA59O8tCqeuSYywYA7uuYJA+sqmOSHJvkW0nOTnLl8P6D9+f79/NXJjmnqupwGxfsVobHJLk9yZ9X1eeq6m1VdVySR7TWvpUkw8ufGa7/qCS3zHr87uEyAGBCWmvfSPLHSb6eQaD7fpLtSb7XWvvpcLXZ++wD+/Ph/d9PsvZwbQh2K8MxSc5Mcnlr7UlJfpR7hl3nMleaN68NAEzQ8Fj45yQ5JckJSY7L4PCpg+3fZ9/v/blgtzLsTrK7tfaZ4e0rMwh6t+0fYh1efnvW+o+e9fgTk3xzTLUCAHN7VpKvtdZub63dleQDSX4xg0OmjhmuM3uffWB/Prz/IUm+c7gGBLsVoLV2a5JbquofDxedk+RLST6Y5EXDZS9K8hfD6x9Mcu7w7NinJvn+/iFbAGBivp7kqVV17PBYuf378+uS/NpwnYP35/v387+W5No2zzdL+OaJFaKqzkjytiSrk9yU5CUZBPP3JTkpgxfLr7fWvjN8sbw5ycYk+5K8pLV2/UQKBwAOqKrfS/LCJD9N8rkkv5nBsXRbkxw/XPavW2t3VtWaJO9O8qQMeupmWms3HXb7gh0AQB8MxQIAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCeOmXQBS23jxo3tmmuumXQZ86nFbmCFPE8A4PAWlQm677Hbs2fPpEsYi2l5nsDKcsUVyfr1yVFHDS6vuGLSFUHfuu+xA2AyrrgiOe+8ZN++we1duwa3k2TTpsnVBT3rvscOgMl4zWvuCXX77ds3WA4sDcEOgCXx9a/fv+XA4gl2ACyJk066f8uBxRPsAFgSl16aHHvsvZcde+xgObA0BDsAlsSmTcnmzcnJJydVg8vNm504AUvJWbEALJlNmwQ5GCc9dgAAnRDsoHMmiAWYHoZioWMmiAWYLnrsoGMmiAWYLoIddMwEsQDTRbCDjpkgFmC6CHbQMRPEAkwXwQ46ZoJYgOnirFjonAliAaaHHjsAgE4IdjCF5pu02KTGACuToViYMvNNWmxSY4CVS48dTJn5Ji02qTHAyiXYwZSZb9JikxoDrFyCHUyZ+SYtNqkxwMol2MGUmW/SYpMaA6xcgh1MmfkmLTapMcDKVa21SdewpDZs2NCuv/76SZcxn1rsBlbI8wQADm9RmUCPHZDE3HUAPTCPHWDuOoBO6LEDzF0H0AnBDjB3HUAnBDvA3HUAnRDsAHPXAXRCsAPMXQfQCWfFAkkGIU6QA1jZ9NgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdGGuwq6rjqupDVfW3VbWjql5YVTdX1e9V1Q1V9YWqevxw3Uuq6t1VdW1VfaWqfmu4/JFV9fGqunG4jWeM8zkAACxX4+6x25jkm621J7bWnpDkmuHyPa21M5NcnuTiWev/fJJ/nuRpSX63qk5I8htJPtxaOyPJE5PceHAjVXVeVV1fVdfffvvtS/h0AACWj3EHuy8keVZVvaGqntFa+/5w+QeGl9uTrJ+1/l+01n7cWtuT5LokT0myLclLquqSJD/XWvvBwY201ja31ja01jasW7duqZ4LAMCyMtZg11r7uyRPziDg/WFV/e7wrjuHl3fn3l9z1u67ifbxJM9M8o0k766qc5ewZACAFWPcx9idkGRfa+2/J/njJGfO85DnVNWaqlqb5Kwk26rq5CTfbq29NcnbF7ANAICpcMz8q4zUzyX5o6r6hyR3Jfl3Sa48zPqfTfKhJCcl+f3W2jer6kVJ/mNV3ZXkh0n02AEAZMzBrrX24SQfPmjx+ln3X59Bz9x+f9daO++gbbwzyTuXqEQAgBXLPHYAAJ0Y91DsgrXWLpl0DQAAK4keOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAw7piiuS9euTo44aXF5xxaQrAuBwjpl0AcDydMUVyXnnJfv2DW7v2jW4nSSbNk2uLgAOTY8dMKfXvOaeULffvn2D5QAsT4IdMKevf/3+LQdg8gQ7YE4nnXT/lgMweYIdMKdLL02OPfbey449drAcgOVJsAPmtGlTsnlzcvLJSdXgcvNmJ04ALGfVWpt0DUuqqq5prW2cdB1LrapuT7Jr0nUAAIuyZzG5pftgBwAwLQzFAgB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATgh0AQCcEOwCATgh2AACdEOwAADoh2AEAdEKwAwDohGAHANAJwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0AnBDgCgE4IdAEAnBDsAgE4IdgAAnRDsAAA6IdgBAHRCsAMA6IRgBwDQCcEOAKATx0y6gKW2cePGds0110y6jPnUYjfw8Ic/vK1fv34EpcDAt++8KUnyMw94zIQrAZge27dv39NaW3ekj+8+2O3Zs2fSJYzF+vXrc/3110+6DDrypp0zSZILTt064UoApkdV7VrM4w3FAgB0QrADAOiEYAcA0Inuj7EDAKbDXXfdld27d+eOO+6YdCnzWrNmTU488cSsWrVqpNsV7ACALuzevTsPfvCDs379+lQtesKJJdNay969e7N79+6ccsopI922oVgAoAt33HFH1q5du6xDXZJUVdauXbskPYuCHQDQjeUe6vZbqjoFOwCATjjGDgDo0s4Xz4x0e6e+Y2ETtl911VV5/vOfny9/+ct5/OMfP9Ia5qPHDgBghLZs2ZKnP/3p2bp1/N/cI9gBAIzID3/4w3zqU5/K29/+dsEOAGAlu/rqq7Nx48acdtppOf7443PDDTeMtX3BDgBgRLZs2ZKZmcGxfTMzM9myZctY23fyBADACOzduzfXXnttduzYkarK3XffnarKZZddNrZpWPTYAQCMwJVXXplzzz03u3btys0335xbbrklp5xySj75yU+OrQY9dgBAlxY6PcmobNmyJa985SvvtewFL3hB3vOe9+QZz3jGWGoQ7AAARuCjH/3ofZZdcMEFY63BUCwAQCcEOwCATgh2AACdEOwAADoh2AEHzOzYmZkdOyddBsxp54tnRv6l7tAbwQ4AoBOmOwEAuvSmnaPt4b3g1IXNi3frrbfmwgsvzLZt2/KABzwg69evzxvf+MacdtppI61nLnrsAABGpLWW5z3veTnrrLPy1a9+NV/60pfyute9LrfddttY2tdjBwAwItddd11WrVqV888//8CyM844Y2zt67EDABiRHTt25MlPfvLE2hfsAAA6IdgBAIzI6aefnu3bt0+sfcEOgBXlXTM7864Z8y2yPJ199tm5884789a3vvXAsm3btuVjH/vYWNp38gRMgf2Tup76joWdqg/Qg4VOTzJKVZWrrroqF154YV7/+tdnzZo1B6Y7GQfBDgBghE444YS8733vm0jbhmIBADoh2AEAdEKwAwDohGAHjNzOF88cOGEDgPER7AAAOiHYAQB0wnQnAECXZnaMdiLrrU84dUHr3Xbbbbnooovy6U9/Og972MOyevXqvOIVr8jznve8kdYzlyPusauqV8+6/tCq+veLKaSqLqmqixezDQCASWqt5bnPfW6e+cxn5qabbsr27duzdevW7N69eyztL2Yo9tWzrj80yZzBrqqOXkQbAAArxrXXXpvVq1fn/PPPP7Ds5JNPzstf/vKxtL+godiqujrJo5OsSfKnSR6T5IFVdWOSLyY5Osljh7f/OsmHkvznJN9KckaSf3KI7b4myblJbklye5Ltw+UfTfKZJP8sg9D40tbaJ6pqfZJ3JzluuIn/0Fr733Ns97wk5yXJSSedtJCnCACwaF/84hdz5plnTqz9hR5j929ba9+pqgcm2Zbkn2YQqs5IkmHgesKs22clecpw2dfm2mBVPTnJTJInDeu4IcNgt7+21tpTqurZGYTEZyX5dpJfaq3dUVWPS7IlyYaDt91a25xkc5Js2LChLfA5AgCM1Mte9rJ88pOfzOrVq7Nt27Ylb2+hwe6Cqtp/xN+jkzxuAY/57KFC3dAzklzVWtuXJFX1wYPu/8DwcnuS9cPrq5K8uarOSHJ3ktMWUAcAndp/cPxCD2qHpXb66afn/e9//4Hbb3nLW7Jnz55s2HCffqglMe8xdsPet2cleVpr7YlJPpfBkOx8frSAdQ7Xm3bn8PLu3BNAL0pyW5InZtBTt3oBbQAAjMXZZ5+dO+64I5dffvmBZfv27Rtb+wvpsXtIku+21vZV1eOTPHW4/K6qWtVauyvJD5I8+H62/fEk76iq1w/r+NUk/3UBtexurf1DVb0og2P7AADuYxI9uVWVq6++OhdddFEuu+yyrFu3Lscdd1ze8IY3jKX9hQS7a5KcX1WfT/L/knx6uHxzks9X1Q2ttU1V9amq2pHkLzM4eeKwWms3VNV7k9yYZFeSTyyglj9L8v6q+vUk12VhvYIAAGPzyEc+Mlu3bp1I2/MGu9banUl+ZY67Pprkd2at9xtz3D/fti9Ncukcy8+adX1PhsfYtda+kuTnZ636qvnaAACYFr55Asi7Zoazs792snUAsDhLHuyqam2Sj8xx1zmttb1L3T6wODtfPJMkOfUdkxlWAGDhljzYDcPbGUvdDgDAtFvMV4oBI/amnTN5086ZSZcBwAol2AEAdMLJEwCwRCb9zRjTfozsgRPDRuTcrYf/O+7duzfnnHNOkuTWW2/N0UcfnXXr1iVJPvvZz2b16qX/XgXBDgBgBNauXZsbb7wxSXLJJZfkQQ96UC6++OKx1mAoFgAm4F0zO0feowSCHUwpOxWA/hiKBYAR2P9Bab7jsEZt/5n0F5y69T7Lnj3WSlgO9NgB3dEbybjsfPHMgRMU9vP6Y5IEOwCAThiKhQmZ1LDNOE37VAvc1zS87lk+pvF1JtixIsx1DAlzE6aYRkv5HiGMciQuueSSibRrKBYAoBOCHQBAJwQ7ABZtrrNDYRJaa5MuYUGWqk7H2MEKNNcxP/uX/eKaiZS0LNwTLP7giLfheKqVz3Gm02vNmjXZu3dv1q5dm6qadDmH1FrL3r17s2bN6N+wBTuYMubXAnp14oknZvfu3bn99tsnXcq81qxZkxNPPHHk2xXsAA5hZscgBG99wtL13uldgtFZtWpVTjnllEmXMVGCHQDcT44nZLly8gSMwJt2zhyYR2sUpu0riWZ27DzQO8bSmdTv+Uhez07GmG6jfk+dJnrsAJaAkzDumTT4oYs4mWWaHRzCl/KQgMO1P1+7cx1OYFL5ydFjx9jM/gQ+DZ/GF/oc5+vNWE69WePsSfSJfXIW8/+5nF6vHJlpeH/umWAHrGjj2An1uKNbjsF5OdY023I6REKA5lAMxdKtcZzRCDCt9r/HGmhfXgQ7ltxKnjh31FNROO7q3g70gl3c767BdCZzm/2/cKAXrN+XwYJ5vcxvvt/RtH+oF+w4pGn/5+DQ5gqoXi/Ly+y/x2ICdM8fRpZTiOptqH+hZg+9O9FiNAQ77mMxx22M40woZ1tNr+UQMg73dW4rPfwYWlv+RnFc3XJ9Dx1H0F5OYX6pCHZM1HLaIU7DPzz3tpxef/M5MHXIawexa9w1zzV0uhJ+bwsxzb3Nvf0tEey4n+Z6AxzVp79xvMH09ga+kt6U9wfn1w6HAw/1NzgwJPXasZR1wKTnXJvr/2i+v29vr+e5TOo5TsPvtgej+Dv19rcW7Dhik9oBH9z+qe/YqrftEAytLZ3FfKBZrkNhPVvoe8SBEzmW8H1t0h8iejTqD7krOewJdlNuGnYwyyH03XNg9JG/ka/kNxoGFvphyI6f5WoxAWr2/mbcZ0JP04dcwY6uCD8cynII+NC7SY/kzHYk+4MeJn0W7KbIwafTL+UObqX3BC506gHzby3/uejGOY/ipF73iwmti6l5uX+Qmuv3stxr5t5mvz6X+5Qwy+XDo2DXiZt+fOe9bi/0QPXlaHZX/1J9vdDsf8D9bTz7CB7LeBiaZKUbR6CcL6Tvr2ElTha/FBbzoXQ57weqtTbpGpZUVV3TWts46TqWWlXdnmTXpOsAABbl5NbauiN9cPfBDgBgWhw16QIAABgNwQ4AoBOCHQBAJwQ7AIBOCHYAAJ0Q7AAAOiHYAQB0QrADAOiEYAcA0In/DxsmzW9pChKJAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 720x360 with 5 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"gi = DNA.gindexer[3]\n",
"chrom = gi.chrom\n",
"start = gi.start\n",
"end = gi.end\n",
"\n",
"plotGenomeTrack([LineTrack(snpcov,\n",
" linestyle=\"None\"), SeqTrack(attr_oct[0])],\n",
" chrom, start, end)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The score difference shows a dip around the site that is indicated as most important from the input attribution as well."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is also possible to export the variant effect predictions as bigwig for further explorations in e.g. IGV.\n",
"To this end, use the `export_to_bigwig` method"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"os.makedirs('./snps', exist_ok=True)\n",
"snpcov.export_to_bigwig('./snps')"
]
}
],
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}