Genomic Grammar Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Genomic Grammar Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jenhan/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:34: 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": [
    "from __future__ import print_function\n",
    "import keras\n",
    "from keras.datasets import mnist\n",
    "from keras.models import Sequential\n",
    "from keras.layers import *\n",
    "from keras.models import Model\n",
    "from keras import backend as K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import Bio\n",
    "from Bio import SeqIO\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import Bio.motifs\n",
    "%matplotlib inline\n",
    "from sklearn import model_selection\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "import sklearn\n",
    "from IPython.display import SVG\n",
    "from keras.utils.vis_utils import model_to_dot\n",
    "import scipy\n",
    "sns.set_context('notebook')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if not os.path.isdir('/home/jtao/analysis/genomic_grammar_analysis/'):\n",
    "    os.mkdir('/home/jtao/analysis/genomic_grammar_analysis')\n",
    "os.chdir('/home/jtao/analysis/genomic_grammar_analysis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def convert_sequences_to_array(sequences):\n",
    "    '''\n",
    "    inputs: sequence of nucleotides represented as a string composed of A, C, G, T\n",
    "    outputs: a list of numpy array representations of a sequence with:\n",
    "             A = [1, 0, 0, 0]\n",
    "             C = [0, 1, 0, 0]\n",
    "             G = [0, 0, 1, 0]\n",
    "             T = [0, 0, 0, 1]\n",
    "             \n",
    "    '''\n",
    "\n",
    "    nucleotide_array_dict = {'A': [1, 0, 0, 0],\n",
    "                             'C': [0, 1, 0, 0],\n",
    "                             'G': [0, 0, 1, 0],\n",
    "                             'T': [0, 0, 0, 1],\n",
    "                             'N': [0.25,0.25,0.25,0.25]}\n",
    "\n",
    "    sequence_array_list = []\n",
    "    for seq in sequences:\n",
    "        seq_array = []\n",
    "        for nuc in seq:\n",
    "            seq_array.append(nucleotide_array_dict[nuc])\n",
    "        seq_array = np.array(seq_array)\n",
    "        sequence_array_list.append(seq_array)\n",
    "    sequence_array_list = np.array(sequence_array_list)\n",
    "    return sequence_array_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def quantile_normalize_df(df_input):\n",
    "    df = df_input.copy()\n",
    "    #compute rank\n",
    "    dic = {}\n",
    "    for col in df:\n",
    "        dic.update({col : sorted(df[col])})\n",
    "    sorted_df = pd.DataFrame(dic)\n",
    "    rank = sorted_df.mean(axis = 1).tolist()\n",
    "    #sort\n",
    "    for col in df:\n",
    "        t = np.searchsorted(np.sort(df[col]), df[col])\n",
    "        df[col] = [rank[i] for i in t]\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract Sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "if [ ! -d ./bed_files ]; then mkdir ./bed_files; else rm ./bed_files/*; fi\n",
    "for peak in ./atac_idr_peaks/*tsv;\n",
    "do echo $peak;\n",
    "new_path=${peak/atac_idr_peaks/bed_files};\n",
    "new_path=${new_path/.tsv/.bed};\n",
    "echo $new_path;\n",
    "pos2bed.pl $peak >$new_path\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "if [ ! -d ./fasta_files ]; then mkdir ./fasta_files; else rm ./fasta_files/*; fi\n",
    "for bed_path in ./bed_files/*bed;\n",
    "do echo $bed_path;\n",
    "new_path=${bed_path/bed_files/fasta_files};\n",
    "new_path=${new_path/.bed/.fa};\n",
    "echo $new_path;\n",
    "/home/jtao/code/tba/model_training/extract_sequences.py $bed_path mm10 $new_path\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Background Coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "if [ ! -d ./background_files ]; then mkdir ./background_files; else rm ./background_files/*; fi\n",
    "for bed_path in ./bed_files/*bed;\n",
    "do echo $bed_path;\n",
    "echo /home/jtao/code/tba/model_training/generate_background_coordinates.py $bed_path ./background_files -genome mm10\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "rm ./make_background.sh\n",
    "script_path=\"./make_background.sh\"\n",
    "if [ ! -d ./background_files/ ]; then mkdir ./background_files/ ; fi\n",
    "for i in ./bed_files/*bed;\n",
    "do \n",
    "    factor=${i##*/};\n",
    "    factor=${factor%.bed};\n",
    "    fasta_path=\"./background_files/${factor}_background.fasta\"\n",
    "    bed_path=\"./background_files/${factor}_background.bed\"\n",
    "\n",
    "    echo \"/home/jtao/code/tba/model_training/generate_background_coordinates.py $i ./background_files/ -genome mm10\" >> $script_path;\n",
    "    echo \"mv ./background_files/background.bed $bed_path\" >> $script_path;\n",
    "    echo \"mv ./background_files/background.fasta $fasta_path\" >> $script_path;\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "rm ./background/*\n",
    "chmod a+x ./*sh\n",
    "bash ./make_background.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# positive_seqRecords = list(SeqIO.parse('./fasta_files/c57bl6_il4-24h_peaks.fa', 'fasta'))\n",
    "# negative_seqRecords = list(SeqIO.parse('./background_files/c57bl6_il4-24h_peaks_background.fasta', 'fasta'))\n",
    "\n",
    "positive_seqRecords = list(SeqIO.parse('/home/jtao/analysis/ap1_fdr_analysis/fasta_files/c57bl6_atf3_veh_idr.fasta', 'fasta'))\n",
    "negative_seqRecords = list(SeqIO.parse('./background_files/c57bl6_il4-24h_background.fasta', 'fasta'))[:len(positive_seqRecords)]\n",
    "\n",
    "fasta_seq = [str(x.seq[:200]) for x in positive_seqRecords] + [str(x[:200].seq) for x in negative_seqRecords]\n",
    "\n",
    "fasta_rc_seq = [str(x[:200].reverse_complement().seq) for x in positive_seqRecords] + \\\n",
    "    [str(x[:200].reverse_complement().seq) for x in negative_seqRecords]\n",
    "\n",
    "sequence_arrays = convert_sequences_to_array(fasta_seq)\n",
    "sequence_arrays = np.array(sequence_arrays)\n",
    "\n",
    "sequence_rc_arrays = convert_sequences_to_array(fasta_rc_seq)\n",
    "sequence_rc_arrays = np.array(sequence_rc_arrays)\n",
    "\n",
    "\n",
    "labels = [1 for x in positive_seqRecords] + [0 for x in negative_seqRecords]\n",
    "labels = np.array(labels)\n",
    "\n",
    "x_train, x_test, x_rc_train, x_rc_test, y_train, y_test = model_selection.train_test_split(sequence_arrays, sequence_rc_arrays, labels, test_size=0.2)\n",
    "\n",
    "num_classes = 2\n",
    "y_train = keras.utils.to_categorical(y_train, num_classes)\n",
    "y_test = keras.utils.to_categorical(y_test, num_classes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Count Params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seq_size = 150\n",
    "num_classes = 2\n",
    "num_motifs = 50\n",
    "motif_size = 20\n",
    "adjacent_bp_pool_size = 10\n",
    "attention_dim = 50 # 350 from A Structured Self-attentive Sentence Embedding\n",
    "attention_hops = 1 # from A Structured Self-attentive Sentence Embedding\n",
    "num_dense_neurons = 500 # 2-layer, 2000 units, from A Structured Self-attentive Sentence Embedding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seq_size = 150\n",
    "num_classes = 2\n",
    "num_motifs = 100\n",
    "motif_size = 20\n",
    "adjacent_bp_pool_size = 1\n",
    "attention_dim = 50 \n",
    "attention_hops = 1 \n",
    "num_dense_neurons = 500 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Top Motif Scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_seq_length = len(fasta_seq[0])\n",
    "\n",
    "input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "# find motifs\n",
    "convolution_layer = Conv1D(filters=num_motifs, \n",
    "    kernel_size=motif_size,\n",
    "    activation='relu',\n",
    "    input_shape=(total_seq_length,4),\n",
    "    name='convolution_layer',\n",
    "    padding = 'same'\n",
    "    )\n",
    "forward_motif_scores = convolution_layer(input_fwd)\n",
    "reverse_motif_scores = convolution_layer(input_rev)\n",
    "print('forward_motif_scores', forward_motif_scores.shape)\n",
    "\n",
    "# crop motif scores to avoid parts of sequence where motif score is computed in only one direction\n",
    "to_crop = int((total_seq_length - seq_size)/2)\n",
    "crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "    name='crop_layer')\n",
    "cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "print('cropped_fwd_scores', cropped_fwd_scores.shape)\n",
    "\n",
    "# # flip motif scores\n",
    "# flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "#     output_shape=(seq_size, num_motifs),\n",
    "#     name='flip_layer')\n",
    "# flipped_rev_scores = flip_layer(cropped_rev_scores)\n",
    "# print('flipped_rev_scores', flipped_rev_scores.shape)\n",
    "\n",
    "# calculate max scores for each orientation\n",
    "seq_pool_layer = MaxPool1D(pool_size=seq_size)\n",
    "max_fwd_scores = seq_pool_layer(cropped_fwd_scores)\n",
    "max_rev_scores = seq_pool_layer(cropped_rev_scores)\n",
    "print('max_fwd_scores', max_fwd_scores.shape)\n",
    "\n",
    "# calculate max score for strand\n",
    "orientation_max_layer = Maximum()\n",
    "max_seq_scores = orientation_max_layer([max_fwd_scores, max_rev_scores])\n",
    "print('max_seq_scores', max_seq_scores.shape)\n",
    "\n",
    "# fully connected layer\n",
    "dense_out = Dense(num_dense_neurons, activation='relu', \n",
    "                 )(max_seq_scores)\n",
    "\n",
    "# drop out\n",
    "drop_out = Dropout(0.25)(dense_out)\n",
    "\n",
    "# make prediction\n",
    "flattened = Flatten()(drop_out)\n",
    "predictions = Dense(num_classes,\n",
    "                    activation = 'softmax', \n",
    "                   )(flattened)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define and compile model\n",
    "convolution_model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "convolution_model.compile(loss=keras.losses.categorical_crossentropy,\n",
    "              optimizer=keras.optimizers.Adam(),\n",
    "              metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SVG(model_to_dot(convolution_model, show_shapes=True, show_layer_names=True).create(prog='dot', format='svg'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "convolution_model.fit([x_train, x_rc_train], y_train,\n",
    "          batch_size=100,\n",
    "          epochs=10,\n",
    "          verbose=1,\n",
    "          validation_data=([x_test, x_rc_test], y_test))\n",
    "score = convolution_model.evaluate([x_test, x_rc_test], y_test, verbose=0)\n",
    "print('Test loss:', score[0])\n",
    "print('Test accuracy:', score[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "probs = convolution_model.predict([x_test, x_rc_test])\n",
    "\n",
    "sklearn.metrics.roc_auc_score([y[1] for y in y_test], probs[:,1], )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attention Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seq_size = 150\n",
    "num_classes = 2\n",
    "num_motifs = 50\n",
    "motif_size = 20\n",
    "adjacent_bp_pool_size = 10\n",
    "attention_dim = 350 \n",
    "attention_hops = 1 \n",
    "num_dense_neurons = 1000 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_seq_length = len(fasta_seq[0])\n",
    "\n",
    "input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "# find motifs\n",
    "convolution_layer = Conv1D(filters=num_motifs, \n",
    "    kernel_size=motif_size,\n",
    "    activation='relu',\n",
    "    input_shape=(total_seq_length,4),\n",
    "    name='convolution_layer',\n",
    "    padding = 'same'\n",
    "    )\n",
    "forward_motif_scores = convolution_layer(input_fwd)\n",
    "reverse_motif_scores = convolution_layer(input_rev)\n",
    "print('forward_motif_scores', forward_motif_scores.shape)\n",
    "\n",
    "# crop motif scores to avoid parts of sequence where motif score is computed in only one direction\n",
    "to_crop = int((total_seq_length - seq_size)/2)\n",
    "crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "    name='crop_layer')\n",
    "cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "print('cropped_fwd_scores', cropped_fwd_scores.shape)\n",
    "\n",
    "# flip motif scores\n",
    "flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "    output_shape=(seq_size, num_motifs),\n",
    "    name='flip_layer')\n",
    "flipped_rev_scores = flip_layer(cropped_rev_scores)\n",
    "print('flipped_rev_scores', flipped_rev_scores.shape)\n",
    "\n",
    "# concatenate motif scores\n",
    "concatenate_layer = keras.layers.Concatenate(axis=2, name='concatenate_layer')\n",
    "concatenated_motif_scores = concatenate_layer([cropped_fwd_scores, flipped_rev_scores])\n",
    "print('concatenated_motif_scores', concatenated_motif_scores.shape)\n",
    "\n",
    "# pool across length of sequence\n",
    "sequence_pooling_layer = MaxPool1D(pool_size=adjacent_bp_pool_size, \n",
    "    strides=adjacent_bp_pool_size,\n",
    "    name='sequence_pooling_layer')\n",
    "pooled_scores = sequence_pooling_layer(concatenated_motif_scores)\n",
    "print('pooled_scores', pooled_scores.shape)\n",
    "\n",
    "# bidirectional LSTM\n",
    "forward_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "    return_sequences=True,\n",
    "    name = 'forward_lstm_layer'\n",
    "    )\n",
    "forward_hidden_states = forward_lstm_layer(pooled_scores)\n",
    "print('forward_hidden_states', forward_hidden_states.shape)\n",
    "\n",
    "reverse_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "    return_sequences=True,\n",
    "    name = 'reverse_lstm_layer',\n",
    "    go_backwards=True\n",
    "    )\n",
    "reverse_hidden_states = reverse_lstm_layer(pooled_scores)\n",
    "print('reverse_hidden_states', reverse_hidden_states.shape)\n",
    "\n",
    "# concatenate lstm hidden states\n",
    "lstm_concatenate_layer = Concatenate(axis=2)\n",
    "bilstm_hidden_states = lstm_concatenate_layer([forward_hidden_states, reverse_hidden_states])\n",
    "print('bilstm_hidden_states', bilstm_hidden_states.shape)\n",
    "\n",
    "# fully connected layer\n",
    "dense_layer = Dense(num_dense_neurons, \n",
    "    activation='relu', \n",
    "    name = 'dense_layer'\n",
    "    )\n",
    "\n",
    "dense_output = dense_layer(bilstm_hidden_states)\n",
    "\n",
    "# drop out\n",
    "drop_out = Dropout(0.25,name='dense_dropout')(dense_output)\n",
    "\n",
    "# make prediction\n",
    "flattened = Flatten(name='flatten')(drop_out)\n",
    "predictions = Dense(num_classes,\n",
    "                    name='predictions',\n",
    "                    activation = 'softmax', \n",
    "                   )(flattened)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define and compile model\n",
    "bilstm_model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "bilstm_model.compile(loss=keras.losses.categorical_crossentropy,\n",
    "              optimizer=keras.optimizers.Adam(),\n",
    "              metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SVG(model_to_dot(bilstm_model, show_shapes=True, show_layer_names=True).create(prog='dot', format='svg'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "bilstm_model.fit([x_train, x_rc_train], y_train,\n",
    "          batch_size=100,\n",
    "          epochs=10,\n",
    "          verbose=1,\n",
    "          validation_data=([x_test, x_rc_test], y_test))\n",
    "score = bilstm_model.evaluate([x_test, x_rc_test], y_test, verbose=0)\n",
    "print('Test loss:', score[0])\n",
    "print('Test accuracy:', score[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "probs = bilstm_model.predict([x_test, x_rc_test])\n",
    "\n",
    "sklearn.metrics.roc_auc_score([y[1] for y in y_test], probs[:,1], )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attention Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seq_size = 150\n",
    "num_classes = 2\n",
    "num_motifs = 25\n",
    "motif_size = 12\n",
    "adjacent_bp_pool_size = 1\n",
    "attention_dim = 200\n",
    "attention_hops = 1\n",
    "num_dense_neurons = 500 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_seq_length = len(fasta_seq[0])\n",
    "\n",
    "input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "### find motifs ###\n",
    "convolution_layer = Conv1D(filters=num_motifs, \n",
    "    kernel_size=motif_size,\n",
    "    activation='relu',\n",
    "    input_shape=(total_seq_length,4),\n",
    "    name='convolution_layer',\n",
    "    padding = 'same'\n",
    "    )\n",
    "forward_motif_scores = convolution_layer(input_fwd)\n",
    "reverse_motif_scores = convolution_layer(input_rev)\n",
    "print('forward_motif_scores', forward_motif_scores.get_shape())\n",
    "\n",
    "### crop motif scores to avoid parts of sequence where motif score is computed in only one direction ###\n",
    "to_crop = int((total_seq_length - seq_size)/2)\n",
    "crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "    name='crop_layer')\n",
    "cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "print('cropped_fwd_scores', cropped_fwd_scores.get_shape())\n",
    "\n",
    "### flip motif scores ###\n",
    "flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "    output_shape=(seq_size, num_motifs),\n",
    "    name='flip_layer')\n",
    "flipped_rev_scores = flip_layer(cropped_rev_scores)\n",
    "print('flipped_rev_scores', flipped_rev_scores.get_shape())\n",
    "\n",
    "### concatenate motif scores ###\n",
    "concatenate_layer = keras.layers.Concatenate(axis=2, name='concatenate_layer')\n",
    "concatenated_motif_scores = concatenate_layer([cropped_fwd_scores, flipped_rev_scores])\n",
    "print('concatenated_motif_scores', concatenated_motif_scores.get_shape())\n",
    "\n",
    "### pool across length of sequence ###\n",
    "sequence_pooling_layer = MaxPool1D(pool_size=adjacent_bp_pool_size, \n",
    "    strides=adjacent_bp_pool_size,\n",
    "    name='sequence_pooling_layer')\n",
    "pooled_scores = sequence_pooling_layer(concatenated_motif_scores)\n",
    "print('pooled_scores', pooled_scores.get_shape())\n",
    "\n",
    "## bidirectional LSTM ###\n",
    "forward_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "    return_sequences=True,\n",
    "    input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "    name = 'forward_lstm_layer'\n",
    "    )\n",
    "forward_hidden_states = forward_lstm_layer(pooled_scores)\n",
    "print('forward_hidden_states', forward_hidden_states.get_shape())\n",
    "\n",
    "reverse_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "    return_sequences=True,\n",
    "    input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "    name = 'reverse_lstm_layer',\n",
    "    go_backwards=True,\n",
    "    )\n",
    "reverse_hidden_states = reverse_lstm_layer(pooled_scores)\n",
    "print('reverse_hidden_states', reverse_hidden_states.get_shape())\n",
    "### concatenate lstm hidden states ###\n",
    "lstm_concatenate_layer = Concatenate(axis=2)\n",
    "bilstm_hidden_states = lstm_concatenate_layer([forward_hidden_states, reverse_hidden_states])\n",
    "\n",
    "# bilstm_layer = Bidirectional(LSTM(\n",
    "#     units=int(seq_size/adjacent_bp_pool_size),\n",
    "#     return_sequences=True,\n",
    "#     input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "#     name = 'bilstm_layer'))\n",
    "# bilstm_hidden_states = bilstm_layer(pooled_scores)\n",
    "print('bilstm_hidden_states', bilstm_hidden_states.get_shape())\n",
    "\n",
    "### attention layer ###\n",
    "# transpose hidden states\n",
    "# transpose_layer = Lambda(lambda x: K.permute_dimensions(x,(0,2,1)),\n",
    "#     name='transpose_layer')\n",
    "# transposed_hidden_states = transpose_layer(bilstm_hidden_states)\n",
    "# print('transposed_hidden_states', transposed_hidden_states.get_shape())\n",
    "\n",
    "# tanh layer\n",
    "attention_tanh_layer = Dense(attention_dim,\n",
    "    activation='tanh',\n",
    "    use_bias=False,\n",
    "    name = 'attention_tanh_layer')\n",
    "attention_tanh_layer_out = attention_tanh_layer(bilstm_hidden_states)\n",
    "print('attention_tanh_layer_out', attention_tanh_layer_out.get_shape())\n",
    "\n",
    "# rotate_layer = Lambda(lambda x: K.permute_dimensions(x,(0,2,1)),\n",
    "# #     name='rotate_layer'\n",
    "#     )\n",
    "# rotated_attention_tanh_layer_out = rotate_layer(attention_tanh_layer_out)\n",
    "\n",
    "# outer layer\n",
    "attention_outer_layer = Dense(attention_hops,\n",
    "    activation='relu',\n",
    "    use_bias=False,\n",
    "    name = 'attention_outer_layer')\n",
    "attention_outer_layer_out = attention_outer_layer(attention_tanh_layer_out)\n",
    "print('attention_outer_layer_out', attention_outer_layer_out.get_shape())\n",
    "\n",
    "# apply softmax\n",
    "softmax_layer = Softmax(axis=1, name='attention_softmax_layer')\n",
    "attention_softmax_layer_out = softmax_layer(attention_outer_layer_out)\n",
    "print('attention_softmax_layer_out', attention_softmax_layer_out.get_shape())\n",
    "\n",
    "# attend to hidden states\n",
    "attending_layer = Dot(axes=(1,1),\n",
    "    name='attending_layer')\n",
    "\n",
    "attended_states = attending_layer([attention_softmax_layer_out, bilstm_hidden_states])\n",
    "print('attended_states', attended_states.get_shape())\n",
    "\n",
    "# # fully connected layer\n",
    "dense_layer = Dense(num_dense_neurons, \n",
    "    activation='relu', \n",
    "    name = 'dense_layer'\n",
    "    )\n",
    "\n",
    "dense_output = dense_layer(attended_states)\n",
    "\n",
    "# drop out\n",
    "drop_out = Dropout(0.25,name='dense_dropout')(dense_output)\n",
    "\n",
    "# make prediction\n",
    "flattened = Flatten(name='flatten')(drop_out)\n",
    "predictions = Dense(num_classes,\n",
    "                    name='predictions',\n",
    "                    activation = 'softmax', \n",
    "                   )(flattened)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define and compile model\n",
    "attention_model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "attention_model.compile(loss=keras.losses.categorical_crossentropy,\n",
    "              optimizer=keras.optimizers.Adam(),\n",
    "              metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "keras.utils.plot_model(attention_model, to_file='attention_model.pdf')\n",
    "SVG(model_to_dot(attention_model, show_shapes=True, show_layer_names=True).create(prog='dot', format='svg'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "parallel_attention_model = keras.utils.multi_gpu_model(attention_model, gpus=2)\n",
    "parallel_attention_model.compile(loss=keras.losses.categorical_crossentropy,\n",
    "              optimizer=keras.optimizers.Adam(),\n",
    "              metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "parallel_attention_model.fit([x_train, x_rc_train], y_train,\n",
    "          batch_size=200,\n",
    "          epochs=50,\n",
    "          verbose=1,\n",
    "          validation_data=([x_test, x_rc_test], y_test))\n",
    "score = parallel_attention_model.evaluate([x_test, x_rc_test], y_test, verbose=0)\n",
    "print('Test loss:', score[0])\n",
    "print('Test accuracy:', score[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "probs = attention_model.predict([x_test, x_rc_test])\n",
    "\n",
    "sklearn.metrics.roc_auc_score([y[1] for y in y_test], probs[:,1], )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "attention_output, full_attention, attended_sequence = get_sequence_attention(attention_model,\n",
    "                      str(positive_seqRecords[2].seq),\n",
    "                      5)\n",
    "frame = pd.DataFrame(full_attention).T\n",
    "frame.columns = list(attended_sequence)\n",
    "fig, ax = plt.subplots(figsize=(30,1)) \n",
    "hm=sns.heatmap(frame, square=True, ax=ax, vmin=0, vmax=0.10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "attention_output, full_attention, attended_sequence = get_sequence_attention(attention_model,\n",
    "                      str(positive_seqRecords[9000].seq),\n",
    "                      5)\n",
    "frame = pd.DataFrame(full_attention).T\n",
    "frame.columns = list(attended_sequence)\n",
    "fig, ax = plt.subplots(figsize=(30,1)) \n",
    "hm=sns.heatmap(frame, square=True, ax=ax, vmin=0, vmax=0.10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "attention_output, full_attention, attended_sequence = get_sequence_attention(attention_model,\n",
    "                      str(positive_seqRecords[999].seq),\n",
    "                      5)\n",
    "frame = pd.DataFrame(full_attention).T\n",
    "frame.columns = list(attended_sequence)\n",
    "fig, ax = plt.subplots(figsize=(30,1)) \n",
    "hm=sns.heatmap(frame, square=True, ax=ax, vmin=0, vmax=0.10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_sequence_attention(model, sequence, adjacent_bp_pool_size):\n",
    "    get_attention = K.function([model.get_layer('input_fwd').input, \n",
    "                                model.get_layer('input_rev').input,\n",
    "                                K.learning_phase()\n",
    "                               ], \n",
    "                               [model.get_layer('attention_softmax_layer').output])\n",
    "    fwd_seq = sequence[:200]\n",
    "    rev_seq = str(Bio.Seq.Seq(fwd_seq).reverse_complement())\n",
    "    \n",
    "    fwd_seq_array = convert_sequences_to_array([fwd_seq])[0]\n",
    "    rev_seq_array = convert_sequences_to_array([rev_seq])[0]\n",
    "\n",
    "    layer_output = get_attention(([fwd_seq_array], [rev_seq_array], 0))[0]\n",
    "    reshaped_output = layer_output.reshape((layer_output.shape[1], layer_output.shape[2]))\n",
    "\n",
    "\n",
    "    full_attention = []\n",
    "    for x in reshaped_output:\n",
    "        for i in range(adjacent_bp_pool_size):\n",
    "            full_attention.append(x)\n",
    "    full_attention = np.array(full_attention)\n",
    "\n",
    "    crop_distance = int((len(fwd_seq) - full_attention.shape[0])/2)\n",
    "\n",
    "    attended_sequence = fwd_seq[crop_distance:-crop_distance]\n",
    "    return layer_output, full_attention, attended_sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def convert_arrays_to_sequences(sequence_arrays):\n",
    "    sequence_list = []\n",
    "    for arr in sequence_arrays:\n",
    "        current_seq = ''\n",
    "        for pos in arr:\n",
    "            if int(pos[0]) == 1:\n",
    "                current_seq += 'A'\n",
    "            elif int(pos[1]) == 1:\n",
    "                current_seq += 'C'\n",
    "            elif int(pos[2]) == 1:\n",
    "                current_seq += 'G'\n",
    "            elif int(pos[3]) == 1:\n",
    "                current_seq += 'T'\n",
    "            else:\n",
    "                current_seq += 'N'\n",
    "        sequence_list.append(current_seq)\n",
    "    return sequence_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train Models For each signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_attention_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        attention_hops,\n",
    "                        num_dense_neurons):\n",
    "    input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "    input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "    ### find motifs ###\n",
    "    convolution_layer = Conv1D(filters=num_motifs, \n",
    "        kernel_size=motif_size,\n",
    "        activation='relu',\n",
    "        input_shape=(total_seq_length,4),\n",
    "        name='convolution_layer',\n",
    "        padding = 'same'\n",
    "        )\n",
    "    forward_motif_scores = convolution_layer(input_fwd)\n",
    "    reverse_motif_scores = convolution_layer(input_rev)\n",
    "    print('forward_motif_scores', forward_motif_scores.get_shape())\n",
    "\n",
    "    ### crop motif scores to avoid parts of sequence where motif score is computed in only one direction ###\n",
    "    to_crop = int((total_seq_length - seq_size)/2)\n",
    "    crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "        name='crop_layer')\n",
    "    cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "    cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "    print('cropped_fwd_scores', cropped_fwd_scores.get_shape())\n",
    "\n",
    "    ### flip motif scores ###\n",
    "    flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "        output_shape=(seq_size, num_motifs),\n",
    "        name='flip_layer')\n",
    "    flipped_rev_scores = flip_layer(cropped_rev_scores)\n",
    "    print('flipped_rev_scores', flipped_rev_scores.get_shape())\n",
    "\n",
    "    ### concatenate motif scores ###\n",
    "    concatenate_layer = keras.layers.Concatenate(axis=2, name='concatenate_layer')\n",
    "    concatenated_motif_scores = concatenate_layer([cropped_fwd_scores, flipped_rev_scores])\n",
    "    print('concatenated_motif_scores', concatenated_motif_scores.get_shape())\n",
    "\n",
    "    ### pool across length of sequence ###\n",
    "    sequence_pooling_layer = MaxPool1D(pool_size=adjacent_bp_pool_size, \n",
    "        strides=adjacent_bp_pool_size,\n",
    "        name='sequence_pooling_layer')\n",
    "    pooled_scores = sequence_pooling_layer(concatenated_motif_scores)\n",
    "    print('pooled_scores', pooled_scores.get_shape())\n",
    "\n",
    "    ## bidirectional LSTM ###\n",
    "    forward_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "        return_sequences=True,\n",
    "        input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "        name = 'forward_lstm_layer'\n",
    "        )\n",
    "    forward_hidden_states = forward_lstm_layer(pooled_scores)\n",
    "    print('forward_hidden_states', forward_hidden_states.get_shape())\n",
    "\n",
    "    reverse_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "        return_sequences=True,\n",
    "        input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "        name = 'reverse_lstm_layer',\n",
    "        go_backwards=True,\n",
    "        )\n",
    "    reverse_hidden_states = reverse_lstm_layer(pooled_scores)\n",
    "    print('reverse_hidden_states', reverse_hidden_states.get_shape())\n",
    "    ### concatenate lstm hidden states ###\n",
    "    lstm_concatenate_layer = Concatenate(axis=2)\n",
    "    bilstm_hidden_states = lstm_concatenate_layer([forward_hidden_states, reverse_hidden_states])\n",
    "\n",
    "    print('bilstm_hidden_states', bilstm_hidden_states.get_shape())\n",
    "\n",
    "    # tanh layer\n",
    "    attention_tanh_layer = Dense(attention_dim,\n",
    "        activation='tanh',\n",
    "        use_bias=False,\n",
    "        name = 'attention_tanh_layer')\n",
    "    attention_tanh_layer_out = attention_tanh_layer(bilstm_hidden_states)\n",
    "    print('attention_tanh_layer_out', attention_tanh_layer_out.get_shape())\n",
    "\n",
    "    # outer layer\n",
    "    attention_outer_layer = Dense(attention_hops,\n",
    "        activation='relu',\n",
    "        use_bias=False,\n",
    "        name = 'attention_outer_layer')\n",
    "    attention_outer_layer_out = attention_outer_layer(attention_tanh_layer_out)\n",
    "    print('attention_outer_layer_out', attention_outer_layer_out.get_shape())\n",
    "\n",
    "    # apply softmax\n",
    "    softmax_layer = Softmax(axis=1, name='attention_softmax_layer')\n",
    "    attention_softmax_layer_out = softmax_layer(attention_outer_layer_out)\n",
    "    print('attention_softmax_layer_out', attention_softmax_layer_out.get_shape())\n",
    "\n",
    "    # attend to hidden states\n",
    "    attending_layer = Dot(axes=(1,1),\n",
    "        name='attending_layer')\n",
    "\n",
    "    attended_states = attending_layer([attention_softmax_layer_out, bilstm_hidden_states])\n",
    "    print('attended_states', attended_states.get_shape())\n",
    "\n",
    "    # # fully connected layer\n",
    "    dense_layer = Dense(num_dense_neurons, \n",
    "        activation='relu', \n",
    "        name = 'dense_layer'\n",
    "        )\n",
    "\n",
    "    dense_output = dense_layer(attended_states)\n",
    "\n",
    "    # drop out\n",
    "    drop_out = Dropout(0.25,name='dense_dropout')(dense_output)\n",
    "\n",
    "    # make prediction\n",
    "    flattened = Flatten(name='flatten')(drop_out)\n",
    "    predictions = Dense(num_classes,\n",
    "                        name='predictions',\n",
    "                        activation = 'softmax', \n",
    "                       )(flattened)\n",
    "    # define and compile model\n",
    "    model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "    model.compile(loss=keras.losses.categorical_crossentropy,\n",
    "                  optimizer=keras.optimizers.Adam(),\n",
    "                  metrics=['accuracy'])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_seq_length = len(fasta_seq[0])\n",
    "seq_size = 150\n",
    "num_classes = 2\n",
    "num_motifs = 50\n",
    "motif_size = 15\n",
    "adjacent_bp_pool_size = 5\n",
    "attention_dim = 200\n",
    "attention_hops = 1\n",
    "num_dense_neurons = 500 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### KLA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "positive_seqRecords = list(SeqIO.parse('./fasta_files/c57bl6_kla-1h_peaks.fa', 'fasta'))\n",
    "negative_seqRecords = list(SeqIO.parse('./background_files/c57bl6_kla-1h_peaks_background.fasta', 'fasta'))\n",
    "\n",
    "fasta_seq = [str(x.seq[:200]) for x in positive_seqRecords] + [str(x[:200].seq) for x in negative_seqRecords]\n",
    "\n",
    "fasta_rc_seq = [str(x[:200].reverse_complement().seq) for x in positive_seqRecords] + \\\n",
    "    [str(x[:200].reverse_complement().seq) for x in negative_seqRecords]\n",
    "\n",
    "sequence_arrays = convert_sequences_to_array(fasta_seq)\n",
    "sequence_arrays = np.array(sequence_arrays)\n",
    "\n",
    "sequence_rc_arrays = convert_sequences_to_array(fasta_rc_seq)\n",
    "sequence_rc_arrays = np.array(sequence_rc_arrays)\n",
    "\n",
    "\n",
    "labels = [1 for x in positive_seqRecords] + [0 for x in negative_seqRecords]\n",
    "labels = np.array(labels)\n",
    "\n",
    "x_train, x_test, x_rc_train, x_rc_test, y_train, y_test = model_selection.train_test_split(sequence_arrays, sequence_rc_arrays, labels, test_size=0.2)\n",
    "\n",
    "num_classes = 2\n",
    "y_train = keras.utils.to_categorical(y_train, num_classes)\n",
    "y_test = keras.utils.to_categorical(y_test, num_classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "kla_model = get_attention_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        attention_hops,\n",
    "                        num_dense_neurons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "parallel_attention_model = keras.utils.multi_gpu_model(kla_model, gpus=2)\n",
    "parallel_attention_model.compile(loss=keras.losses.binary_crossentropy,\n",
    "              optimizer=keras.optimizers.Adam(),\n",
    "              metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "parallel_attention_model.fit([x_train, x_rc_train], y_train,\n",
    "          batch_size=200,\n",
    "          epochs=50,\n",
    "          verbose=1,\n",
    "          validation_data=([x_test, x_rc_test], y_test))\n",
    "score = parallel_attention_model.evaluate([x_test, x_rc_test], y_test, verbose=0)\n",
    "print('Test loss:', score[0])\n",
    "print('Test accuracy:', score[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# serialize model to JSON\n",
    "model_json = kla_model.to_json()\n",
    "with open(\"kla_model.json\", \"w\") as json_file:\n",
    "    json_file.write(model_json)\n",
    "# serialize weights to HDF5\n",
    "kla_model.save_weights(\"kla_model.h5\")\n",
    "print(\"Saved model to disk\")\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Restore Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load json and create model\n",
    "json_file = open('kla_model.json', 'r')\n",
    "loaded_model_json = json_file.read()\n",
    "json_file.close()\n",
    "loaded_model = keras.models.model_from_json(loaded_model_json)\n",
    "# load weights into new model\n",
    "loaded_model.load_weights(\"kla_model.h5\")\n",
    "print(\"Loaded model from disk\")\n",
    "loaded_model.compile(loss=keras.losses.categorical_crossentropy,\n",
    "              optimizer=keras.optimizers.Adam(),\n",
    "              metrics=['accuracy'])\n",
    "kla_model = loaded_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "arg1_enhancer = 'ggtagccgacgagagaccagctcatcttcaataaggaagtcagagagcagaaggctttgtcagcagggcaagactatactttgttaggaagtgaggcattgttcagacttccttatgctttcttatgaacaggctgtattagccaacagtcctgtc'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "arg1_enhancer = 'AAAGTGGCACAACTCACGTACAGACAGGACTGTTGGCTAATACAGCCTGTTCATAAGAAAGCATAAGGAAGTCTGAACAATGCCTCACTTCCTAACAAAGTATAGTCTTGCCCTGCTGACAAAGCCTTCTGCTCTCTGACTTCCTTATTGAAGATGAGCTGGTCTCTCGTCGGCTACCACCCTCCGTGACCTTATGCAGA'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tnf_enhancer = 'CTAAGCTGTGTCACGGGAGCTGGCAGCACGCTGGCGGATATGCCTTGCCATGGGCCAATTTTGGTTTCAATCTCAGTTTTAGAGGTTGTGTGAAATTCAGTTTCTCTCTTGGGGAGGCCAACAGCTGTCTGGGACTTTCCCCGGGGGGGAGGGCTGATGACTAGGAGTCTTGTGCATCGTCTATAACCACTCTCAGGAAG'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "index = 0\n",
    "\n",
    "\n",
    "attention_output, full_attention, attended_sequence = get_sequence_attention(kla_model,\n",
    "                      str(positive_seqRecords[index].seq),\n",
    "                      5)\n",
    "frame = pd.DataFrame(full_attention).T\n",
    "frame.columns = list(attended_sequence)\n",
    "fig, ax = plt.subplots(figsize=(30,1)) \n",
    "hm=sns.heatmap(frame, square=True, ax=ax, vmin=0, vmax=0.20)\n",
    "plt.ylabel('KLA')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multilabel Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_attention_multilabel_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        attention_hops,\n",
    "                        num_dense_neurons,\n",
    "                        dropout_rate=0.25):\n",
    "    input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "    input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "    ### find motifs ###\n",
    "    convolution_layer = Conv1D(filters=num_motifs, \n",
    "        kernel_size=motif_size,\n",
    "        activation='relu',\n",
    "        input_shape=(total_seq_length,4),\n",
    "        name='convolution_layer',\n",
    "        padding = 'same'\n",
    "        )\n",
    "    forward_motif_scores = convolution_layer(input_fwd)\n",
    "    reverse_motif_scores = convolution_layer(input_rev)\n",
    "    print('forward_motif_scores', forward_motif_scores.get_shape())\n",
    "\n",
    "    ### crop motif scores to avoid parts of sequence where motif score is computed in only one direction ###\n",
    "    to_crop = int((total_seq_length - seq_size)/2)\n",
    "    crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "        name='crop_layer')\n",
    "    cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "    cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "    print('cropped_fwd_scores', cropped_fwd_scores.get_shape())\n",
    "\n",
    "    ### flip motif scores ###\n",
    "    flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "        output_shape=(seq_size, num_motifs),\n",
    "        name='flip_layer')\n",
    "    flipped_rev_scores = flip_layer(cropped_rev_scores)\n",
    "    print('flipped_rev_scores', flipped_rev_scores.get_shape())\n",
    "\n",
    "    ### concatenate motif scores ###\n",
    "    concatenate_layer = keras.layers.Concatenate(axis=2, name='concatenate_layer')\n",
    "    concatenated_motif_scores = concatenate_layer([cropped_fwd_scores, flipped_rev_scores])\n",
    "    print('concatenated_motif_scores', concatenated_motif_scores.get_shape())\n",
    "\n",
    "    ### pool across length of sequence ###\n",
    "    sequence_pooling_layer = MaxPool1D(pool_size=adjacent_bp_pool_size, \n",
    "        strides=adjacent_bp_pool_size,\n",
    "        name='sequence_pooling_layer')\n",
    "    pooled_scores = sequence_pooling_layer(concatenated_motif_scores)\n",
    "    print('pooled_scores', pooled_scores.get_shape())\n",
    "\n",
    "    ## bidirectional LSTM ###\n",
    "    forward_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "        return_sequences=True,\n",
    "        input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "        name = 'forward_lstm_layer'\n",
    "        )\n",
    "    forward_hidden_states = forward_lstm_layer(pooled_scores)\n",
    "    print('forward_hidden_states', forward_hidden_states.get_shape())\n",
    "\n",
    "    reverse_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "        return_sequences=True,\n",
    "        input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "        name = 'reverse_lstm_layer',\n",
    "        go_backwards=True,\n",
    "        )\n",
    "    reverse_hidden_states = reverse_lstm_layer(pooled_scores)\n",
    "    print('reverse_hidden_states', reverse_hidden_states.get_shape())\n",
    "    ### concatenate lstm hidden states ###\n",
    "    lstm_concatenate_layer = Concatenate(axis=2)\n",
    "    bilstm_hidden_states = lstm_concatenate_layer([forward_hidden_states, reverse_hidden_states])\n",
    "\n",
    "    print('bilstm_hidden_states', bilstm_hidden_states.get_shape())\n",
    "\n",
    "    # tanh layer\n",
    "    attention_tanh_layer = Dense(attention_dim,\n",
    "        activation='tanh',\n",
    "        use_bias=False,\n",
    "        name = 'attention_tanh_layer')\n",
    "    attention_tanh_layer_out = attention_tanh_layer(bilstm_hidden_states)\n",
    "    print('attention_tanh_layer_out', attention_tanh_layer_out.get_shape())\n",
    "\n",
    "    # outer layer\n",
    "    attention_outer_layer = Dense(attention_hops,\n",
    "        activation='linear',\n",
    "        use_bias=False,\n",
    "        name = 'attention_outer_layer')\n",
    "    attention_outer_layer_out = attention_outer_layer(attention_tanh_layer_out)\n",
    "    print('attention_outer_layer_out', attention_outer_layer_out.get_shape())\n",
    "\n",
    "    # apply softmax\n",
    "    softmax_layer = Softmax(axis=1, name='attention_softmax_layer')\n",
    "    attention_softmax_layer_out = softmax_layer(attention_outer_layer_out)\n",
    "    print('attention_softmax_layer_out', attention_softmax_layer_out.get_shape())\n",
    "\n",
    "    # attend to hidden states\n",
    "    attending_layer = Dot(axes=(1,1),\n",
    "        name='attending_layer')\n",
    "\n",
    "    attended_states = attending_layer([attention_softmax_layer_out, bilstm_hidden_states])\n",
    "    print('attended_states', attended_states.get_shape())\n",
    "\n",
    "    # # fully connected layer\n",
    "    dense_layer = Dense(num_dense_neurons, \n",
    "        activation='relu', \n",
    "        name = 'dense_layer'\n",
    "        )\n",
    "\n",
    "    dense_output = dense_layer(attended_states)\n",
    "    print('dense_output', dense_output.shape)\n",
    "    \n",
    "    # drop out\n",
    "    drop_out = Dropout(dropout_rate,name='dense_dropout')(dense_output)\n",
    "    print('drop_out', drop_out.shape)\n",
    "    \n",
    "    # make prediction\n",
    "    flattened = Flatten(name='flatten')(drop_out)\n",
    "    print('flattened', flattened.shape)\n",
    "    \n",
    "    predictions = Dense(num_classes,\n",
    "                        name='predictions',\n",
    "                        activation = 'sigmoid', \n",
    "                       )(flattened)\n",
    "    print('predictions', predictions.shape)\n",
    "    \n",
    "    # define and compile model\n",
    "    model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "    model.compile(loss=keras.losses.binary_crossentropy,\n",
    "                  optimizer=keras.optimizers.Adam(),\n",
    "                  metrics=['binary_accuracy'])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_normed_multilabel_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        attention_hops,\n",
    "                        num_dense_neurons,\n",
    "                        dropout_rate=0.25):\n",
    "    input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "    input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "    ### find motifs ###\n",
    "    convolution_layer = Conv1D(filters=num_motifs, \n",
    "        kernel_size=motif_size,\n",
    "        activation='relu',\n",
    "        input_shape=(total_seq_length,4),\n",
    "        name='convolution_layer',\n",
    "        padding = 'same'\n",
    "        )\n",
    "    forward_motif_scores = convolution_layer(input_fwd)\n",
    "    reverse_motif_scores = convolution_layer(input_rev)\n",
    "    print('forward_motif_scores', forward_motif_scores.get_shape())\n",
    "\n",
    "    ### crop motif scores to avoid parts of sequence where motif score is computed in only one direction ###\n",
    "    to_crop = int((total_seq_length - seq_size)/2)\n",
    "    crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "        name='crop_layer')\n",
    "    cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "    cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "    print('cropped_fwd_scores', cropped_fwd_scores.get_shape())\n",
    "\n",
    "    ### flip motif scores ###\n",
    "    flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "        output_shape=(seq_size, num_motifs),\n",
    "        name='flip_layer')\n",
    "    flipped_rev_scores = flip_layer(cropped_rev_scores)\n",
    "    print('flipped_rev_scores', flipped_rev_scores.get_shape())\n",
    "\n",
    "    ### concatenate motif scores ###\n",
    "    concatenate_layer = keras.layers.Concatenate(axis=2, name='concatenate_layer')\n",
    "    concatenated_motif_scores = concatenate_layer([cropped_fwd_scores, flipped_rev_scores])\n",
    "    print('concatenated_motif_scores', concatenated_motif_scores.get_shape())\n",
    "\n",
    "    ### pool across length of sequence ###\n",
    "    sequence_pooling_layer = MaxPool1D(pool_size=adjacent_bp_pool_size, \n",
    "        strides=adjacent_bp_pool_size,\n",
    "        name='sequence_pooling_layer')\n",
    "    pooled_scores = sequence_pooling_layer(concatenated_motif_scores)\n",
    "    print('pooled_scores', pooled_scores.get_shape())\n",
    "    \n",
    "    ### normalize motif scores ###\n",
    "    motif_score_norm_layer = BatchNormalization(name='motif_score_norm_layer')\n",
    "    normed_pooled_scores = motif_score_norm_layer(pooled_scores)\n",
    "    print('normed_pooled_scores', normed_pooled_scores.shape)\n",
    "    \n",
    "    ### bidirectional LSTM ###\n",
    "    forward_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "        return_sequences=True,\n",
    "        input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "        name = 'forward_lstm_layer'\n",
    "        )\n",
    "    forward_hidden_states = forward_lstm_layer(normed_pooled_scores)\n",
    "    print('forward_hidden_states', forward_hidden_states.get_shape())\n",
    "\n",
    "    reverse_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "        return_sequences=True,\n",
    "        input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "        name = 'reverse_lstm_layer',\n",
    "        go_backwards=True,\n",
    "        )\n",
    "    reverse_hidden_states = reverse_lstm_layer(normed_pooled_scores)\n",
    "    print('reverse_hidden_states', reverse_hidden_states.get_shape())\n",
    "    \n",
    "    ### concatenate lstm hidden states ###\n",
    "    lstm_concatenate_layer = Concatenate(axis=2)\n",
    "    bilstm_hidden_states = lstm_concatenate_layer([forward_hidden_states, reverse_hidden_states])\n",
    "\n",
    "    print('bilstm_hidden_states', bilstm_hidden_states.get_shape())\n",
    "    \n",
    "    ### normalize lstm states ###\n",
    "    lstm_norm_layer = BatchNormalization(name='lstm_norm_layer')\n",
    "    normed_bilistm_hidden_states = lstm_norm_layer(bilstm_hidden_states)\n",
    "    \n",
    "    ### attention tanh layer ###\n",
    "    attention_tanh_layer = Dense(attention_dim,\n",
    "        activation='tanh',\n",
    "        use_bias=False,\n",
    "        name = 'attention_tanh_layer')\n",
    "    attention_tanh_layer_out = attention_tanh_layer(normed_bilistm_hidden_states)\n",
    "    print('attention_tanh_layer_out', attention_tanh_layer_out.get_shape())\n",
    "\n",
    "    ### outer layer ###\n",
    "    attention_outer_layer = Dense(attention_hops,\n",
    "        activation='linear',\n",
    "        use_bias=False,\n",
    "        name = 'attention_outer_layer')\n",
    "    attention_outer_layer_out = attention_outer_layer(attention_tanh_layer_out)\n",
    "    print('attention_outer_layer_out', attention_outer_layer_out.get_shape())\n",
    "\n",
    "    ### apply softmax ###\n",
    "    softmax_layer = Softmax(axis=1, name='attention_softmax_layer')\n",
    "    attention_softmax_layer_out = softmax_layer(attention_outer_layer_out)\n",
    "    print('attention_softmax_layer_out', attention_softmax_layer_out.get_shape())\n",
    "\n",
    "    ### attend to hidden states ###\n",
    "    attending_layer = Dot(axes=(1,1),\n",
    "        name='attending_layer')\n",
    "\n",
    "    attended_states = attending_layer([attention_softmax_layer_out, normed_bilistm_hidden_states])\n",
    "    print('attended_states', attended_states.get_shape())\n",
    "    \n",
    "    ### normalize attended states ###\n",
    "    attention_norm_layer = BatchNormalization(name='attention_norm_layer')\n",
    "    normed_attended_states = attention_norm_layer(attended_states)\n",
    "    print('normed_attended_states', normed_attended_states.shape)\n",
    "    \n",
    "    ### fully connected layer ###\n",
    "    dense_layer = Dense(num_dense_neurons, \n",
    "        activation='relu', \n",
    "        name = 'dense_layer'\n",
    "        )\n",
    "\n",
    "    dense_output = dense_layer(normed_attended_states)\n",
    "    print('dense_output', dense_output.shape)\n",
    "    \n",
    "    # drop out\n",
    "    drop_out = Dropout(dropout_rate,name='dense_dropout')(dense_output)\n",
    "    print('drop_out', drop_out.shape)\n",
    "    \n",
    "    # make prediction\n",
    "    flattened = Flatten(name='flatten')(drop_out)\n",
    "    print('flattened', flattened.shape)\n",
    "    \n",
    "    predictions = Dense(num_classes,\n",
    "                        name='predictions',\n",
    "                        activation = 'sigmoid', \n",
    "                       )(flattened)\n",
    "    print('predictions', predictions.shape)\n",
    "    \n",
    "    # define and compile model\n",
    "    model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "    model.compile(loss=keras.losses.binary_crossentropy,\n",
    "                  optimizer=keras.optimizers.Adam(),\n",
    "                  metrics=['binary_accuracy'])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "negative_seqRecords = list(SeqIO.parse('./background.fasta', 'fasta'))\n",
    "negative_seq = [str(x.seq[:200]) for x in negative_seqRecords]\n",
    "negative_rc_seq = [str(x.seq[:200]) for x in negative_seqRecords]\n",
    "\n",
    "negative_sequence_arrays = convert_sequences_to_array(negative_seq)\n",
    "negative_sequence_arrays = np.array(negative_sequence_arrays)\n",
    "\n",
    "negative_sequence_rc_arrays = convert_sequences_to_array(negative_rc_seq)\n",
    "negative_sequence_rc_arrays = np.array(negative_sequence_rc_arrays)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "positive_seqRecords = list(SeqIO.parse('./merged_atac_peaks_filtered_resized.fasta', 'fasta'))\n",
    "positive_seqRecords = [x for x in positive_seqRecords]\n",
    "fasta_seq = [str(x.seq[:200]) for x in positive_seqRecords]\n",
    "\n",
    "fasta_rc_seq = [str(x[:200].reverse_complement().seq) for x in positive_seqRecords ]\n",
    "\n",
    "seq_ids = [x.name for x in positive_seqRecords]\n",
    "\n",
    "sequence_arrays = convert_sequences_to_array(fasta_seq)\n",
    "sequence_arrays = np.array(sequence_arrays)\n",
    "\n",
    "sequence_rc_arrays = convert_sequences_to_array(fasta_rc_seq)\n",
    "sequence_rc_arrays = np.array(sequence_rc_arrays)\n",
    "\n",
    "index_seqArray_dict = dict(zip(seq_ids, zip(sequence_arrays, sequence_rc_arrays)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "summary_frame = pd.read_csv('./group_atac_summary.tsv' , sep='\\t', low_memory=False)\n",
    "summary_frame = summary_frame.fillna('0')\n",
    "for col in summary_frame.columns[5:]:\n",
    "    floatValues = []\n",
    "    for val in summary_frame[col].values.astype(str):\n",
    "        if ',' in val:\n",
    "            maxVal = np.mean([float(x) for x in val.split(',')])\n",
    "            floatValues.append(maxVal)\n",
    "        else:\n",
    "            floatValues.append(float(val))\n",
    "    summary_frame[col] = floatValues\n",
    "# summary_frame.index = summary_frame['ID'].values\n",
    "summary_frame.index = summary_frame['chr'] + ':' + (summary_frame['start'] - 1).astype(str) + '-' + summary_frame['end'].astype(str)\n",
    "\n",
    "# remove peaks in unknown/random chromosomes\n",
    "summary_frame = summary_frame[~summary_frame['chr'].str.contains('random')]\n",
    "summary_frame = summary_frame[~summary_frame['chr'].str.contains('Un')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "c57bl6_ifng-24h    33333\n",
       "c57bl6_il13-24h    23066\n",
       "c57bl6_il1b-24h    13650\n",
       "c57bl6_il23-24h    33482\n",
       "c57bl6_il4-24h     33518\n",
       "c57bl6_il5-24h      2091\n",
       "c57bl6_il6-24h     20308\n",
       "c57bl6_kla-1h      21776\n",
       "c57bl6_tgfb-24h    29540\n",
       "c57bl6_tnfa-24h    31600\n",
       "c57bl6_veh-24h     11500\n",
       "c57bl6_veh         31550\n",
       "dtype: int64"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(summary_frame[[x for x in summary_frame.columns if 'c57' in x]]>0).sum(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "c57bl6_ifng-24h    0.608588\n",
       "c57bl6_il13-24h    0.421135\n",
       "c57bl6_il1b-24h    0.249219\n",
       "c57bl6_il23-24h    0.611309\n",
       "c57bl6_il4-24h     0.611966\n",
       "c57bl6_il5-24h     0.038177\n",
       "c57bl6_il6-24h     0.370780\n",
       "c57bl6_kla-1h      0.397583\n",
       "c57bl6_tgfb-24h    0.539337\n",
       "c57bl6_tnfa-24h    0.576948\n",
       "c57bl6_veh-24h     0.209965\n",
       "c57bl6_veh         0.576035\n",
       "dtype: float64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(summary_frame[[x for x in summary_frame.columns if 'c57' in x]]>0).mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jenhan/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "peak_count_threshold = 20000\n",
    "peak_counts = (summary_frame[[x for x in summary_frame.columns if 'c57' in x]]>0).sum(axis=0)\n",
    "to_filter_conditions = peak_counts[peak_counts < peak_count_threshold].index.values\n",
    "label_frame = summary_frame[[x for x in summary_frame.columns if 'c57' in x]]\n",
    "for cond in to_filter_conditions:\n",
    "    label_frame.drop(cond, axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "label_frame = summary_frame[['c57bl6_il4-24h', 'c57bl6_kla-1h', 'c57bl6_veh']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "quantiled_label_frame = quantile_normalize_df(label_frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(54771, 3)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "label_frame.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f85589a5208>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAFcCAYAAADh1zYWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8HWd99/3PzFm177K8yvvEjmNnj0MSEkhYbqDc7C19\nWppSmkIDD9CHQoCGu5RSbrjL2qQt0BJSCAFuoBAgG04cZ3ES7DhxbMcZb/Im2dq3o6Ozzjx/HB1H\nXiQdHWl0dOTv+/XKK9LMnJmfxjr6neua6/pdhuu6iIiISPEwCx2AiIiITI6St4iISJFR8hYRESky\nSt4iIiJFRslbRESkyCh5i4iIFBm/lye3LGsd8Cvg67Zt33HGvpuAfwLSwP22bX/By1hERETmCs9a\n3pZllQH/AjwyxiHfAt4JXAO83rKstV7FIiIiMpd42W0eB94EtJ25w7Ks5UCPbdvHbNt2gPuBGz2M\nRUREZM7wLHnbtp2ybXt4jN1NQOeo7zuA+V7FIiIiMpd4+sx7EoyJDkil0q7f75uJWDz1jd/9csx9\nH3vd22YwEhERKQLnzI+FSt5tZFrfWQs5R/f6aL29UU8DmknxWPKc2zs7B2c4ktmtoaFC9yQHuk+5\n073Kne5Vbry+Tw0NFefcXpCpYrZtHwYqLctaalmWH3gL8HAhYhERESk2nrW8Lcu6DPgqsBRIWpb1\nLuA+oMW27f8GPgTcO3L4T2zb3udVLCIiInOJZ8nbtu3ngBvG2f84cLVX1xcREZmrVGFNRESkyCh5\ni4iIFBklbxERkSKj5C0iIlJkZkuRFhERkbN8+MO38Dd/80mWL195atvQUIQvfekL9Pb24Dhpqqqq\n+exnP09FxbnnRM9FSt4iIueJx15ozfnYivIwg5HYuMfccPHCqYaUl5/85EesXXshf/zH7wPg+9//\nDx5++AHe+c73FCSeQlDyFhERz9x//6959tmtDA0N0dnZwXve88csWrSYb3/7Tvx+P42N8/jUp/4O\nwzD44hf/ns7ODoaHh3n/+2/hmmuuO3WeoaEIH/vYrXz6058jEhkklUqd2nfzzR849fU999zNY489\ngmGYfPCDH+bSSy/npz+9l0ceydQBu+666/mTP7mZL37x7/H7AwwM9PEP//C/+cpXvkhbWyupVIoP\nfOCDXHbZFTN3k/Kg5C0iIp5qaTnE9753D5FIhJtvfi81NTV885v/RmVlFf/6r99k8+ZNXHHFVVx5\n5Ub+x/94C62tx7n99ttOJW/XhX/8x7/n/e+/heXLV/COd7yHj3/8wzzzzFNceeXV3Hjj61m1ajXH\njh3lscce4dvf/j5tba388Iffp6lpPg888Gu++93/AuCWW/6M17zmJgAqKyv51Kc+y4MP/pa6uno+\n/enP0dfXx0c/+kHuvvvHBbpbuVHyFhERT1188aX4/X6qq6spKyvj6NEjfOYzfwtALBajqqqaiopK\n9u7dw333/QLDMBkY6D/1+rvu+g7z5s3j6quvAWDRosXce+/P2bFjO88++zQf+9iH+NCH/l9KSkpZ\nu3YdpmmyaNFibrvtdrZseZQLL7wIvz+T7i66aAMHDmQKeq5deyEAu3e/yM6dz/Piiy8AEI/HSSaT\nBAKBGbtHk6XkLSIinnIc99TXhmFSV1fPHXd857RjHnjgNwwMDHDnnf/BwMAAH/jAn57aV1FRybZt\nz9Lf30dVVTXxeIxQKMyVV27kyis3cu21r+Z73/sO73zne0671sgVcd1XtiWTSQwjM9HK7w+c+v/7\n3vd+Xve6N07zT+4dTRUTERFP7dnzIul0mr6+PqLRIUzTpKXlEAA/+9mPOXBgP319fcyfvwDTNNmy\n5VGSyVdWX3z3u/+IP/7j9/GNb/wzAB/72K1s2/bsqf2dnR0sWLAQy1rDrl07SaVS9PR08+lPf4LV\nqy12795FKpUilUrx0kt7WL3aOi2+tWvX8eSTWwDo7e3h29++0+tbMmVqeYuIiKeamhZw++230dp6\njFtu+Wvmz1/IP/3T5wkEAtTXN/DWt76DsrIybrvtb3jppd28+c1vpbGxkbvu+u6pc7z5zW/l0Uc3\n8eSTW/jMZ/4XX/val/n+9/8Dn89HeXkFn/jEbdTW1vGGN7yJD3/4FlzX5a/+6lbmz1/AW9/6dj7y\nkVtwHJc/+IP/SVPT/NPie+1rb2LHjm188IPvJ51O8/733zLTt2jSjNHdCbNZZ+dgcQQ6gXteeGTM\n9bzfv7F4umxmgtYTzo3uU+50r3I3Xffq/vt/zaFDB/nwhz82DVHNPjOwnrdxru3qNhcRESky6jYX\nERHPvOlNf1DoEOYktbxFRESKjJK3iIhIkVHyFhERKTJK3iIiIkVGyVtERIrOY489AmSmom3Zspkd\nO7bzd3/3yQJHNbFDhw7w4Q9PfR65RpuLiJwnnmx9JudjKwbCDA6OvyTotQs3TjWkvJw40camTQ9x\nww03nhrNvmPH9oLEUihK3iIi4omhoQif+cwnSSTiXHbZFTz00P383/97H+961x/wX//1E0pLS7nj\njm+wfPkKrr/+NXz+83/H8PAwsViMj3/8b1m7dh1/+Idv461vfTtbtz5JIpHgm9/8V772tS+zd+8e\n7rrruziOQ3V1NcuWrTh13S1bHuXHP/4hPp8fy1rDRz7y8dPi+uIX/566unr27XuZ9vaTfO5z/4hl\nXTDh0qHXXPNqXnhhB319fbS0HOKWWz7E448/wr59+/nc5/6RCy9cx89//lM2bXoQwzC57robeO97\n/4SOjnZuv/02AoEAK1eunpZ7q25zERHxxIMP3s+qVav5t3/7T5YuXcZ4FT27u7t5y1vexr/8y7f5\n4Ac/zD333A1AOp2muXkZd975XRYsWMD27dt473v/lIsvvpQ///O/POs80WiUu+/+T775zX/njju+\nQ0dH+6nVwkZLJpN87Wt38O53/xEPPvhb2tpaeeCBX3Pnnd/lzju/y6OP/o7W1uNAZunQL37x/wBw\n7NhRvvzlr/Gnf3ozP/zh97nzzjv50z+9mU2bHqKtrZXHHnuEf/3X/+TOO7/Lli2PcvLkSX72sx9z\n442v5447vkN9ff103Fq1vEVExBtHjrRwySWXAZz6/1hqa+u4++7/4N57f0AymSQcDp/at2HDJQA0\nNMxjaChCeXn5mOdpaTlEe/tJ/uZvPgxkWv8nT55k/frTjxt9zpde2sP+/faES4cCXHDBWgzDoK6u\nnhUrVuHz+aipqWNoaCd79+7h+PFjfOQjfwVANDrEyZNtHD7ccmoN8UsuuZxnntk6/o3LgZK3iIh4\nwnXBMDKluX2+V9JNdhtAKpUC4Kc//RH19Y3cfvsXePnll7jjjm+cOsbn84065/jLXAQCma7yr33t\njnGPO/ucEy8deubrzjyH3x/g6quv4ZOf/Oxp17rnnrtPnct1nXHjypW6zUVExBPNzc289NJuALZv\nf2UJz9LSMrq7u0in0+zZswuA/v4+Fi5cBMCWLZtPJfVzMU2TdDp9zn1Llizl8OEWent7APjP//w2\nnZ0dE8aay9KhE7GsNezY8RyxWAzXdfnGN/6ZeDzGkiXNvPzyS8D0DaxTy1tERDzxhje8mc985hPc\neutfsn79xae2v/Od7+FTn/o4S5Y0s2zZcgDe+MY384//+L/YvHkT73zne9i06WF++9v7znne5uZl\n2PbLfOtbX6Ws7PQu9HA4zEc/+v/xiU98lGAwwKpVFvX1DRPGmsvSoRNpamriPe95L7fe+peYpsmr\nX30DoVCYd7/7vdx++208/vhmVqxYNalzjkVLgs4wLQmaOy3fmBvdp9zpXuVuuu9VNBrlfe/7Q372\ns19P2zlnAy0JKiIiIjlR8hYREc+VlpbOuVZ3ISl5i4iIFBklbxERkSKj5C0iIlJklLxFRESKjJK3\niIhIkVHyFhERKTJK3iIiIkVGyVtERKTIKHmLiIgUGSVvERGRIqPkLSIiUmSUvEVERIqMkreIiEiR\nUfIWEREpMkreIiIiRUbJW0REpMgoeYuIiBQZJW8REZEio+QtIiJSZJS8RUREioySt4iISJHxe3ly\ny7K+DmwEXOCjtm1vG7XvVuBPgDSw3bbtj3kZi4iIyFzhWcvbsqzrgVW2bV8N/AXwrVH7KoG/Ba6z\nbftaYK1lWRu9ikVERGQu8bLb/EbglwC2be8FakaSNkBi5L9yy7L8QCnQ42EsIiIic4aXybsJ6Bz1\nfefINmzbjgGfBw4BR4Bnbdve52EsIiIic4anz7zPYGS/GGmBfwZYDQwAj1qWtcG27Z1jvbimphS/\n3+d9lDMgFA6cc3tDQ8UMRzL76Z7kRvcpd7pXudO9yk0h7pOXybuNkZb2iAXAiZGv1wCHbNvuArAs\n6wngMmDM5N3bG/UozJkXjyXPub2zc3CGI5ndGhoqdE9yoPuUO92r3Ole5cbr+zTWBwMvu80fBt4F\nYFnWpUCbbdvZn/AwsMayrJKR7y8H9nsYi4iIyJzhWcvbtu2tlmU9Z1nWVsABbrUs62ag37bt/7Ys\n6/8Amy3LSgFbbdt+wqtYRERE5hJPn3nbtn3bGZt2jtr3beDbXl5fRERkLlKFNRERkSKj5C0iIlJk\nlLxFRESKjJK3iIhIkVHyFhERKTJK3iIiIkVGyVtERKTIKHmLiIgUGSVvERGRIqPkLSIiUmSUvEVE\nRIqMkreIiEiRUfIWEREpMkreIiIiRUbJW0REpMgoeYuIiBQZJW8REZEio+QtIiJSZJS8RUREioyS\n9wxJppO4rlvoMEREZA7wFzqA80HaSfP3z3yFimA5DeZyDAKFDklERIqYkvcM6Bzuoi/eT1+8nxNG\nO0uDF1Llqyt0WCIiUqTUbT4D2obaAWgKLiHtOhyI76QzebzAUYmISLFS8p4BbZGTABx5sZGFySsx\n8XEydbTAUYmISLFS8p4BJ4YyydsZLufo3ipKjHISbgzHdQocmYiIFCMl7xlwfPAkbiqAmQoTjxmk\nomUAJNxYgSMTEZFipOTtsUQ6SXesGydazk2XLyYUcol0lwIQd6MFjk5ERIqRkrfH2qMduLi4w+Vc\nsKSGlReAG8u0vOPOcIGjExGRYqTk7bETIyPNneEKls6vYF4TlPjDAAzGlLxFRGTylLw9lh1pXk4N\n1eUhDAOa6jLJe1gtbxERyYOSt8eO9rcB0Fy94NS20rAfNxkgaSh5i4jI5Cl5e6w1chI3EWLFvPpT\n28JhFzdeiuMbVr1zERGZNCVvDw2nYkTSAzjD5SybX3lqeyAIbrwUDFfTxUREZNKUvD10cmSwmjtc\nTnNTxanthgG+dHa6mLrORURkcpS8PZQdrFZKDRWlwdP2BdwSAIZTSt4iIjI5St4eOtTTCsDCivln\n7QsZmeQ9lFTyFhGRyVHy9tCRkZHmVv2is/aV+DLJW4VaRERkspS8PdSV6MSJh1k5/+y1u0vCfty0\njwRK3iIiMjlK3h6JJIdIMow7XEFzU+VZ+0vCBm6slJQZ1XQxERGZFCVvj/TF+gEIOuWUhv1n7c/O\n9cZ0SLqJmQ5PRESKmJK3R/rjAwCU+MrOud8fABKaLiYiIpOn5O2RrmgmeZf7z528Afzp7KA1LQ0q\nIiK5U/L2SOdgLwBVobOfd2cFjEzLO6q53iIiMglK3h7pGc60vGtKxk7eJeZIoZa0kreIiOROydsj\nffFBAOrLqsY8piQYxHVMEnrmLSIik6Dk7ZFIMpO8myqrxzymJAxuvISkocVJREQkd0reHommo7ip\nALXlYw9YC5WAmwzimklc15nB6EREpJgpeXsk4UZxk0GqK0JjHhMOu7jJzP4UyZkKTUREipyStwfS\nTpq0GYdkiIrSwJjH+XxgOpn9SVfJW0REcqPk7YHBZAQAv1uCaRjjHusjs1Ro0lGVNRERyY2Stwey\n1dXCI/O4xxMwMi3vWFLJW0REcnN20e1pZFnW14GNgAt81LbtbaP2LQbuBYLADtu2P+hlLDOpcyhT\n17x0nOpqWUFfkAQQS6rbXEREcuNZy9uyrOuBVbZtXw38BfCtMw75KvBV27avBNKWZS3xKpaZ1jGQ\nqa5WGSyf8NiwP9NtHk8reYuISG687Da/EfglgG3be4Eay7IqASzLMoHrgPtG9t9q2/ZRD2OZUdm6\n5lXhsaurZZUERwas6Zm3iIjkyMtu8ybguVHfd45sGwAagEHg65ZlXQo8Ydv2p8c7WU1NKX6/z6tY\np9WQMwTAkoZGGhoqztofCr8yAr2yDEhB2kie89jzne5JbnSfcqd7lTvdq9wU4j55+sz7DMYZXy8E\nvgkcBn5rWdabbdv+7Vgv7u0tnpW3uiN9AJQZJXR2Dp61Px57pYvccF1cxyBN/JzHns8aGip0T3Kg\n+5Q73avc6V7lxuv7NNYHAy+7zdvItLSzFgAnRr7uAo7Ytn3Qtu008AhwoYexzKihVATXhflVNRMe\nGwwauMkQjqlucxERyY2Xyfth4F0AI13jbbZtDwLYtp0CDlmWtWrk2MsA28NYZlTMiUIqSG1FyYTH\nGgYY6SD4NWBNRERy41nytm17K/CcZVlbyYw0v9WyrJsty3r7yCEfA+4a2d8P/NqrWGZa0hjGTQap\nLA3mdLyRDoCZJp5W61tERCbm6TNv27ZvO2PTzlH7DgDXenn9Qkimk7hmEp9ThWmOX10ty3SCOEDf\n8ADzyuu9DVBERIqeKqxNs4FEZuBCkImrq2X53EwLvX2wz5OYRERkblHynmadkZHqar6Jq6tl+YxM\n8u4aqcwmIiIyngmTt2VZH7QsS5P9cnRypLpaWWDi6mpZ2frm3SPFXURERMaTS8t7PfCiZVl3W5Z1\nndcBFbuuoUzXd1Uo9887ATPT8u6LaU6liIhMbMLkbdv2XwMrgLuB/8eyrK2WZX3SsqyJJzGfh7qH\nM63nupKJS6NmBX2ZlvdAXMlbREQmltMzb9u2HeAgcJzMKmCXAU9YlvU/PYytKGWXA60vq875NdnF\nSSLJIU9iEhGRuWXCqWKWZb0P+HOgHvgu8Drbtnsty6oGtgC/8jbE4hJJDoEB8ypz75gIBzIt7+G0\nkreIiEwsl3nerwdut237ydEbbdvusyzrG96EVbyG0xFcn8H8qtxb3qGgiRv1E2fYw8hERGSuyKXb\n/EvAm7PfWJZ1l2VZ6wBs277Lq8CKVYJhSAapLg/l/BqfD9xUkJQR8zAyERGZK3JJ3ncA94/6/nvA\nnd6EU/xSZgwjHcJn5j6FPlvf3DFjOK7jYXQiIjIX5JJh/LZtP5H9ZvTXcrp4OgFmioA78YIkZzLS\nQTAgmlTXuYiIjC+XZ979lmV9CHiMTLJ/I6A5TeeQneMdNnMvjZplupn65t3RfsqDuVdnExGR808u\nLe8/JzM17KfAvcCqkW1yho6BTPIuMSeffH1uZsR5h+qbi4jIBCZsedu23Ql8YAZiKXrd0UyHRFlg\n8i3vACGSQNeQSqSKiMj4cpnn/V7gk0AtcGqNS9u2l3gYV1HqHc4k73y6vQNmpuXdM6zkLSIi48vl\nmffnybS8j3gcS9Hrj0UAqArnvihJVnCkvnm/6puLiMgEckne+23bftzzSOaAwUQmedeEJ78IW9A/\nUt88oeQtIiLjyyV5b7Us65/IjDZPZTfatv2oV0EVq6FkFICasskn72x982gqOq0xiYjI3JNL8r5p\n5P9Xj9rmAkreZxhOR8EHDeVVk35tKODDdQxijpK3iIiML5fR5q8BsCzLsG3b9T6k4hV3Yrgm1Jfn\n0W0eNGAwSMJUkRYRERnfhPO8LcvaYFnWdmDvyPe3W5Z1leeRFaEkMUgFKC8JTvq1Ph+QDpJS8hYR\nkQnkWtv8/cCJke9/AnzNs4iKWNqIYzhBTMOY+OBz8DthXDNF0klNfLCIiJy3ckneSdu2X8x+Y9v2\nPkYNXJMM13VxfQn8bjjvcwSMzGuHElrXW0RExpZL8k5ZlrWMzCA1LMv6H4wq1iIZkUQUDJcA+Sfv\nsJlZ0ERV1kREZDy5jDb/BPArwLIsqx84DPyZl0EVo85IJuGGzPyTd6m/lL6Rc62sm6bARERkzsll\ntPmLwHrLshqAuG3bahaeQ1ekH4AS3+TrmmeVBUrBhZ4hFWoREZGx5VLb/AeMdJmPfA+Abdvv8y6s\n4pNNuPksSpJVGSqHGPSpRKqIiIwjl27zTaO+DgKvAVq8Cad49Y7UNa+Ywlrc1SWZ5J2tkS4iInIu\nuXSb333Gpu9alvUbj+IpWgMjreXKPBYlyaotrYBeiGi0uYiIjCOXbvMzR6QvBlZ5E07xiiQyZU1r\nSydfXS2rvqwSUH1zEREZXy7d5ikyz7yz08P6gS97FlGRGkpFwZxa8m6srAZgOK0qayIiMrZcus1z\nmQt+3htOZ5J3Y3l13ueoLSvFdQwSbmwaIxMRkbkml27zfxhvv23bn5u+cIpXwo3hulBfkX/L2+/z\nYaSDpAwlbxERGVsurerFwJuAEiAEvBVYCqRH/hNGFiVJBwgHAlM6j+mEcMzENEUlIiJzUS7PvOuA\njbZtpyCzqhjwC83zPl3aiGM6k19N7Ex+QiR8A6TSKfy+XP55RETkfJNLy3tBNnED2LadAOZ7F1Lx\ncRxnyouSZAWNTH3zniHN9RYRkXPLpWm3w7KsZ4AnR75/FfDiOMefd/qHoximS9CYevIOm2EiQOdQ\n/6nR5yIiIqNN2PK2bfsW4LNAG5k1vT9PZn1vGdExUtd8KouSZJX6M+VVuyMqkSoiIueW6zSwMJCw\nbfurwAEP4ylKXSMripX68y+NmlU+Ul61Z1jJW0REzm3C5G1Z1peBvwD+fGTTHwPf8jKoYtMTzSTa\n8iksSpKVrY2u+uYiIjKWXFre19u2/Q5gAMC27S8Al3oaVZHpG576oiRZ1eHMPPHBuJK3iIicWy7J\nO1ur0wWwLMtHbgPdzhuDiUyirQrnX6Alq2akvGokqcVJRETk3HJJ3lsty7oLWGBZ1t8AW4DHPI2q\nyGRXAauZQl3zrPqyzDmGU6pvLiIi55bLaPPPAr8FHgEWAV+zbftTXgdWTIZGVgHLJt6pqK+oAiDm\nKHmLiMi55VLb/Dbbtv838LMZiKcoxdPD4IeG8qopn6umtAzXRYuTiIjImHLpNl9nWdZKzyMpYvGR\nRFs3DS1vv6nFSUREZHxjtrwty1pg23YbmYpqey3L6gYSZNb1dm3bXjJDMc56KSMGaf+01SLPLE4S\nn5ZziYjI3DNetrnPsqxryLTOLUaS9qj/ywjHjGM6oWk7X4AQMd8gyXSKgBYnERGRM4yXGQ4BQ2SS\n9/5R27PJ2+dhXEUjmUpnFiVJl0/bOYNGmLgB3UMRmlTfXEREzjBm8rZt+z0AlmV917btv5y5kIpL\nz9BQZlGS9NTrmmeFzRIGga7BASVvERE5Sy5TxZS4x9ExmKlrXuIrmbZznlqcZGhg2s4pIiJzh6cP\nVC3L+jqwkUw3+0dt2952jmO+BFxt2/YNXsbilZ6RBFvin3pd86yyQBkkoCeqEqkiInK2XFcVmzTL\nsq4HVtm2fTWZhU3OWszEsqy1wKu9imEmZFf/Kg9Mva55VmUo8/y8P66VxURE5GyeJW/gRuCXALZt\n7wVqLMuqPOOYr5JZK7xo9Y0k76rw9A1YqyrJnGswrvrmIiJyNi+7zZuA50Z93zmybQDAsqybydRJ\nP5zLyWpqSvH7Z98A9ziZYioLautoaMitSEsoHDjn9uzrlzTWQ0fm3Lmec64633/+XOk+5U73Kne6\nV7kpxH2ayUnERvYLy7JqyawPfhOwMJcX9/ZGPQpranqG+sGEMCE6O3Pr5o7Hkufcnn19yAkCMBCL\n5HzOuaihoeK8/vlzpfuUO92r3Ole5cbr+zTWBwMvu83byLS0sxYAJ0a+fi3QADwB/Ddw6cjgtqIT\nHVn9q778zCcC+WvMLk6S1uIkIiJyNi+T98PAuwAsy7oUaLNtexDAtu2f2ba91rbtjcDbgR22bX/c\nw1g8Ex9Z/atxGhYlyaoKZwa/JVwlbxEROZtnydu27a3Ac5ZlbSUz0vxWy7Jutizr7V5dsxAS7jC4\nUDmNA9b8Pj+kA6QM1TcXEZGzefrM27bt287YtPMcxxwGbvAyDi+ljTiGE8Q0pvdzkM8JkTIT03pO\nERGZG7zsNp/zXNfF8cXxTeOiJFl+QuBLkEimp/3cIiJS3JS8p2AolgB/kgDTV9c8K2iEMUyX7oiq\nrImIyOmUvKegKzKIYbiEzOmra54VNjPlVjsj/dN+bhERKW5K3lPQGckuSjJ9dc2zXlmcRMlbRERO\np+Q9BT1DmYn5ZYHpT94Vwczo9WztdBERkSwl7ynoGc60vMsD0zdNLCu7OElfTMlbREROp+Q9Bf3x\nzGCy6mmc451VE86UxBtMaHESERE5nZL3FERGEmt16fQXpa8ry5RbHUoqeYuIyOmUvKcgm1jry6av\nrnlWfVmm3Go0reQtIiKnU/KeguGRhUMaprGueVZjRTXwSu10ERGRrJlcEnTOyS4cUlc6/S3v8lAY\n1zFJjqwXLlKMHnuhdcx9N1yc02rAInIOanlPQZIYOCZh//SXRzUMAzMdImUoeYuIyOmUvKcgbcYx\nnRCGYXhyfp8bwvUlcF3Xk/OLiEhxUvLOUzLlgC+B353+VndW0CjB8KXpj0Y9u4aIiBQfJe889Q0N\nY/hTBIzpX5QkKzxSM/3EQJ9n1xARkeKj5J2njsFMzfHsAiJeKPOXAdA1UkNdREQENNo8b11DmYRa\n5vcueZcHyyCuxUlkdhtvRLmIeEMt7zz1REfqmgfLPLtG9UiJVC1OIiIioyl556l/ZMGQ7AIiXqgp\nycwfHxipoS4iIgJK3nkbiI/UNQ9Pf13zrGx980hSyVtERF6h5J2nbF3zWg8WJclqHCm7Gk1pqpiI\niLxCA9byFE1FwQd1HtQ1z2qsyJxb9c2lGA3HU6TSDgG/ScBv4jPVVhCZLkreeYqNJNTG8umva55V\nFigF1yCBkrcUj0g0yc6DXRxqHWB0bcDq8iDrV9bTPK/cs6qEIucLJe88ZRcl8XLAmmEYGOkQaTPu\n2TVEpksq7fCc3cn+Y304LlSVB6mvCpNMOcSTaTp6h3n8hTZqKkJcsqq+0OGKFDUl7zyljDik/fhN\nb29hgDArxVwkAAAgAElEQVRxX4REMk0w4PP0WiL5cl2XJ188wdH2CBWlATasrGfp/ArMUS3sgaEE\nLx7spqVtgEd3tFIa9vO265afdoyI5EYPofLguC6uL47Pw7rmWUEjjOFP0RPRoDWZvZ7f38XR9gjz\nakp467VLWb6g8qykXFkW5Nr183nLNUupKA3wm61H+Pdf7SGRTBcoapHipZZ3HoaGk+BPEHC8G2me\nVeIrJQK0D/TTVOP99UQm62BrP7sP9VBRGuD6SxaeNTDtUGzX6S8IwNpL4Mj+MNtf7qC7P8bH37OB\n8pLADEYtUtyUvPPQMzSEYbqEXO8WJckq85fRmYKuSD+wyPPriUzGwdZ+nt59kmDA5MbLFhEO5vZo\nJxCAN9wU4qlnEhw4NMAX7nmKN94UIhA4vbV+7cKNXoQtUvTUbZ6HzsHMKl8lPu/qmmdVBDMD4rqj\nWpxEZhfXdfnRpn04Llx/8QIqy4KTer3PZ3Ddq4KsWOajs8vh0S1x0mmtXS+SCyXvPHRHM6VRywLe\n1TXPqhqp4NYXU31zmV127Ouk5cQgzU0VzK/L771gGJkEvnihSesJh8efSuA4SuAiE1HyzkPvSCu4\nwsNFSbLqRiq4qb65zCZpx+EXjx/CNIwpT/syTYPXvDrEvEaTliNptu1ITlOUInOXknceekdawdUl\n3s3xzqovy1RZi4yUYxWZDbbuPsmJ7ijXrp8/6e7yc/H7DW56TYiqSoM9e1McOJSahihF5i4l7zz0\nxzPJu6Gs2vNrZUukDqc1VUxmh2Qqza+ebCHgN3nrNUun7byhoMFNN4QIBuCppxN0dWsKmchYNNo8\nD5FkBHwwv7LW82tln3nHXZVIldlh8/Nt9AzEeeNVS6itzH/Gxb5jfefcvvICg5d2Bdj0WILXLI9T\nVe59PQWRYqOWdx6GnUwXdmNFjefXytQ3hyQxz68lMhHHddm0/RjBgMmbNjZP+XwpN0FXqo0DsZ08\nH93C3tg24pUtrLt0mGjU4d9/tYe040xD5CJzi1reeUgQBRcqg94/8zYNE9MJkTLjpB1HKzNJQe09\n0ktXf4xr18+fclGV3lQHLYk9uCPLlwSNMFEnQtQZBP8hqtfPw96zjvuePMzbX718OsIXmTOUvCfJ\ncV3S5jA+J4zPnJla4wHCpANRItGkuhBlxj32Quupr7e80AZARWngtO2TNZDuoSWxBwOTBYGlVPsa\nCJulpNwkA+luIr52Ommn9KJBfrszirWkmrVLvX9MJVIslLwnKTKchECcgOvdUqBnCpolxI1++ob0\n/E8KJ5ZIcax9kOqR1cLyNZQe4GB8F2CwMrSeCt8rj5/8RoBafxM17jz8/hZOcJjgmme447Eh/mDN\ndZSEMn+ybrh44VR/HJGipuQ9SR39gxi+NCWO93O8s0p9pQw60DHYT/O8mfvQIOePXFrRh1oHcFxY\nuajqrPW4z6pfPoa4E2V/fCcOaZYHLzotcY9mGAYLgssp91VzYHgPzuIXefCgy/plDRgG+FuPqXSq\nnNf0AHWSTgz0AFAW8P55d1b5SCW3zsi5R+eKeM11XfYf78c0DJYvqMr7PMeTB0mTZEnAosbfMOHx\nlb5arJKLMRw/yXm7OXiyK+9ri8wlSt6T1BHpBaAqOHMt4KpQ5lrd0f4Zu6bIaJ19w/QPJVgyrzzn\nxUfONJQeoC/dSZlZSb1/Qc6vK/NVsCK4ARyTvopdnOjvyev6InOJkvckZRNobUn+rY/JqivJFIPp\njSl5S2HsP5753Vu1OP/f+7bkIQAWBlac1e0+kepQFYvc9YBBq7mLjkG9F+T8puQ9SX2xTF3zbNnS\nmdBYnnkuOJDUymIy81JphyMnBykvCdBUm99KeoPpXgacHirMmjGfc0+kqbKGqshaDF+azUeeJpJQ\n1UE5fyl5T1IkmSmN2jQD1dWymioy1xpKaXESmXltXUOk0i7NTRWTbjFD5nl566hW91SsnNeIr2cp\nTiDKPz/9PRxXBVzk/KTkPUnRdKa6WtMMVFfLyra8Y64WJ5GZd7Q986GxeV5+gzQHnB6GnH6qfPWU\n+aY2VsQwYE3DMhiopzN9lP968VdTOp9IsVLynqQ4ma666tDMDVgrC5SCY5I0lLxlZqUdl2MdEUrD\nfurynNvdkTwOwILAsmmJKRQyeNWCK3FipWzrfpodJ/ZMy3lFionmeU+CO1JdzXQCBHxTKw05GYZh\n4EuXkvLFcF03r65LkXyc7I6STDmsWFiZ1+9dwokz4HRTZlZSalZMW1zNC8Mc67yMo86T3LXnXrri\nr6U08MrzeM0Bl7lOLe9JGIqlRqqr5TdoZyqClEIgzmAsPuPXlvPXkfbMGI/mefkl3u70CQDq/POn\nLSbIrEhWXxPA12HhmAke3Pcs9tHeMVcqE5lr1PKehK6BIQx/krAz88m71Cxj2IDW3h4qS3KfIysy\n2pOtz5xz+6FYH8vDF522zXFcjrVHCAd9NNSUTPparuvSnTqBgUmtb15e8Y7HNOGCpvns7ushXt3B\nseHDLCmdnq55kdlOyXsSssUhyvwzV10tqyJYSXcaWvu6WbNAyVvGNlaCnqyO3mHiyTSrF1dh5tFl\nHnH6ibvD1Pqa8Bne/KkpKTVYHFnDsfgAncEWqpM1QLUn1xKZTTxN3pZlfR3YCLjAR23b3jZq32uA\nLwFpwAY+YNv2rJ730TFSnrRyBqurZdWWVHE4Au2Dqi4lMyPbZb4k3y7zVGYFsvpp7jI/U1Ojj96W\ndQw1bufA8EtclNaHW5n7PHvmbVnW9cAq27avBv4C+NYZh3wHeJdt29cAFcAbvYplunRFM8m7tmTm\nk3djWc1pMYh4yXVdjrZHCAbMvAqzpN0UvekOgkaYctP7lvDqJZUYHStx/TEea9mO67qeX1OkkLwc\nsHYj8EsA27b3AjWWZY3OepfZtn185OtOoM7DWKZF70h1tbrSme+Wm1+VuT19cVVZE+91D8QYjqdY\n1FCOaU6+y7w33YGDQ71//ozMjvD5wGpYgjNQQ0/6JA8efNzza4oUkpfJu4lMUs7qHNkGgG3bAwCW\nZc0HXg/c72Es02IwMVJdbQYLtGQtqs4k76HU4IxfW84/rZ2ZmgKLGvMb39GdOglAnc/bLvPRysth\nvnMhbjLAb448wPGBthm7tshMm8kBa2d9/LYsqxH4NfDXtm13j/fimppS/P78VjOaLjE3U6BlzdLF\nNFTlP2c1FD73HPGGhrHPWVUbht9DjKFxj5trzqefdSpG36eKgckXUwmFA1SUv/K6E91RDANWN9cS\nCkz8vgvxyu90wokTifZR4a+monRm//2WLwuQPHwJPVW/544dP+Df3v45gv7gacfodyp3ule5KcR9\n8jJ5tzGqpQ0sAE5kvxnpQn8A+Kxt2w9PdLLe3sIvQjCczpSJdKM+OhP5t4DjseQ5t3d2jn9OIx0k\nZURp7xjIa/RvsWloqJjwnsjZ92lwMDbpc8RjSQbJvG44nqKjd5h5NSUk4kkS8XP/vp75+qyOZKbF\nW2U0jPm77qUbNizkvueXMVjTwlc23c1fXfZHp/bpdyp3ule58fo+jfXBwMtu84eBdwFYlnUp0Gbb\n9uif8KvA123bftDDGKaN67qkzGFwfIR9oYLEkC3UMjCUKMj15fzQ1pXpMl+YZ5d5b7oDgBpf47TF\nNBmhkMEHr3gXTrScF/t3sK1tV0HiEPGSZy1v27a3Wpb1nGVZWwEHuNWyrJuBfuAh4H3AKsuyPjDy\nkh/Ztv0dr+KZqmg8Bf44AbekYOVJS81y4kYfJ/oGqC5vKEgMUnhnzuOuGAjn1doey/Hs8+6Gskm/\nNunGiTh9lJtVBM3CfMjdd6yP5eEEy1LXc9h5gLv3/ITOE0FKfeW8+3UXFCQmkenm6TNv27ZvO2PT\nzlFfF+adnaeewRgE4oTcmVvH+0yVwUp6kyOFWhYpecv0cxyXtq4hysJ+qsqCE7/gDL2pzBjVQrW6\nR7u8eQXt+9YRa3iRLd3388aGdxc6JJFpo9rmOTrZ34thQJlv5qurZdWEMx8c2gdUqEW80dE3TDLl\nsKixPK8ept50OwDVsyB5m4bBa5dci9vfQMR/gp392yZ+kUiRUPLOUftgLwAVwcKNvmwoy8wv74z2\nFiwGmduyU8QW1k++yzzhxIk4/ZSb1QXrMj9TWUmAy8teh5sMsj/1LHvbWwodksi0UPLOUbayWbb1\nWwjzK1WoRbzV2hnBZxo01U2+qlpfgQeqjWV5YwNNQ1eD6fClR79DPKWV+aT4KXnnqHc4kzDrywqZ\nvGsBiKhQi3ggMpykL5Kgqa4Uv2/yfxp6RrrMa/yFH49xKLbrtP8WNPkwe5YQM/v4l2d+WujwRKZM\nq4rlqC/eDyXQVFFbsBiqR1r9MWeoYDHI3NXamaljkE+XedwZZsgZoMKsIWDMji7z0UwTLqhfxp6h\nHlrYxb9t3sSamjWnHXPDxQsLFJ3I5KnlnaOBVKbbfEl14boEywNl4JqkfVESyXTB4pC56ZUpYpMf\nlJkdqFbrn/51u6dLSdjH8tB6XMdkd3Iz/bH+Qockkjcl7xw4rkvcGATXoK6kcC1v0zAJuiUQiNM7\nqOd2Mr5YKs7hgaN0D/eQdsdfbTedhpPdUarKg5SXnrt873h6Uh0YGFT7Ct9lPp759RU0Dl8K/iSP\ndv8Wx5nVqxCLjEnd5jnojyQgOETALcNnFra+eomvnDgddPUPMy+PpRrl/HByqIOnT2wjls58yDMN\nk7pwDZfPu5jq0NnjNvr7DNKOm1eXeX+qm2E3QpWvHr8x+cQ/065beDX3tR4nVXaSx9ue4IZF1xc6\nJJFJU8s7B8e7+zCCiRlZl3giFYEKDNOlrV/TxeRsjuvwYtdLbD7+JPF0gjW1q1lVvZyqYCWdw91s\nPvYkA+eoy9/bk/lTkE+X+bGEDUCtb/Z2mY/mM01e2/gmSITpKNmJ3X2w0CGJTJpa3jlo6R5Z3jA0\n80uBnqkmXMXxhAq1yLnt6HiR/X2HKAuU8qr5V1I/6jHP/t6DbO/YyaPHnuCmxddTHsy0sl3Xpafb\nJOA3aawpmdT1XNflaGIfJj6qfPXT+rN45VAsU+t8kXkhx9jBztTvSPb14W89BsC1CzcWMjyRnKjl\nnYPWgcz81aZZUE+8sSzzx7gjquQtp+sa7mF/3yEqgxW8sfm1pyVugFU1K7i4YR3DqRiPHn+CaGoY\ngP4Bl3jMYEFdKaY5uapqvel2hpx+qn31+IzCPlKarHmVVVREV2AE4rw8+BLJlJ5/S/FQ8s5B53Bm\nqfElNYXvFmyuzqyy2psYd/lzOc84rsO29ucBuGLeJQR9565LvqZ2Nevq1jCUjPLMiedwXZdjrZmZ\nCwvz6DI/Gs90mdfM4lHm41lVtxhftB63rIff7X6p0OGI5EzJOwf9ycw0seZZkLwXVWZiiDh65i2v\nsHsP0BfvZ0XVUhpLx+++Xld3AQvK5tEe7eBgfwvHj2eT9+QGq6XcJEcSLxM0Sqg0CzcLYypM02BN\n9RpIhOkL72P7vhOFDkkkJ0reE3Bdl2Ey1dUaJvijOBPqS+rANYibA7iuW+hwZBYYjEfY1bWXkC/E\nhoZ1Ex5vGAZXzLuUgBlgR8cu2vsilFc4lIQmNwTmSHwvCTfGitBFmEbx/ikJ+QMsC14IGOxLPMfO\nI8cLHZLIhIr3HTdDBqNJ3OAQfqeE0BhdkTPJZ/oIOOUQHmIolip0ODILPHv8BdJumksaLsr5d7Q0\nUMJljetJu2n8y3ZTXTu5572u67Iv9jwmPlaG1+cT9qxSW1rFivCFGIEE333xHnoGo4UOSWRcSt4T\nONETwQgOU2oWrqb5mcqNGgx/kuM9GrR2vhtMRGjpPUZNqIqllYsn9dqllUsIxxvxVfbgqz82qdee\nSLYQcfpYErQIm5OfGz4bXdm8igqnCbesmy8/+hNSaQ1gk9lLyXsCh7pOYhhQGyz8NLGsmlDm+eKh\n7rYCRyKFZvcewMXlgtrVk15/23UhemAtpAJ0+fYTSedeLnRfLDM4bnX4kkldczYzDIPXrbqMoFNB\npGIvdzz6UKFDEhmTkvcEjvdnponNKyv88+6sRRWZQWtHek8WOBIppHg6zqH+I5QHS1lSMflFNTo6\nHeJDYepi63Bw2Db0u5zGUfSmOuhMHWeefwlV/tnzvpgOIX+Ij13+F+D42cfj/Oq5nYUOSeSclLwn\n0DHcBcDi6sKPNM9aWZ/5Q31yqKPAkUgh7e9rIe2mWdd4QV4Dxo4cy4wyv2DeYqp9DXSlWjkQnzhZ\n7YvtAOZWq3u05uoF/OHKd2P40jzU8Qt2HNQANpl9lLwn0JfITMlqrmkqcCSvWNmwAIC+lJ55n6/S\nTpp9vQcJmAEuaFgx6de7rsvR42n8flg438eSoEXQCLMr+hSD6bGnIR6Jv8zRhE2Vr555geap/Aiz\n2quXXsJVtddihIb5z933cLxroNAhiZxGyXsCUTfzpp1N3eaVwQoMJ0DCN6ClQc9TLQNHiafjrKxe\nRtA3+cVA+vpdBgddFi3w4fMZBIwgl5a+hjQpnhj85TkTeG+qne1Dm/AbQa4uf9Okn7HPdvuO9bHv\nWB/fe+ZBvvfMgySHA4STDVDRzf9+7AdsOvh0oUMUOUXJexyR4SROIILpBCkLzJ4VvAzDoJRqjFCU\n1q5IocORGea6Lvt6D2JisLpm8q1ugKMjXeZLFr9S0nRRcBUXlmxkyBng0YH/S0/qlTEVMWeIpyK/\nwSHNxrI3UuGbPQM4vWIYBhdUrsGXKsOtPcIDzx1QCVWZNZS8x9HRG8UIDVNqVBY6lLPUheowTBe7\nXSPOzzedw930JwZYVLGQUv/kFhLJOnIsjWHA4oWvJG/DMFhbchWXld5Iwo3x2MDPeSbyAFsGfsGm\n/h8z7ERYV/Iq5geXTdePMuv5DD9W2UUYjp/h2t3c+dDjKo4ks4KS9zhaujswTIfqWTRNLGvhSJnU\nlh6Vczzf7O87BMCq6uV5vT4adejqdmhqNAmFzu76Xh5exzXlbwHgWGIfHaljpEiwMrSBC8KX5x94\nkSrxlXLtwisxDJd9/k385KkXCx2SiJYEHc+xvnaACWtFF8KKugU83Q0nIu2FDkVm0EBikOODrVQF\nK2koqcvrHEePZ7p+R3eZn2lBcDlvqf4ASTdO2CzFZ5zffyoWVTaxPraeF3tfZEv/L1mwu5pXr5u7\nA/Zk9lPLexztkcw0sUVVhV8K9ExLTq0uphHn55On27bh4LKyelneA8aOHM2U1V2yaPwlPINmiDJf\n5XmfuLMubFzJVfUbMUuGuPfgj9hzuLPQIcl5TMl7HD3xTGKcTdPEshpK6sGFpH+QyHCy0OHIDHBc\nhyfbnsVn+FhauSSvc0SHXdpOOjTUm1RU6O2fq+xI9MXpK6lzl2FW9HLH9h/y8y0HeOyF1kKHJ+ch\nvXvHkEo7DJD5ZL2kckGBozlb0BcgbFRghodo7dSI8/PBS902PbFellYuzmt6GEDL4RSuCyuWjd/q\nlnMzDIPra99EmdOIWXuCR9seYTCaKHRYch5S8h7D4ZMDGGV9hN0qSmfRNLHRaoK1GME4LR3qOj8f\nPNGamWe8Ms+BagAHWzKjzJc1qys8H4diuzgS38vKstX40qUw7yAPHNzCwwc0B1xmlt7BY3jhWAuG\nL83C8KJChzKmBeWNnOg5wsHuE0B+831lYk+2PjPmvmsXbpyRGDqiXezptllWuYTacHVe5+jvz4wy\nX7TApKRkbhVYmWl+I8iasvXsiT6Hu+Al7n+6hFctvJzykvx6REQmSy3vMbzc3QLAmsb8WzleW1ab\n6c5vi2i62Fz3eOtWXFxuWHRN3uc42JIZqLZimT6zT4eQWcrqkvUYmCSanudLv/wd0Viq0GHJeULJ\n+xxc16U9nhmEctG82duiXVmTKZbR65xQ4Yg5LJaK8XTbdqqCFVzceFFe53Bdl4MtmVrm400Rk8kp\n91WxIrQOw3TprX+Cr/xqM7GEErh4T8n7HHoG4qRCPZhugAXls2+kedbC8iZ8bhC3rJvu/lihwxGP\nPHtyB7F0jOsWXo3fzK/V3NnlMBhxaV7sIxBQl/l0qvbXc/X8yzFMh47qx/jKLx5nOK4ELt5S8j6H\nPcdOYpYMUeubl9dSizPFNEzqfAsxw8PsPaEyqXOR4zpsOb4Vv+HjmoVX5X2egy2ZWubL1WXuiaVV\ni3mv9Q6MQJKT1Zv58i+2qAtdPDV7M1MBvXjiIAAramZ/BaUVlZmu810d+wociXjB7jlAe7SDy+Zd\nTGWwIq9zJJMuh1pShMOwcL7e8l7Yd6yPVNciNpRcjxGM01G7mdvv2cRD244WOjSZo/ROPocjkcwb\nbsP8VQWOZGJXNa8F4NBAS4EjES9sPv4kANcvelXe5zhwKEU8ARes9mOa6jL30uqSi7m49AaMQILh\nRU/wwAu76B2MFzosmYOUvM8QT6SJ0AEUScu7dhGmE2TI305/RH8k5pKjg8fZ0/0yy6uaaa5cnNc5\nXNdlz94UPhPWrNY0ppmwKryBS0tfi+FPkmjeyhd+9jAnuocKHZbMMUreZzjY1odR3kcJVZQHygod\nzoRMw6QxsAgzNMzWfYcKHY5Mo98eehiANy17Xd7nOHo8zcCgy4rlPs3t9tih2K5T/xlAc3ANhi9F\nbNFT/NOv7udga3+hQ5Q5RMn7DDuPH84UZymdvcVZzrShaTUA21tfLnAkMl1a+o+wu/tlVlQt44Ka\n/B/f7HkpM2jqwjVqdc+0+sB8VoYvwuczSDdv4ysP/4qn95wsdFgyRyh5n2Ffz2EALmycvfO7z3TZ\nwjUAtMWOkkimCxyNTIffjLS6/2D56/NePayrO83JDoeFC0xqqvVWL4QqXz03LbmOEl8J/ubdfH/n\nz/nJ5n04juoyyNToHT1KKu3QkTwOwNqG2VtZ7Uzzy+bhJwTl3ew90lvocGSKDvS18HLvfi6oWcWq\nmvw/RO4eaXWvU6u7oOpKavnUlR+hPtSAv+kIm/t+wdd/8XutBihTouQ9yrP2cdzKE4TcilldnOVM\npmHSXNaMGYrx+wMadV7MXNfl14ceBOAty1+f93l6+xxajqSpqTZYoOlhBbXvWB8v7Ytzbem7aPIv\nx1fVw4Gy3/KpH9zPvmN9hQ5PipTe1aM8sG8rhs/hVfOvnNXFWc7lkvkXALC7e59KpRaxrW2/50Bf\nC+vq1rCsKr/ZDq7rsvXZBK4Ll18SzLvbXaZXwAhybcVbWBu+CjMYw12xlX9+7Kf88omDpB2n0OFJ\nkVG5pRGtnRG6/fswXZM3rMx/8YdC2dB4IT/b/2sSVQc5cnKQpfMrCx2STFJHtIufHfg1Jf4S/sh6\ne97nOXAoTXuHQ/NiH4sXvVLHXK28wjMMgwtLN9IYWMTTgw8SX3iAB3u6eO6eq/jL11/Jknn5FeKR\n809xNS89dN/zOzBLIywtWUVFsLzQ4UxabbiGFaVrMEsj/G7fc4UORyYp7aS5+6Ufk0gn+CPr7dTk\nuexnPO6y7bkEfh9cdYWedc8Wo6eRHYrtYjDdywUll7KodCG+ij665z/MF3/3I376mE1cg04lB0re\nQCyRYtfA8wC8xbquwNHk723WTQC8MPAs0ZgGwxSTh448yuGBo1w+72Iun3dx3ufZ/nySWBwu2RCg\nvExv79nMbwS4dtGVXLdgI+WBcvwLDrI5ei+fuvcXPLWrDUePv2Qc6jYHtuw6AtUnKKWKC2pnf0nU\nsSyvXUSD2Uxn2RF+8NQz/NWNxftBZDaIpxMc6DvEjo4XGUwMknRSJJ0UJgaVoQqqgpXUhWtZVbM8\n79W+AJ5qfZb7WzZRHariD1e/Le/zHGtNY+9PUV1lcOEavbWLgWEYLKpYwNtWvpn7Dj7E461PkVy0\njR+27OP+3RfzR1dtZN2yWo1bkLOc9+9w13XZdOhpjDqHVy/eWPRvkvde9Ea+tfPbvND/DK1dl7Kw\nfvZXiZtNXNfl5d79PHbsKV7u2UfKPb0LM2AGSLtpeuKZ58c7u/ZQ4g+zvv5CLmm8iLW1Fj4zt/Wy\nXdflNy0P8+DhRygLlHLLRe+jNFCaV9wdnWk2b4nj88G1VwdVw7xIZMchpDq7aOQy3lC1nBciT3Gy\n/CB95Vu4c9dOan+/lrdt2MjlF8zDLPK/TzJ9zvvkfd+2PQxV7cV0fbymeWOhw5kyq24F84ILaa9u\n5b82/57b3nVD0X8gmQnJdJJt7c+z+diTtA1lqmAtLJ/P2lqLlJOiLlyD3/RjGAaO6zKUHKI/PoDf\n5+eFjt08e/I5nj35HOWBMq5ouoSrmi5jUfmCMe99T6yX+w4+yLb256kvqePWDe+nsbQhr9h7etL8\n7tE4aQduvD5IY0NuHx5k9jgU23Xq64WhZqqdGrqcY3RVtNNfsZW7Wl7g3p3LuW7JFdy4fjlV5aEC\nRiuzwXmdvLfvb+XBzp9jliR565L/SXlwbrRS32bdxLd33c3Rkid45uXlXL1m9i+wUigDiUEeP/40\nT7Q+TSQ5hGmYXD7vYl67+LpTi4E82frMaa8xDYOKYDkVwXKuXbiRd6x8C0cGjrGt/QW2j3wA2Hzs\nSaqCFayuWcWqmmWU+EswMIilYmxvfwG79wAuLs2Vi/nQ+j/Pe5DkwIDDg48MEk9kWtxLFp/Xb+k5\no8ys5JLma+iPD3BioIcXul8g0bibTdE9/O6hRhb5V3Pd8g1csXoBJSH9m5+Pztt/9cMn+/ne7nsx\nq4a4tOZKXl+E08PGclH9Wi6vu4LtbOOHh35Aafj9bFi2oNBhzSrHB9vYfOxJtrc/T8pNU+ov4XVL\nbuD6Ra+a9Ehv0zBZVtXMsqpm3rHyzezpfpkdHS9i9xxgW/sOtrXvOOs1y6uWcvX8y7li3iUEfJMf\nFe66Lnv3pdj+XJJUGi6/JMDqleft23lOynaprwxfx+Kayzk0vJcD0d3Eatppo50ftz/FvfvrmOdf\nysXzLuDKFcuYX1emnrbzhKfvdsuyvg5sBFzgo7Ztbxu17ybgn4A0cL9t21/wMpbRBqIJvr7l5xgN\nHaPMmboAAAnTSURBVCwILeXmDfnPqZ2NDMPgz9a/k56twxxiN/++6y7e3PEe3nzlyvP6jT2UjLKz\ncze/P7mD/X2ZFdgaS+t5zaLruGr+ZYR8wSlfw2/62dCwjg0N63Bdl7ahkxzuP0rKTWeK5xiwpmYV\n88oa875Gb6/DM9sTnDjpEArCTa8tY0GTRibPVdku9ZAvwIUVlxB1Bukc7mSQbuLVnXTSye+GtvHQ\ntjD+WD2NoSZW1jSzYeEyls2rVct8jvLsX9WyrOuBVbZtX21Z1hrge8DVow75FvAGoBXYYlnWz23b\nfsmreEZr7RwiZgxQa9bz8StvznmAUTExDZOPv+pP+NftP2IvL3J///d59v7VvGnVq7liRTN+39yf\nRpRyUhwbbOVg/2Hs3gPYPQdIjwxAW12zkhsXX8faOsuzanqGYbCwfD4Ly+dP+Vz9Qwl22B08uesE\nLSf+//buNUauuozj+PfMzM7uzm5Ll20LFNpC0T5owHhDSQrIVQQbEcGQSCFaDBi8QIjxDQqCL7wg\nXiKEmAii+MIQUcBAgCCJlgBGifEC8mihtLSFwnbZ28zOzuUcX5yzdHahW6TdOTM7v08yu/85l5ln\nd8+e5/zPnPP8ywCsPCLLuhPyHLK8m/Hx8n6/h7SHQmYRq/sWAWtY3rWKbcUt7ChvZTS3k/Cg7bzM\ndl6e/CuPbYbw6R5ytX76s0sYzA+yvLCUFYuWcdiSQVYsWcxB/T26CG4vwihs6Uqb83lIdjpwD4C7\n/9vMBsxssbuPmdkaYNjdXwQwsweS5ZuSvN+1eoAfLL+C3p5cS/9x9lcmyHDFBz/Dvf9dxqNbH2O4\n92nu3PYMv3p2gIPzS1nes5ylfQMM9vXTny/QncvTlc2Sy2boymbIZbPkMhmCICAIMgQEBMSnUWLR\nnq9RQztoaCfzImb2DCOi1xeYOW/Puq+Fo4y8VpzxftPtMAqphlWqYZWpeoVirUipWmK8OsFweZih\n8m52l3dTi2qvr31Y4VDec/BxHDt4HAd3DxBFMDrxxvvhZ5eXnSjOKl3ZMHtoZPLNJs/42d9sfuOT\nMIqo1kIq1ZCpap3xUoWxYoXRYoWdQ0W27hpnZKICQBDAysMzrH1njlVHZDv6TIrAK9Vt9OSzHJ1f\nQxQdRTmcZHRqjFqmTLleZiI3Qr17iDGGGGMzW8pAGXgVojCAWp5s2EM26iFPD/lMN/lsnp5cN4sL\nBYJ6lnwmT1c2F+8bMlnySbsrm6MrlyOfzZIJsmSTbTEIMgQBBARkgiBuB/HeI25nCIi35UzDstOm\n243bduP8Rm93+w+jkHpUpx7VqdQrlOqTlGolRqZG2FHcwfbidir1Clce9xUW5+NqlWEYUQ9DavWI\ncqVOqVylNFWj9I+XeOb53WzbNU69HnH9xuMp9Mx/gaT5TN6HAo2lvl5Npo0l319tmPcK0NQxOPt6\n9/8UaTvIBBnOW3sW56w5lQeefYInd/2Zib4hhoNhhsP/wDjxY4GJ6lmich/hxBLC8SWEEwM8X+nl\neeCeA3iMeBdPHLDX2ptCIWDl4RkOPSTL0UflKBSUsOWNgiCgN1ugtzDzdsMwqjNRnaQ4NUmpXmIq\nLFOLKtSDCmGmQr1rgjA7ShUoNq5Yamb0rSWqdFMfG+T62/4G4VtLk309OY5ZNUBXrjkdwmZ+GDLX\nHmefe6NlyxYtiD3WVWe+/SIc++uyw9ZzGetTe39J2dsv3PYWHDOfLy4is8znIcJO4h72tBXAS3uZ\nd3gyTURERPZhPpP3w8AFAGb2fmCnu48DuPsLwGIzO9LMcsD6ZHkRERHZh2A+x342s+8AJwMh8EXg\nfcCou//OzE4Gvpssere7f3/eAhEREVlA5jV5i4iIyIG3cO+TEhERWaCUvEVERNqM6uY10VzlYmUP\nM/secBLx9vltd/9tyiG1NDPrBf4FfMvd70g5nJZlZhcBXwNqwLXufn/KIbUcM+sHfgkMAN3A9e7+\nULpRtRYzOxa4F/ihu99sZiuBO4Es8R1VF7v71HzHoZ53kzSWiwUuJS4PK7OY2anAscnv6WPAj1IO\nqR18HRhOO4hWZmaDwHXAicR3t5ybbkQt67OAu/upxHcL/TjdcFqLmfUBPwH+0DD5BuAWdz8J2Axs\nbEYsSt7NM6NcLDBgZovTDakl/Qn4dNIeAfrMbOEVnz9AzOwY4N2AepFzOwN4xN3H3f0ld78s7YBa\n1BAwmLQHkueyxxRwDjPrkpwC3Je0f0+8rc07Je/mmV0SdrpcrDRw97q7T1dpvJR4xLl6mjG1uJuA\nq9MOog0cCRTM7D4z22Rmp6cdUCty918Dq8xsM/GB9FdTDqmluHvN3SdnTe5rOE3+CrD/IxG9BUre\n6VkQ5V7ni5mdS5y8v5R2LK3KzC4BnnD3LWnH0gYC4h7lp4hPDf/czPQ/OIuZbQC2ufs7gNOAm1MO\nqd00bZtS8m6eucrFSgMzOwu4Bjjb3UfTjqeFfRw418yeBD4PfMPMmnLKrg3tAh5Pek7PEQ/Hsyzl\nmFrROuAhAHf/O7BCH1vt00Ry0Sg0sdS3knfz7LVcrOxhZgcBNwLr3V0XYc3B3S909+Pd/QTgZ8RX\nmz+Sdlwt6mHgNDPLJBev9aPPc9/MZuDDAGa2GpjQx1b79AhwftI+H3iwGW+qW8WaxN0fN7OnzOxx\n9pSLlTe6EFgK3GVm09Mucfdt6YUk7c7dd5jZb4Ank0lfdvdwrnU61E+B283sj8T54Qspx9NSzOwD\nxNeZHAlUzewC4CLgDjO7HNgK/KIZsag8qoiISJvRaXMREZE2o+QtIiLSZpS8RURE2oySt4iISJtR\n8hYREWkzSt4iHSypqCUibUbJW6RDJZWzrk07DhH5/6lIi0jnuh1YbWYPExcvmR6sYzuwwd2rZrYR\nuIp4IJ1NwBnufqKZXQlsAErJY4O77276TyDSodTzFulc1xEn5XOIE/BJ7r4OWAKclQxZeyNwpruf\nDqxtWPcG4hK2HyEec31FUyMX6XBK3iIdzt1rQB3YlJTFfC9xidq1wFZ335UsenfDarcBD5rZNcAW\nd/9nM2MW6XRK3iIdzszWARuBjyY96U3JrAxxHf5prw9Q4e5XA58EhoF7zOzsJoUrIih5i3SyEOgC\nDgFecPdiMpLUCUA38BxwtJkNJMufB2BmA2b2TeBFd78VuAX4ULODF+lkGphEpEOZWQ54CqgRH8gX\ngaeBvxBfhX4G8AniXvnWZNkT3f0UM7sJOAV4DagCl7p7U8YxFhElbxGZg5ldDNzv7sNmdjVg7n55\n2nGJdDrdKiYic+kHHjWzUeIe9udSjkdEUM9bRESk7eiCNRERkTaj5C0iItJmlLxFRETajJK3iIhI\nm1HyFhERaTNK3iIiIm3mf3Una7/Jg34kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f85589a5320>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(np.log2(np.array([x for x in label_frame.values.flatten()])+1), label = 'peakScore')\n",
    "sns.distplot(np.log2(np.array([x for x in quantiled_label_frame.values.flatten()])+1), label = 'quantile normed')\n",
    "plt.ylim(0,1)\n",
    "plt.xlabel('tags')\n",
    "plt.ylabel('frequency')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "score_threshold = 1\n",
    "labels = (quantiled_label_frame >= score_threshold + 0).values\n",
    "target_indices = quantiled_label_frame[quantiled_label_frame.max(axis=1) >= score_threshold].index.values\n",
    "index_label_dict = dict(zip(target_indices, labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "c57bl6_il4-24h    33518\n",
       "c57bl6_kla-1h     21776\n",
       "c57bl6_veh        31550\n",
       "dtype: int64"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(quantiled_label_frame >= score_threshold + 0).sum(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_classes = quantiled_label_frame.shape[1]\n",
    "labels = [index_label_dict[x] for x in target_indices] + [[0]*num_classes] * len(negative_seq)\n",
    "labels = np.array(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "target_seq_arrays = np.array([index_seqArray_dict[x][0] for x in target_indices])\n",
    "target_seq_rc_arrays = np.array([index_seqArray_dict[x][1] for x in target_indices])\n",
    "\n",
    "seq_arrays = np.concatenate([target_seq_arrays, negative_sequence_arrays])\n",
    "seq_rc_arrays = np.concatenate([target_seq_rc_arrays, negative_sequence_rc_arrays])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "x_train, x_test, x_rc_train, x_rc_test, y_train, y_test = model_selection.train_test_split(seq_arrays, seq_rc_arrays, labels, test_size=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "label counts [25795 16140 24768]\n",
      "naive accuracy 0.2595853879440475\n"
     ]
    }
   ],
   "source": [
    "print('label counts', labels.sum(axis=0))\n",
    "print('naive accuracy', max(labels.sum(axis=0))/len(labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_convolution_multilabel_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size):\n",
    "    input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "    input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "    ### find motifs ###\n",
    "    convolution_layer = Conv1D(filters=num_motifs, \n",
    "        kernel_size=motif_size,\n",
    "        activation='relu',\n",
    "        input_shape=(total_seq_length,4),\n",
    "        name='convolution_layer',\n",
    "        padding = 'same'\n",
    "        )\n",
    "    forward_motif_scores = convolution_layer(input_fwd)\n",
    "    reverse_motif_scores = convolution_layer(input_rev)\n",
    "    print('forward_motif_scores', forward_motif_scores.get_shape())\n",
    "\n",
    "    ### crop motif scores to avoid parts of sequence where motif score is computed in only one direction ###\n",
    "    to_crop = int((total_seq_length - seq_size)/2)\n",
    "    crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "        name='crop_layer')\n",
    "    cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "    cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "    print('cropped_fwd_scores', cropped_fwd_scores.get_shape())\n",
    "\n",
    "    ### flip motif scores ###\n",
    "    flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "        output_shape=(seq_size, num_motifs),\n",
    "        name='flip_layer')\n",
    "    flipped_rev_scores = flip_layer(cropped_rev_scores)\n",
    "    print('flipped_rev_scores', flipped_rev_scores.get_shape())\n",
    "\n",
    "    ### pool across length of sequence ###\n",
    "    sequence_pooling_layer = MaxPool1D(pool_size=seq_size, \n",
    "        strides=adjacent_bp_pool_size,\n",
    "        name='sequence_pooling_layer')\n",
    "    pooled_fwd_scores = sequence_pooling_layer(cropped_fwd_scores)\n",
    "    pooled_rev_scores = sequence_pooling_layer(cropped_rev_scores)\n",
    "    print('pooled_fwd_scores', pooled_fwd_scores.get_shape())\n",
    "    print('pooled_rev_scores', pooled_rev_scores.get_shape())\n",
    "\n",
    "    ### concatenate motif scores ###\n",
    "    concatenate_layer = keras.layers.Concatenate(axis=1, name='concatenate_layer')\n",
    "    concatenated_motif_scores = concatenate_layer([pooled_fwd_scores, pooled_rev_scores])\n",
    "    print('concatenated_motif_scores', concatenated_motif_scores.get_shape())\n",
    "\n",
    "    ## pool across forward and reverse strand ###\n",
    "    strand_pooling_layer = MaxPool1D(pool_size=2, \n",
    "        strides=2,\n",
    "        name='strand_pooling_layer',\n",
    "        )\n",
    "    pooled_strand_scores = strand_pooling_layer(concatenated_motif_scores)\n",
    "    print('pooled_strand_scores', pooled_strand_scores.shape)\n",
    "    \n",
    "    # make prediction\n",
    "    flattened = Flatten(name='flatten')(pooled_strand_scores)\n",
    "    print('flattened', flattened.shape)\n",
    "    \n",
    "    predictions = Dense(num_classes,\n",
    "                        name='predictions',\n",
    "                        activation = 'relu', \n",
    "                       )(flattened)\n",
    "    print('predictions', predictions.shape)\n",
    "    \n",
    "    # define and compile model\n",
    "    model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "    model.compile(loss=keras.losses.binary_crossentropy,\n",
    "                  optimizer=keras.optimizers.adam(),\n",
    "                  metrics=['binary_accuracy'])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_classes = len(labels[0])\n",
    "total_seq_length = 200\n",
    "seq_size = 150\n",
    "num_motifs = 50\n",
    "motif_size = 20\n",
    "adjacent_bp_pool_size = 5\n",
    "attention_dim = 200\n",
    "attention_hops = 1\n",
    "num_dense_neurons = 200 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "signal_convolution_model = get_convolution_multilabel_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size\n",
    "                                                           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "signal_convolution_model.fit([x_train, x_rc_train], y_train,\n",
    "          batch_size=128,\n",
    "          epochs=50,\n",
    "          verbose=1,\n",
    "          validation_data=([x_test, x_rc_test], y_test))\n",
    "score = signal_convolution_model.evaluate([x_test, x_rc_test], y_test, verbose=0)\n",
    "print('Test loss:', score[0])\n",
    "print('Test accuracy:', score[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "signal_convolution_model.predict([x_train[0:100], x_rc_train[0:100]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attention Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_classes = len(labels[0])\n",
    "total_seq_length = 200\n",
    "seq_size = 150\n",
    "num_motifs = 50\n",
    "motif_size = 20\n",
    "adjacent_bp_pool_size = 5\n",
    "attention_dim = 200\n",
    "attention_hops = 1\n",
    "num_dense_neurons = 200 \n",
    "\n",
    "# signal_model = get_attention_multilabel_model(total_seq_length,\n",
    "#                         seq_size,\n",
    "#                         num_motifs, \n",
    "#                         motif_size,\n",
    "#                         adjacent_bp_pool_size,\n",
    "#                         attention_dim,\n",
    "#                         attention_hops,\n",
    "#                         num_dense_neurons,\n",
    "#                         dropout_rate=0.25)\n",
    "\n",
    "signal_model = get_normed_multilabel_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        attention_hops,\n",
    "                        num_dense_neurons,\n",
    "                        dropout_rate=0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "signal_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# parallel_attention_model = keras.utils.multi_gpu_model(signal_model, gpus=2)\n",
    "# parallel_attention_model.compile(loss=keras.losses.binary_crossentropy,\n",
    "#               optimizer=keras.optimizers.Adam(),\n",
    "#               metrics=['acc'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "signal_model.fit([x_train, x_rc_train], y_train,\n",
    "          batch_size=100,\n",
    "          epochs=10,\n",
    "          verbose=1,\n",
    "          validation_data=([x_test, x_rc_test], y_test))\n",
    "score = signal_model.evaluate([x_test, x_rc_test], y_test, verbose=0)\n",
    "print('Test loss:', score[0])\n",
    "print('Test accuracy:', score[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# # serialize model to JSON\n",
    "# model_json = signal_model.to_json()\n",
    "# with open(\"signal_model.json\", \"w\") as json_file:\n",
    "#     json_file.write(model_json)\n",
    "# # serialize weights to HDF5\n",
    "# signal_model.save_weights(\"signal_model.h5\")\n",
    "# print(\"Saved model to disk\")\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Neural Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "negative_seqRecords = list(SeqIO.parse('./background.fasta', 'fasta'))\n",
    "negative_seq = [str(x.seq[:200]) for x in negative_seqRecords]\n",
    "negative_rc_seq = [str(x.seq[:200]) for x in negative_seqRecords]\n",
    "\n",
    "negative_sequence_arrays = convert_sequences_to_array(negative_seq)\n",
    "negative_sequence_arrays = np.array(negative_sequence_arrays)\n",
    "\n",
    "negative_sequence_rc_arrays = convert_sequences_to_array(negative_rc_seq)\n",
    "negative_sequence_rc_arrays = np.array(negative_sequence_rc_arrays)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "positive_seqRecords = list(SeqIO.parse('./merged_atac_peaks_filtered_resized.fasta', 'fasta'))\n",
    "positive_seqRecords = [x for x in positive_seqRecords]\n",
    "fasta_seq = [str(x.seq[:200]) for x in positive_seqRecords]\n",
    "\n",
    "fasta_rc_seq = [str(x[:200].reverse_complement().seq) for x in positive_seqRecords ]\n",
    "\n",
    "seq_ids = [x.name for x in positive_seqRecords]\n",
    "\n",
    "sequence_arrays = convert_sequences_to_array(fasta_seq)\n",
    "sequence_arrays = np.array(sequence_arrays)\n",
    "\n",
    "sequence_rc_arrays = convert_sequences_to_array(fasta_rc_seq)\n",
    "sequence_rc_arrays = np.array(sequence_rc_arrays)\n",
    "\n",
    "index_seqArray_dict = dict(zip(seq_ids, zip(sequence_arrays, sequence_rc_arrays)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "summary_frame = pd.read_csv('./group_atac_summary.tsv' , sep='\\t', low_memory=False)\n",
    "summary_frame = summary_frame.fillna('0')\n",
    "for col in summary_frame.columns[5:]:\n",
    "    floatValues = []\n",
    "    for val in summary_frame[col].values.astype(str):\n",
    "        if ',' in val:\n",
    "            maxVal = np.mean([float(x) for x in val.split(',')])\n",
    "            floatValues.append(maxVal)\n",
    "        else:\n",
    "            floatValues.append(float(val))\n",
    "    summary_frame[col] = floatValues\n",
    "# summary_frame.index = summary_frame['ID'].values\n",
    "summary_frame.index = summary_frame['chr'] + ':' + (summary_frame['start'] - 1).astype(str) + '-' + summary_frame['end'].astype(str)\n",
    "\n",
    "# remove peaks in unknown/random chromosomes\n",
    "summary_frame = summary_frame[~summary_frame['chr'].str.contains('random')]\n",
    "summary_frame = summary_frame[~summary_frame['chr'].str.contains('Un')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_normed_regression_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        attention_hops,\n",
    "                        num_dense_neurons,\n",
    "                        dropout_rate=0.25,\n",
    "                               ):\n",
    "    input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "    input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "    ### find motifs ###\n",
    "    convolution_layer = Conv1D(filters=num_motifs, \n",
    "        kernel_size=motif_size,\n",
    "        activation='relu',\n",
    "        input_shape=(total_seq_length,4),\n",
    "        name='convolution_layer',\n",
    "        padding = 'same'\n",
    "        )\n",
    "    forward_motif_scores = convolution_layer(input_fwd)\n",
    "    reverse_motif_scores = convolution_layer(input_rev)\n",
    "    print('forward_motif_scores', forward_motif_scores.get_shape())\n",
    "\n",
    "    ### crop motif scores to avoid parts of sequence where motif score is computed in only one direction ###\n",
    "    to_crop = int((total_seq_length - seq_size)/2)\n",
    "    crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "        name='crop_layer')\n",
    "    cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "    cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "    print('cropped_fwd_scores', cropped_fwd_scores.get_shape())\n",
    "\n",
    "    ### normalize motif scores ###\n",
    "    motif_score_norm_layer = BatchNormalization(name='motif_score_norm_layer', axis=2)\n",
    "    normed_cropped_fwd_scores = motif_score_norm_layer(cropped_fwd_scores)\n",
    "    normed_cropped_rev_scores = motif_score_norm_layer(cropped_rev_scores)\n",
    "    print('normed_cropped_fwd_scores', normed_cropped_fwd_scores.shape)\n",
    "    \n",
    "    ### flip motif scores ###\n",
    "    flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "        output_shape=(seq_size, num_motifs),\n",
    "        name='flip_layer')\n",
    "    flipped_rev_scores = flip_layer(normed_cropped_fwd_scores)\n",
    "    print('flipped_rev_scores', flipped_rev_scores.get_shape())\n",
    "\n",
    "    ### concatenate motif scores ###\n",
    "    concatenate_layer = keras.layers.Concatenate(axis=2, name='concatenate_layer')\n",
    "    concatenated_motif_scores = concatenate_layer([normed_cropped_fwd_scores, flipped_rev_scores])\n",
    "    print('concatenated_motif_scores', concatenated_motif_scores.get_shape())\n",
    "\n",
    "    ### pool across length of sequence ###\n",
    "    sequence_pooling_layer = MaxPool1D(pool_size=adjacent_bp_pool_size, \n",
    "        strides=adjacent_bp_pool_size,\n",
    "        name='sequence_pooling_layer')\n",
    "    pooled_scores = sequence_pooling_layer(concatenated_motif_scores)\n",
    "    print('pooled_scores', pooled_scores.get_shape())\n",
    "    \n",
    "    ### bidirectional LSTM ###\n",
    "    forward_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "        return_sequences=True,\n",
    "        input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "        name = 'forward_lstm_layer'\n",
    "        )\n",
    "    forward_hidden_states = forward_lstm_layer(pooled_scores)\n",
    "    print('forward_hidden_states', forward_hidden_states.get_shape())\n",
    "\n",
    "    reverse_lstm_layer = LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "        return_sequences=True,\n",
    "        input_shape = (int(seq_size/adjacent_bp_pool_size), 2*num_motifs),\n",
    "        name = 'reverse_lstm_layer',\n",
    "        go_backwards=True,\n",
    "        )\n",
    "    reverse_hidden_states = reverse_lstm_layer(pooled_scores)\n",
    "    print('reverse_hidden_states', reverse_hidden_states.get_shape())    \n",
    "    \n",
    "    ### concatenate lstm hidden states ###\n",
    "    lstm_concatenate_layer = Concatenate(axis=2)\n",
    "    bilstm_hidden_states = lstm_concatenate_layer([forward_hidden_states, reverse_hidden_states])\n",
    "    print('bilstm_hidden_states', bilstm_hidden_states.get_shape())\n",
    "    \n",
    "#     bilstm_layer = Bidirectional(\n",
    "#         LSTM(units=int(seq_size/adjacent_bp_pool_size),\n",
    "#              return_sequences=True)\n",
    "#     )\n",
    "#     bilstm_hidden_states = bilstm_layer(pooled_scores)\n",
    "    \n",
    "    ### normalize lstm states ###\n",
    "    lstm_norm_layer = BatchNormalization(name='lstm_norm_layer', axis=2)\n",
    "    normed_bilstm_hidden_states = lstm_norm_layer(bilstm_hidden_states)\n",
    "    \n",
    "    \n",
    "    ### attention tanh layer ###\n",
    "    attention_tanh_layer = Dense(attention_dim,\n",
    "        activation='tanh',\n",
    "        use_bias=False,\n",
    "        name = 'attention_tanh_layer')\n",
    "    attention_tanh_layer_out = attention_tanh_layer(normed_bilstm_hidden_states)\n",
    "    print('attention_tanh_layer_out', attention_tanh_layer_out.get_shape())\n",
    "\n",
    "    ### outer layer ###\n",
    "    attention_outer_layer = Dense(attention_hops,\n",
    "        activation='linear',\n",
    "        use_bias=False,\n",
    "        name = 'attention_outer_layer')\n",
    "    attention_outer_layer_out = attention_outer_layer(attention_tanh_layer_out)\n",
    "    print('attention_outer_layer_out', attention_outer_layer_out.get_shape())\n",
    "\n",
    "    ### apply softmax ###\n",
    "    softmax_layer = Softmax(axis=1, name='attention_softmax_layer')\n",
    "    attention_softmax_layer_out = softmax_layer(attention_outer_layer_out)\n",
    "    print('attention_softmax_layer_out', attention_softmax_layer_out.get_shape())\n",
    "\n",
    "    ### attend to hidden states ###\n",
    "    attending_layer = Dot(axes=(1,1),\n",
    "        name='attending_layer')\n",
    "\n",
    "    attended_states = attending_layer([attention_softmax_layer_out, normed_bilstm_hidden_states])\n",
    "    print('attended_states', attended_states.get_shape())\n",
    "            \n",
    "    ### fully connected layer ###\n",
    "    dense_layer = Dense(num_dense_neurons, \n",
    "        activation='relu', \n",
    "        name = 'dense_layer'\n",
    "        )\n",
    "\n",
    "    dense_output = dense_layer(attended_states)\n",
    "    print('dense_output', dense_output.shape)\n",
    "    \n",
    "    # drop out\n",
    "    drop_out = Dropout(dropout_rate,name='dense_dropout')(dense_output)\n",
    "    print('drop_out', drop_out.shape)\n",
    "    \n",
    "    # make prediction\n",
    "    flattened = Flatten(name='flatten')(drop_out)\n",
    "    print('flattened', flattened.shape)\n",
    "    \n",
    "    predictions = Dense(num_classes,\n",
    "                        name='predictions',\n",
    "                        activation = 'relu', \n",
    "                       )(flattened)\n",
    "    print('predictions', predictions.shape)\n",
    "    \n",
    "    # define and compile model\n",
    "    model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# normed_h3k27ac_meanTag_frame = pd.read_csv('./normed_h3k27ac_reads.tsv', sep='\\t', index_col=0)\n",
    "normed_h3k27ac_meanTag_frame = pd.read_csv('./h3k27ac_reads.tsv', sep='\\t', index_col=0)\n",
    "\n",
    "normed_h3k27ac_meanTag_frame = normed_h3k27ac_meanTag_frame[normed_h3k27ac_meanTag_frame.max(axis=1)>=0]\n",
    "\n",
    "target_indices = normed_h3k27ac_meanTag_frame.index.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "target_seq_arrays = np.array([index_seqArray_dict[x][0] for x in target_indices])\n",
    "target_seq_rc_arrays = np.array([index_seqArray_dict[x][1] for x in target_indices])\n",
    "\n",
    "seq_arrays = target_seq_arrays\n",
    "seq_rc_arrays = target_seq_rc_arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "labels = normed_h3k27ac_meanTag_frame.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_train, x_test, x_rc_train, x_rc_test, y_train, y_test = model_selection.train_test_split(seq_arrays, seq_rc_arrays, labels, test_size=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "54771"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(target_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /home/jenhan/anaconda3/lib/python3.6/site-packages/tensorflow/python/util/deprecation.py:497: calling conv1d (from tensorflow.python.ops.nn_ops) with data_format=NHWC is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "`NHWC` for data_format is deprecated, use `NWC` instead\n",
      "forward_motif_scores (?, 200, 50)\n",
      "cropped_fwd_scores (?, 150, 50)\n",
      "normed_cropped_fwd_scores (?, 150, 50)\n",
      "flipped_rev_scores (?, 150, 50)\n",
      "concatenated_motif_scores (?, 150, 100)\n",
      "pooled_scores (?, 30, 100)\n",
      "forward_hidden_states (?, ?, 30)\n",
      "reverse_hidden_states (?, ?, 30)\n",
      "bilstm_hidden_states (?, ?, 60)\n",
      "attention_tanh_layer_out (?, 30, 50)\n",
      "attention_outer_layer_out (?, 30, 1)\n",
      "attention_softmax_layer_out (?, 30, 1)\n",
      "attended_states (?, 1, 60)\n",
      "dense_output (?, 1, 50)\n",
      "drop_out (?, 1, 50)\n",
      "flattened (?, ?)\n",
      "predictions (?, 3)\n"
     ]
    }
   ],
   "source": [
    "num_classes = len(labels[0])\n",
    "total_seq_length = 200\n",
    "seq_size = 150\n",
    "num_motifs = 50\n",
    "motif_size = 20\n",
    "adjacent_bp_pool_size = 5\n",
    "attention_dim = 50\n",
    "attention_hops = 1\n",
    "num_dense_neurons = 50 \n",
    "\n",
    "\n",
    "signal_regression_model = get_normed_regression_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        attention_hops,\n",
    "                        num_dense_neurons,\n",
    "                        dropout_rate=0.75)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_fwd (InputLayer)          (None, 200, 4)       0                                            \n",
      "__________________________________________________________________________________________________\n",
      "convolution_layer (Conv1D)      (None, 200, 50)      4050        input_fwd[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "crop_layer (Cropping1D)         (None, 150, 50)      0           convolution_layer[0][0]          \n",
      "__________________________________________________________________________________________________\n",
      "motif_score_norm_layer (BatchNo (None, 150, 50)      200         crop_layer[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "flip_layer (Lambda)             (None, 150, 50)      0           motif_score_norm_layer[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_layer (Concatenate) (None, 150, 100)     0           motif_score_norm_layer[0][0]     \n",
      "                                                                 flip_layer[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "sequence_pooling_layer (MaxPool (None, 30, 100)      0           concatenate_layer[0][0]          \n",
      "__________________________________________________________________________________________________\n",
      "forward_lstm_layer (LSTM)       (None, 30, 30)       15720       sequence_pooling_layer[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "reverse_lstm_layer (LSTM)       (None, 30, 30)       15720       sequence_pooling_layer[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_1 (Concatenate)     (None, 30, 60)       0           forward_lstm_layer[0][0]         \n",
      "                                                                 reverse_lstm_layer[0][0]         \n",
      "__________________________________________________________________________________________________\n",
      "lstm_norm_layer (BatchNormaliza (None, 30, 60)       240         concatenate_1[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "attention_tanh_layer (Dense)    (None, 30, 50)       3000        lstm_norm_layer[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "attention_outer_layer (Dense)   (None, 30, 1)        50          attention_tanh_layer[0][0]       \n",
      "__________________________________________________________________________________________________\n",
      "attention_softmax_layer (Softma (None, 30, 1)        0           attention_outer_layer[0][0]      \n",
      "__________________________________________________________________________________________________\n",
      "attending_layer (Dot)           (None, 1, 60)        0           attention_softmax_layer[0][0]    \n",
      "                                                                 lstm_norm_layer[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "dense_layer (Dense)             (None, 1, 50)        3050        attending_layer[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "dense_dropout (Dropout)         (None, 1, 50)        0           dense_layer[0][0]                \n",
      "__________________________________________________________________________________________________\n",
      "flatten (Flatten)               (None, 50)           0           dense_dropout[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "predictions (Dense)             (None, 3)            153         flatten[0][0]                    \n",
      "==================================================================================================\n",
      "Total params: 42,183\n",
      "Trainable params: 41,963\n",
      "Non-trainable params: 220\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "keras.utils.plot_model(signal_regression_model, show_layer_names=False, to_file='regression_model.pdf')\n",
    "signal_regression_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "signal_regression_model.compile(loss=keras.losses.mean_squared_logarithmic_error,\n",
    "              optimizer=keras.optimizers.RMSprop(),\n",
    "              )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "parallel_attention_model = keras.utils.multi_gpu_model(signal_regression_model, gpus=2)\n",
    "parallel_attention_model.compile(loss=keras.losses.mean_squared_logarithmic_error,\n",
    "              optimizer=keras.optimizers.RMSprop())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 43816 samples, validate on 10955 samples\n",
      "Epoch 1/10\n",
      "43816/43816 [==============================] - 25s 581us/step - loss: 2.5484 - val_loss: 1.8419\n",
      "Epoch 2/10\n",
      "43816/43816 [==============================] - 15s 348us/step - loss: 1.8636 - val_loss: 1.9882\n",
      "Epoch 3/10\n",
      "43816/43816 [==============================] - 15s 346us/step - loss: 1.7879 - val_loss: 1.7017\n",
      "Epoch 4/10\n",
      "43816/43816 [==============================] - 15s 350us/step - loss: 1.7479 - val_loss: 1.6989\n",
      "Epoch 5/10\n",
      "43816/43816 [==============================] - 15s 351us/step - loss: 1.6827 - val_loss: 1.7323\n",
      "Epoch 6/10\n",
      "43816/43816 [==============================] - 15s 350us/step - loss: 1.6549 - val_loss: 1.6614\n",
      "Epoch 7/10\n",
      "43816/43816 [==============================] - 15s 351us/step - loss: 1.6383 - val_loss: 1.6855\n",
      "Epoch 8/10\n",
      "43816/43816 [==============================] - 15s 349us/step - loss: 1.6195 - val_loss: 1.7919\n",
      "Epoch 9/10\n",
      "43816/43816 [==============================] - 15s 350us/step - loss: 1.5977 - val_loss: 1.7352\n",
      "Epoch 10/10\n",
      "43816/43816 [==============================] - 15s 351us/step - loss: 1.5904 - val_loss: 1.6140\n"
     ]
    }
   ],
   "source": [
    "parallel_attention_model.fit([x_train, x_rc_train], y_train,\n",
    "          batch_size=200,\n",
    "          epochs=10,\n",
    "          verbose=1,\n",
    "          validation_data=([x_test, x_rc_test], y_test))\n",
    "score = signal_regression_model.evaluate([x_test, x_rc_test], y_test, verbose=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predictions = signal_regression_model.predict([seq_arrays, seq_rc_arrays])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predictions_frame = pd.DataFrame(predictions, \n",
    "                                 index = normed_h3k27ac_meanTag_frame.index.values, \n",
    "                                 columns=[x + ' predictions' for x in normed_h3k27ac_meanTag_frame.columns.values])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "il4-24h 0.3931757844046414 0.40827040896917205\n",
      "kla-1h 0.5128383276025223 0.527032632155776\n",
      "veh 0.41502417144727793 0.4376147095697875\n"
     ]
    }
   ],
   "source": [
    "for cond in predictions_frame.columns:\n",
    "    pearson, _  = scipy.stats.spearmanr(predictions_frame[cond],normed_h3k27ac_meanTag_frame[cond.replace(' predictions','')] )\n",
    "    log_pearson, _= scipy.stats.pearsonr(np.log2(predictions_frame[cond]+1),\n",
    "                     np.log2(normed_h3k27ac_meanTag_frame[cond.replace(' predictions','')]+1) )\n",
    "    print(cond.replace(' predictions',''), pearson, log_pearson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "test_predictions = signal_regression_model.predict([x_test, x_rc_test])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "il4-24h 0.34543139745022206 0.362525637963554\n",
      "kla-1h 0.47771161744012963 0.49453613436608584\n",
      "veh 0.36828097421186684 0.3911475028381343\n"
     ]
    }
   ],
   "source": [
    "conditions = normed_h3k27ac_meanTag_frame.columns\n",
    "for i in range(y_test.shape[1]):\n",
    "    pearson, _  = scipy.stats.spearmanr(test_predictions[:,i],\n",
    "                                      y_test[:,i])\n",
    "    log_pearson, _= scipy.stats.pearsonr(np.log2(test_predictions[:,i]+1),\n",
    "                                      np.log2(y_test[:,i]+1))\n",
    "\n",
    "    print(conditions[i], pearson, log_pearson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 855,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "large_indices = normed_h3k27ac_meanTag_frame[normed_h3k27ac_meanTag_frame.max(axis=1) >64].index.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 856,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "diff_seq_arrays = np.array([index_seqArray_dict[x][0] for x in large_indices])\n",
    "diff_seq_rc_arrays = np.array([index_seqArray_dict[x][1] for x in large_indices])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 857,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "large_labels = normed_h3k27ac_meanTag_frame.loc[large_indices,:].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 858,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "differential_predictions = signal_regression_model.predict([diff_seq_arrays, diff_seq_rc_arrays])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 859,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "il4-24h 0.08789134500858443 0.09397435451392305\n",
      "kla-1h 0.4166379672267237 0.4010839445884486\n",
      "veh 0.08716051303096252 0.09405555669399457\n"
     ]
    }
   ],
   "source": [
    "for i in range(y_test.shape[1]):\n",
    "    pearson, _  = scipy.stats.spearmanr(differential_predictions[:,i],\n",
    "                                      large_labels[:,i])\n",
    "    log_pearson, _= scipy.stats.pearsonr(np.log2(differential_predictions[:,i]+1),\n",
    "                                      np.log2(large_labels[:,i]+1))\n",
    "\n",
    "    print(conditions[i], pearson, log_pearson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 778,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predictions_frame.to_csv('./predictions.tsv', sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Attention Only Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_attentionOnly_regression_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        num_dense_neurons,\n",
    "                        dropout_rate=0.25,\n",
    "                               ):\n",
    "    input_fwd = Input(shape=(total_seq_length,4), name='input_fwd')\n",
    "    input_rev = Input(shape=(total_seq_length,4), name='input_rev')\n",
    "\n",
    "    ### find motifs ###\n",
    "    convolution_layer = Conv1D(filters=num_motifs, \n",
    "        kernel_size=motif_size,\n",
    "        activation='relu',\n",
    "        input_shape=(total_seq_length,4),\n",
    "        name='convolution_layer',\n",
    "        padding = 'same'\n",
    "        )\n",
    "    forward_motif_scores = convolution_layer(input_fwd)\n",
    "    reverse_motif_scores = convolution_layer(input_rev)\n",
    "    print('forward_motif_scores', forward_motif_scores.get_shape())\n",
    "\n",
    "    ### crop motif scores to avoid parts of sequence where motif score is computed in only one direction ###\n",
    "    to_crop = int((total_seq_length - seq_size)/2)\n",
    "    crop_layer = Cropping1D(cropping=(to_crop, to_crop), \n",
    "        name='crop_layer')\n",
    "    cropped_fwd_scores = crop_layer(forward_motif_scores)\n",
    "    cropped_rev_scores = crop_layer(reverse_motif_scores)\n",
    "    print('cropped_fwd_scores', cropped_fwd_scores.get_shape())\n",
    "\n",
    "    ### normalize motif scores ###\n",
    "    motif_score_norm_layer = BatchNormalization(name='motif_score_norm_layer', axis=2)\n",
    "    normed_cropped_fwd_scores = motif_score_norm_layer(cropped_fwd_scores)\n",
    "    normed_cropped_rev_scores = motif_score_norm_layer(cropped_rev_scores)\n",
    "    print('normed_cropped_fwd_scores', normed_cropped_fwd_scores.shape)\n",
    "    \n",
    "    ### flip motif scores ###\n",
    "    flip_layer = Lambda(lambda x: K.reverse(x,axes=0),\n",
    "        output_shape=(seq_size, num_motifs),\n",
    "        name='flip_layer')\n",
    "    flipped_rev_scores = flip_layer(normed_cropped_fwd_scores)\n",
    "    print('flipped_rev_scores', flipped_rev_scores.get_shape())\n",
    "\n",
    "    ### concatenate motif scores ###\n",
    "    concatenate_layer = keras.layers.Concatenate(axis=2, name='concatenate_layer')\n",
    "    concatenated_motif_scores = concatenate_layer([normed_cropped_fwd_scores, flipped_rev_scores])\n",
    "    print('concatenated_motif_scores', concatenated_motif_scores.get_shape())\n",
    "\n",
    "    ### pool across length of sequence ###\n",
    "    sequence_pooling_layer = MaxPool1D(pool_size=adjacent_bp_pool_size, \n",
    "        strides=adjacent_bp_pool_size,\n",
    "        name='sequence_pooling_layer')\n",
    "    pooled_scores = sequence_pooling_layer(concatenated_motif_scores)\n",
    "    print('pooled_scores', pooled_scores.get_shape())\n",
    "    \n",
    "    attention_dot_layer = Dot(axes=(2,2))\n",
    "    \n",
    "    dot_products = attention_dot_layer([pooled_scores, pooled_scores])\n",
    "    print('dot_products', dot_products.shape)\n",
    "    \n",
    "    scaling_layer = Lambda(lambda x: x/(int(dot_products.shape[1])**-2),\n",
    "        name='scaling_layer')\n",
    "    scaled_dot_products = scaling_layer(dot_products)\n",
    "    \n",
    "    ### apply softmax ###\n",
    "    softmax_layer = Softmax(axis=1, name='attention_softmax_layer')\n",
    "    attention_softmax_layer_out = softmax_layer(scaled_dot_products)\n",
    "    print('attention_softmax_layer_out', attention_softmax_layer_out.get_shape())\n",
    "    \n",
    "    \n",
    "    sum_layer = Lambda(lambda x: K.sum(x,axis=2), name='sum_layer')\n",
    "\n",
    "    attention = sum_layer(attention_softmax_layer_out)\n",
    "    print('attention', attention.shape)\n",
    "    \n",
    "    repeat_layer = RepeatVector(n=num_motifs*2)\n",
    "    repeated_attention = repeat_layer(attention)\n",
    "    \n",
    "    permute_layer = Permute((2,1))\n",
    "    permuted_repeated_attention = permute_layer(repeated_attention)\n",
    "    print('permuted_repeated_attention', permuted_repeated_attention.shape)\n",
    "\n",
    "    ### attend to hidden states ###\n",
    "    attending_layer = Multiply(name='attending_layer')\n",
    "\n",
    "    attended_states = attending_layer([permuted_repeated_attention, pooled_scores])\n",
    "    print('attended_states', attended_states.get_shape())\n",
    "            \n",
    "    ### fully connected layer ###\n",
    "    dense_layer = Dense(num_dense_neurons, \n",
    "        activation='relu', \n",
    "        name = 'dense_layer'\n",
    "        )\n",
    "\n",
    "    dense_output = dense_layer(attended_states)\n",
    "    print('dense_output', dense_output.shape)\n",
    "    \n",
    "    # drop out\n",
    "    drop_out = Dropout(dropout_rate,name='dense_dropout')(dense_output)\n",
    "    print('drop_out', drop_out.shape)\n",
    "    \n",
    "    # make prediction\n",
    "    flattened = Flatten(name='flatten')(drop_out)\n",
    "    print('flattened', flattened.shape)\n",
    "    \n",
    "    predictions = Dense(num_classes,\n",
    "                        name='predictions',\n",
    "                        activation = 'relu', \n",
    "                       )(flattened)\n",
    "    print('predictions', predictions.shape)\n",
    "    \n",
    "    # define and compile model\n",
    "    model = Model(inputs=[input_fwd, input_rev], outputs=predictions)\n",
    "\n",
    "\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "forward_motif_scores (?, 200, 25)\n",
      "cropped_fwd_scores (?, 150, 25)\n",
      "normed_cropped_fwd_scores (?, 150, 25)\n",
      "flipped_rev_scores (?, 150, 25)\n",
      "concatenated_motif_scores (?, 150, 50)\n",
      "pooled_scores (?, 30, 50)\n",
      "dot_products (?, 30, 30)\n",
      "attention_softmax_layer_out (?, 30, 30)\n",
      "attention (?, 30)\n",
      "permuted_repeated_attention (?, 30, 50)\n",
      "attended_states (?, 30, 50)\n",
      "dense_output (?, 30, 150)\n",
      "drop_out (?, 30, 150)\n",
      "flattened (?, ?)\n",
      "predictions (?, 3)\n"
     ]
    }
   ],
   "source": [
    "num_classes = len(labels[0])\n",
    "total_seq_length = 200\n",
    "seq_size = 150\n",
    "num_motifs = 25\n",
    "motif_size = 20\n",
    "adjacent_bp_pool_size = 5\n",
    "attention_dim = 50\n",
    "num_dense_neurons = 150 \n",
    "\n",
    "\n",
    "signal_attentionOnly_regression_model = get_attentionOnly_regression_model(total_seq_length,\n",
    "                        seq_size,\n",
    "                        num_motifs, \n",
    "                        motif_size,\n",
    "                        adjacent_bp_pool_size,\n",
    "                        attention_dim,\n",
    "                        num_dense_neurons,\n",
    "                        dropout_rate=0.75)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_fwd (InputLayer)          (None, 200, 4)       0                                            \n",
      "__________________________________________________________________________________________________\n",
      "convolution_layer (Conv1D)      (None, 200, 25)      2025        input_fwd[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "crop_layer (Cropping1D)         (None, 150, 25)      0           convolution_layer[0][0]          \n",
      "__________________________________________________________________________________________________\n",
      "motif_score_norm_layer (BatchNo (None, 150, 25)      100         crop_layer[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "flip_layer (Lambda)             (None, 150, 25)      0           motif_score_norm_layer[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_layer (Concatenate) (None, 150, 50)      0           motif_score_norm_layer[0][0]     \n",
      "                                                                 flip_layer[0][0]                 \n",
      "__________________________________________________________________________________________________\n",
      "sequence_pooling_layer (MaxPool (None, 30, 50)       0           concatenate_layer[0][0]          \n",
      "__________________________________________________________________________________________________\n",
      "dot_59 (Dot)                    (None, 30, 30)       0           sequence_pooling_layer[0][0]     \n",
      "                                                                 sequence_pooling_layer[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "scaling_layer (Lambda)          (None, 30, 30)       0           dot_59[0][0]                     \n",
      "__________________________________________________________________________________________________\n",
      "attention_softmax_layer (Softma (None, 30, 30)       0           scaling_layer[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "sum_layer (Lambda)              (None, 30)           0           attention_softmax_layer[0][0]    \n",
      "__________________________________________________________________________________________________\n",
      "repeat_vector_16 (RepeatVector) (None, 50, 30)       0           sum_layer[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "permute_16 (Permute)            (None, 30, 50)       0           repeat_vector_16[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "attending_layer (Multiply)      (None, 30, 50)       0           permute_16[0][0]                 \n",
      "                                                                 sequence_pooling_layer[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "dense_layer (Dense)             (None, 30, 150)      7650        attending_layer[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "dense_dropout (Dropout)         (None, 30, 150)      0           dense_layer[0][0]                \n",
      "__________________________________________________________________________________________________\n",
      "flatten (Flatten)               (None, 4500)         0           dense_dropout[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "predictions (Dense)             (None, 3)            13503       flatten[0][0]                    \n",
      "==================================================================================================\n",
      "Total params: 23,278\n",
      "Trainable params: 23,228\n",
      "Non-trainable params: 50\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "keras.utils.plot_model(signal_attentionOnly_regression_model, show_layer_names=False, to_file='regression_model.pdf')\n",
    "signal_attentionOnly_regression_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "signal_attentionOnly_regression_model.compile(loss=keras.losses.mean_squared_logarithmic_error,\n",
    "              optimizer=keras.optimizers.RMSprop(),\n",
    "              )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "parallel_attention_model = keras.utils.multi_gpu_model(signal_attentionOnly_regression_model, gpus=2)\n",
    "parallel_attention_model.compile(loss=keras.losses.mean_squared_logarithmic_error,\n",
    "              optimizer=keras.optimizers.RMSprop())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 43816 samples, validate on 10955 samples\n",
      "Epoch 1/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6504 - val_loss: 1.6876\n",
      "Epoch 2/10\n",
      "43816/43816 [==============================] - 2s 38us/step - loss: 1.6463 - val_loss: 1.7601\n",
      "Epoch 3/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6402 - val_loss: 1.6764\n",
      "Epoch 4/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6351 - val_loss: 1.8655\n",
      "Epoch 5/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6325 - val_loss: 1.7243\n",
      "Epoch 6/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6302 - val_loss: 1.7349\n",
      "Epoch 7/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6302 - val_loss: 1.7109\n",
      "Epoch 8/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6230 - val_loss: 1.6794\n",
      "Epoch 9/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6146 - val_loss: 1.6743\n",
      "Epoch 10/10\n",
      "43816/43816 [==============================] - 2s 37us/step - loss: 1.6130 - val_loss: 1.6909\n"
     ]
    }
   ],
   "source": [
    "parallel_attention_model.fit([x_train, x_rc_train], y_train,\n",
    "          batch_size=200,\n",
    "          epochs=10,\n",
    "          verbose=1,\n",
    "          validation_data=([x_test, x_rc_test], y_test))\n",
    "score = signal_regression_model.evaluate([x_test, x_rc_test], y_test, verbose=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}