# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + [markdown] colab_type="text" id="view-in-github"
# <a href="https://colab.research.google.com/github/kundajelab/dragonn/blob/master/paper_supplement/PrimerTutorial%205%20-%20Functional%20variant%20characterization%20for%20non-coding%20SNPs%20within%20the%20SPI1%20motif.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + [markdown] colab_type="text" id="XIu-8nkCV4zm"
# # How to train your DragoNN tutorial 5: 
# ## Functional variant characterization for non-coding SNPs within the SPI1 motif 
#
# This tutorial is a supplement to the DragoNN manuscript. 
#
# This tutorial will take 2 - 3 hours if executed on a GPU.
#
# ## Outline<a name='outline'>
# <ol>
#     <li><a href=#1>Input data: SPI1 ChiP-seq and experimental bQTL data</a></li>
#     <li><a href=#2>Genomewide classification and regression labels for SPI1 TF ChiPseq</a></li>
#     <li><a href=#3>Optional: Download pre-generated models and test-set predictions</a></li>
#     <li><a href=#4>Genome-wide classification for SPI1</a></li>
#     <li><a href=#5>Genome-wide regression for SPI1</a></li> 
#     <li><a href=#6>Genome-wide interpretation of true positive predictions in SPI1, with DeepLIFT</a></li>
#     <li><a href=#7>Recovering bQTL effect sizes: Classification vs Regression</a></li>
#     <li><a href=#8>Model-predicted SNP effect sizes vs bQTL effect sizes</a></li>
#     <li><a href=#a>Kat's architecture: Classification Model</a></li>
#     <li><a href=#b>Kat's architecture: Regression Model</a></li> 
#     <li><a href=#9>Conclusions</a></li>    
#     <li><a href=#10>Save tutorial outputs</a></li>
# </ol>
# Github issues on the [dragonn repository](https://github.com/kundajelab/dragonn) with feedback, questions, and discussion are always welcome.
#

# + colab={} colab_type="code" id="72XlYRZBluGr"
# If you don't have bedtools installed in your environment (i.e. Google Colab), uncomment and run the command below 
# #!apt-get install bedtools
# #!pip install pybedtools

# + colab={} colab_type="code" id="FmftiCCDV4zo"
#uncomment the lines below if you are running this tutorial from Google Colab 
# #!pip install dragonn>=0.2.2

# + colab={} colab_type="code" id="fyLzeiF5V4zq"
# Making sure our results are reproducible
from numpy.random import seed
seed(1234)
from tensorflow import set_random_seed
set_random_seed(1234)

# + colab={} colab_type="code" id="8M6gdfuJV4zu"
#load dragonn tutorial utilities 
# %reload_ext autoreload
# %autoreload 2
# %matplotlib inline
import warnings
warnings.filterwarnings('ignore')


from dragonn.tutorial_utils import *

# + [markdown] colab_type="text" id="djDLAi21V4zy"
# ## Input data <a name='1'>
# <a href=#outline>Home</a>
#
# This tutorial uses the same in vivo SPI1 transcription factor CHiP-seq dataset that was used in [Tutorial 4](https://colab.research.google.com/github/kundajelab/dragonn/blob/keras_2.2_tensorflow_1.6_purekeras/paper_supplement/PrimerTutorial%204%20-%20Interpreting%20predictive%20sequence%20features%20in%20in-vivo%20TF%20binding%20events.ipynb). Our goal is to compare predicted variant effect sizes from classification and regression models against experimental bQTL data. The bQTL data in this way serves as a "gold-standard" validation that in silico mutagenesis on the deep learning inputs leads to correct variant effect size prediction.  We  will use bQTL data  that has been intersected with SPI1 CISBP genome motif annotations. 

# + colab={"base_uri": "https://localhost:8080/", "height": 413} colab_type="code" id="O707uf21V4zy" outputId="f5ca2dbc-9594-4a62-aa67-97190945d622"
# SPI1, optimal IDR thresholded peaks, Myers lab, hg19
# https://www.encodeproject.org/experiments/ENCSR000BGQ/
# !wget -O SPI1.narrowPeak.gz http://mitra.stanford.edu/kundaje/projects/dragonn/dragonn_gm12878_pipeline/spi1_ENCSR000BGQ/cromwell-executions/chip/bb0c3c5a-3889-43fe-a218-05851cecc74a/call-reproducibility_idr/execution/optimal_peak.regionPeak.gz

#Fold change bigWig track for the SPI1 dataset: 
# !wget -O SPI1.pooled.fc.bigWig http://mitra.stanford.edu/kundaje/projects/dragonn/dragonn_gm12878_pipeline/spi1_ENCSR000BGQ/cromwell-executions/chip/bb0c3c5a-3889-43fe-a218-05851cecc74a/call-macs2_pooled/execution/ENCFF000OBU.Rep1.merged.nodup.pooled_x_ENCFF000OCW.Control.Rep1.merged.nodup.fc.signal.bigwig
    
## Download the hg19 chromsizes file (We only use chroms 1 -22, X, Y for training)
# !wget http://mitra.stanford.edu/kundaje/projects/dragonn/hg19.chrom.sizes
    
## Download the hg19 fasta reference genome (and corresponding .fai index)
# !wget http://mitra.stanford.edu/kundaje/projects/dragonn/hg19.genome.fa.gz
# !wget http://mitra.stanford.edu/kundaje/projects/dragonn/hg19.genome.fa.gz.fai 
# !wget http://mitra.stanford.edu/kundaje/projects/dragonn/hg19.genome.fa.gz.gzi 

    




# + colab={"base_uri": "https://localhost:8080/", "height": 215} colab_type="code" id="-YwnqCV-V4z2" outputId="85791b16-8647-45ff-cd8b-cf281d618350"
# Download bQTL experimental data for SPI1 loci 
# !wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.bQTLs.txt.gz

# + [markdown] colab_type="text" id="sp9mi-6_V4z4"
# ## Generating genome-wide classification and regression labels <a name='2'>
# <a href=#outline>Home</a>

# + [markdown] colab_type="text" id="Zmt5OJP_V4z5"
# We will use the *genomewide_labels* function from the  [seqdataloader](https://github.com/kundajelab/seqdataloader) package to generate positive and negative labels for the TF-ChIPseq peaks across the genome. We will treat each sample as a task for the model and compare the performance of the model on SPI1 task in the single-tasked and multi-tasked setting.

# + colab={} colab_type="code" id="SLGpH2rOV4z6"
from seqdataloader import * 

# + colab={"base_uri": "https://localhost:8080/", "height": 35} colab_type="code" id="u2wpRugxV4z7" outputId="c695d444-b6d0-408b-f294-6afcdc6b1033"
## seqdataloader accepts an input file, which we call SPI1.tasks.tsv, with task names in column 1, corresponding
## peak files in column 2, and the signal track in column 3. In this tutorial, the task file will have a single task entry for the SPI1 TF CHiP-seq
with open("SPI1.task.tsv",'w') as f: 
    f.write("SPI1\tSPI1.narrowPeak.gz\tSPI1.pooled.fc.bigWig\n")
f.close() 
# !cat SPI1.task.tsv

# + [markdown] colab_type="text" id="-pqz2oVGV4z_"
# With the parameter configuration below, seqdataloader splits the genome into 1kb regions, with a stride of 50. Each 1kb region is centered at a 200 bp bin, with a left flank of 400 bases and a right flank of 400 bases. 
#
# * In the classification case, each 200 bp bin is labeled as positive if a narrowPeak summit overlaps with it. The bin is labeled negative if there is no overlap with the narrowPeak. 
# * In the regression case, the asinh(mean coverage) in the 200 bp bin is computed. 
#

# + [markdown] colab_type="text" id="e--f8QWuV4z_"
# **Note**: The label generation may take 10 - 15 minutes to complete. If you prefer not to wait, you can download the 
# pre-generated classification and regression labels for the training, validation, and test sets by uncommenting the code below: 

# + colab={"base_uri": "https://localhost:8080/", "height": 1205} colab_type="code" id="LMG9IzPnV40A" outputId="aca1e05d-dc62-416e-b22f-b68a53aaf3f7"
## Classification labels 
# ! wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.train.classification.hdf5
# ! wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.valid.classification.hdf5
# ! wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.test.classification.hdf5

## Regression labels 
# ! wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.train.regression.hdf5
# ! wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.valid.regression.hdf5
# ! wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.test.regression.hdf5


# + [markdown] colab_type="text" id="_lleRdAaV40B"
# If you prefer to generate the labels from scratch, execute the two code cell below: 

# + colab={} colab_type="code" id="erh3B4hIV40D"
#  Generate genome-wide classification labels 

#1) Training set: all chromosomes with the exception of 1,2, and 19 in our training set. Also, the dataset does not
# include chromosome Y, so we exclude it as well. 

train_set_params={
    'task_list':"SPI1.task.tsv",
    'outf':"SPI1.train.classification.hdf5",
    'output_type':'hdf5',
    'chrom_sizes':'hg19.chrom.sizes',
    'chroms_to_exclude':['chr1','chr2','chr19','chrY'],
    'bin_stride':50,
    'left_flank':400,
    'right_flank':400,
    'bin_size':200,
    'threads':4,
    'subthreads':4,
    'allow_ambiguous':False,
    'labeling_approach':'peak_summit_in_bin_classification'
    }
genomewide_labels(train_set_params)

#2) Validation set: Chromosome 1
valid_set_params={'task_list':"SPI1.task.tsv",
    'outf':"SPI1.valid.classification.hdf5",
    'output_type':'hdf5',
    'chrom_sizes':'hg19.chrom.sizes',
    'chroms_to_keep':'chr1',
    'bin_stride':50,
    'left_flank':400,
    'right_flank':400,
    'bin_size':200,
    'threads':1,
    'subthreads':4,
    'allow_ambiguous':False,
    'labeling_approach':'peak_summit_in_bin_classification'
    }
genomewide_labels(valid_set_params)

#3) Test set: Chromosomes 2, 19 
test_set_params={
    'task_list':"SPI1.task.tsv",
    'outf':"SPI1.test.classification.hdf5",
    'output_type':'hdf5',
    'chrom_sizes':'hg19.chrom.sizes',
    'chroms_to_keep':['chr2','chr19'],
    'bin_stride':50,
    'left_flank':400,
    'right_flank':400,
    'bin_size':200,
    'threads':2,
    'subthreads':4,
    'allow_ambiguous':False,
    'labeling_approach':'peak_summit_in_bin_classification'
    }
genomewide_labels(test_set_params)


# + colab={} colab_type="code" id="PBqbysuAV40G"
# Generate regression labels genome-wide 

#1) Training set: all chromosomes with the exception of 1,2, and 19 in our training set 

train_set_params={
    'task_list':"SPI1.task.tsv",
    'outf':"SPI1.train.regression.hdf5",
    'output_type':'hdf5',
    'chrom_sizes':'hg19.chrom.sizes',
    'chroms_to_exclude':['chr1','chr2','chr19','chrY'],
    'bin_stride':50,
    'left_flank':400,
    'right_flank':400,
    'bin_size':200,
    'threads':4,
    'subthreads':4,
    'allow_ambiguous':False,
    'labeling_approach':'all_genome_bins_regression'
    }
genomewide_labels(train_set_params)

#2) Validation set: Chromosome 1
valid_set_params={'task_list':"SPI1.task.tsv",
    'outf':"SPI1.valid.regression.hdf5",
    'output_type':'hdf5',
    'chrom_sizes':'hg19.chrom.sizes',
    'chroms_to_keep':'chr1',
    'bin_stride':50,
    'left_flank':400,
    'right_flank':400,
    'bin_size':200,
    'threads':1,
    'subthreads':4,
    'allow_ambiguous':False,
    'labeling_approach':'all_genome_bins_regression'
    }
genomewide_labels(valid_set_params)

#3) Test set: Chromosomes 2, 19 
test_set_params={
    'task_list':"SPI1.task.tsv",
    'outf':"SPI1.test.regression.hdf5",
    'output_type':'hdf5',
    'chrom_sizes':'hg19.chrom.sizes',
    'chroms_to_keep':['chr2','chr19'],
    'bin_stride':50,
    'left_flank':400,
    'right_flank':400,
    'bin_size':200,
    'threads':2,
    'subthreads':4,
    'allow_ambiguous':False,
    'labeling_approach':'all_genome_bins_regression'
    }
genomewide_labels(test_set_params)


# + [markdown] colab_type="text" id="x1m9HsgXV40J"
# Let's examine the files that were generated: 

# + colab={"base_uri": "https://localhost:8080/", "height": 390} colab_type="code" id="Q0SgKnZtV40J" outputId="5e5501e6-61e3-44e1-f6ca-94e70f3ff63f"
#The code generates bed file outputs with a label of 1 or 0 for each 1kb
# genome bin for each task. Note that the bins are shifted with a stride of 50.
pd.read_hdf("SPI1.train.classification.hdf5",start=1000000,stop=1000010)

# + colab={"base_uri": "https://localhost:8080/", "height": 390} colab_type="code" id="weH35tmhV40N" outputId="a8aea293-f560-400f-cb2f-2d0085625934"
pd.read_hdf("SPI1.train.regression.hdf5",start=1000000,stop=1000010)

# + [markdown] colab_type="text" id="FKHBBFpRV40Q"
# ## Optional: Download pre-generated models and test-set predictions <a name='3'>
# <a href=#outline>Home</a>
#
# Next, we will train classification and regression models to predict TF CHiP-seq peaks for SPI1. If you want to skip straight to model interpretation and bQTL analysis, you can download the pre-trained models by uncommenting the 
# block of code below. 

# + colab={} colab_type="code" id="CqyIROINV40R"
from keras.models import load_model 

## Download classification model 
# #! wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.classification.model.hdf5
spi1_classification_model=load_model("SPI1.kat.classification.model.hdf5")

## Download regression model 
# #! wget http://mitra.stanford.edu/kundaje/projects/dragonn/SPI1.regression.model.hdf5
spi1_regression_model=load_model("SPI1.kat.regression.model.hdf5")


## Get test set classification model and regression model predictions 
#import h5py
#test_set_predictions=h5py.File("SPI1.test.predictions.hdf5")
#spi1_test_classification_predictions=test_set_predictions['classification'].value 
#spi1_test_regression_predictions=test_set_predictions['regression'].value 

# + [markdown] colab_type="text" id="FZcwz5AmV40U"
# ## Genome-wide classification model <a name='4'>
# <a href=#outline>Home</a>
#

# + colab={} colab_type="code" id="ZeBBukYGV40V"
#To prepare for model training, we import the necessary functions and submodules from keras
from keras.models import Sequential
from keras.layers.core import Dropout, Reshape, Dense, Activation, Flatten
from keras.layers.convolutional import Conv2D, MaxPooling2D
from keras.optimizers import Adadelta, SGD, RMSprop;
import keras.losses;
from keras.constraints import maxnorm;
from keras.layers.normalization import BatchNormalization
from keras.regularizers import l1, l2
from keras.callbacks import EarlyStopping, History
from keras import backend as K 
K.set_image_data_format('channels_last')

# + colab={} colab_type="code" id="CXtJhYf2V40Z"
from concise.metrics import tpr, tnr, fpr, fnr, precision, f1
def initialize_classification_model(ntasks=1):
    #Define the model architecture in keras (regularized, 3-layer convolution model followed by 1 dense layer)
    model=Sequential() 
    
    model.add(Conv2D(filters=15,kernel_size=(1,10),input_shape=(1,1000,4)))
    model.add(Activation('relu'))
    model.add(Dropout(0.2))
    model.add(MaxPooling2D(pool_size=(1,35)))

    model.add(Conv2D(filters=15,kernel_size=(1,10)))
    model.add(Activation('relu'))
    model.add(Dropout(0.2))

    model.add(Conv2D(filters=15,kernel_size=(1,10)))
    model.add(Activation('relu'))
    model.add(Dropout(0.2))

    model.add(Flatten())
    model.add(Dense(ntasks))
    model.add(Activation("sigmoid"))

    ##compile the model, specifying the Adam optimizer, and binary cross-entropy loss. 
    model.compile(optimizer='adam',loss='binary_crossentropy',
                  metrics=[tpr,
                           tnr,
                           fpr,
                           fnr,
                           precision,
                           f1])
    return model


# + [markdown] colab_type="text" id="XY-g6ik9V40c"
# We create generators for the training and validation data: 

# + colab={"base_uri": "https://localhost:8080/", "height": 35} colab_type="code" id="5PrSMLLiV40d" outputId="d7ae2438-6e8e-4fc0-c1b3-6bde39667604"
#create the generators, upsample positives to ensure they constitute 30% of each batch 
from dragonn.generators import * 
spi1_train_classification_gen=DataGenerator("SPI1.train.classification.hdf5","hg19.genome.fa.gz",upsample_ratio=0.3, batch_size=256)
spi1_valid_classification_gen=DataGenerator("SPI1.valid.classification.hdf5","hg19.genome.fa.gz",upsample_ratio=0.3, batch_size=256)


# + colab={} colab_type="code" id="NLRfBhebV40e"
#Train the SPI1 classification model 
spi1_classification_model=initialize_classification_model()

## use the keras fit_generator function to train the model with early stopping after 3 epochs 
history_classification=spi1_classification_model.fit_generator(spi1_train_classification_gen,
                                                  validation_data=spi1_valid_classification_gen,
                                                  steps_per_epoch=10000,
                                                  validation_steps=5000,
                                                  epochs=150,
                                                  verbose=1,
                                                  use_multiprocessing=True,
                                                  workers=40,
                                                  max_queue_size=100,
                                                  callbacks=[EarlyStopping(patience=3,restore_best_weights=True),History()])

# + colab={} colab_type="code" id="DN2bFz_aV40h" outputId="323bc266-fc71-4e30-e277-b2db3c3e9a0a"
## Plot the learning curves for SPI1  
from dragonn.tutorial_utils import plot_learning_curve
plot_learning_curve(history_classification)

# + [markdown] colab_type="text" id="YoP336y_V40k"
# We now measure how well the model performed by calculating performance metrics on the test splits across the whole genome. 

# + colab={} colab_type="code" id="ZF6VNZH3V40k" outputId="101ca0a8-db49-44cc-c133-e799b070e1e0"
from dragonn.generators import * 
spi1_test_classification_gen=DataGenerator("SPI1.test.classification.hdf5",
                                       "hg19.genome.fa.gz",
                                         upsample=False,
                                         add_revcomp=False,
                                         batch_size=1000,
                                         tasks=['SPI1'])
spi1_test_classification_predictions=spi1_classification_model.predict_generator(spi1_test_classification_gen,
                                                               max_queue_size=5000, 
                                                               workers=40, 
                                                               use_multiprocessing=True, 
                                                               verbose=1)
spi1_test_classification_truth=spi1_test_classification_gen.data

# + colab={} colab_type="code" id="1QGwsaLsV40n" outputId="829e46ff-12f4-42c9-c18c-2d9c66a22e97"
spi1_test_classification_predictions.shape

# + colab={} colab_type="code" id="9LAq8j_MV40r" outputId="d1b31200-f1a6-493b-ad99-6c0c62b65645"
spi1_test_classification_truth.shape

# + colab={} colab_type="code" id="X83Mt6VXV40t" outputId="b7b4fe79-812d-463a-f2f1-611122a6fa0e"
## Generate a ClassificationResult object to print performance metrics on held-out test set 
from dragonn.metrics import ClassificationResult
print(ClassificationResult(spi1_test_classification_truth.values.astype(bool),spi1_test_classification_predictions))

# + colab={} colab_type="code" id="jBvkHg7zV40y"
#save the models 
spi1_classification_model.save("SPI1.classification.model.hdf5")
# -



# + [markdown] colab_type="text" id="b0sibNElV403"
# ## Genome-wide regression model <a name='5'>
# <a href=#outline>Home</a>

# + colab={} colab_type="code" id="nM0-1aa9V404"
def initialize_regression_model(ntasks=1):
    #Define the model architecture in keras (regularized, 3-layer convolution model followed by 1 dense layer)
    model=Sequential() 
    
    model.add(Conv2D(filters=15,kernel_size=(1,10),input_shape=(1,1000,4)))
    model.add(Activation('relu'))
    model.add(Dropout(0.2))
    model.add(MaxPooling2D(pool_size=(1,35)))

    model.add(Conv2D(filters=10,kernel_size=(1,10)))
    model.add(Activation('relu'))
    model.add(Dropout(0.2))

    model.add(Conv2D(filters=5,kernel_size=(1,10)))
    model.add(Activation('relu'))
    model.add(Dropout(0.2))

    model.add(Flatten())
    model.add(Dense(ntasks))

    ##compile the model, specifying the Adam optimizer, and binary cross-entropy loss. 
    model.compile(optimizer='adam',loss='mse')
    return model


# + colab={} colab_type="code" id="zO-uKK8GV407" outputId="04c1a178-b8ee-4b14-c944-328c17ecec9c"
#we want to determine a threshold for upsampling the non-zero bins in a given batch 
# extract 5 million datapoints from the training data and observe the distribution of non-zero signal values  
sample=pd.read_hdf("SPI1.train.regression.hdf5",start=0,stop=5000000)
nonzero_sample=sample[sample.max(axis=1)>0]
print(nonzero_sample.shape)
nonzero_sample.hist(bins=100)

# + [markdown] colab_type="text" id="F5Raqyz3V40_"
# This suggests that 0.1 is a reasonable threshold for upsampling signal bins in regression 

# + colab={} colab_type="code" id="JSjGIlYTV41B"
#create the generators, no upsampling of positives is used for regression. 
from dragonn.generators import * 
spi1_train_regression_gen=DataGenerator("SPI1.train.regression.hdf5","hg19.genome.fa.gz",upsample_ratio=0.3,upsample_thresh=0.01)
spi1_valid_regression_gen=DataGenerator("SPI1.valid.regression.hdf5","hg19.genome.fa.gz",upsample_ratio=0.3,upsample_thresh=0.01)


# + colab={} colab_type="code" id="P3pSs4TYV41E" outputId="bb50b654-2336-4b6d-d16c-a8a15672f379"
#Train the SPI1 regression model 
spi1_regression_model=initialize_regression_model()

## use the keras fit_generator function to train the model with early stopping after 3 epochs 
history_regression=spi1_regression_model.fit_generator(spi1_train_regression_gen,
                                                  validation_data=spi1_valid_regression_gen,
                                                  steps_per_epoch=10000,
                                                  validation_steps=5000,
                                                  epochs=150,
                                                  verbose=1,
                                                  use_multiprocessing=True,
                                                  workers=40,
                                                  max_queue_size=100,
                                                  callbacks=[EarlyStopping(patience=3,restore_best_weights=True),History()])

# + colab={} colab_type="code" id="EFxneNJqV41H" outputId="472b64dc-f85d-49e6-8bd0-4d023374e04a"
plot_learning_curve(history_regression)

# + colab={} colab_type="code" id="_o3I4gN7V41K" outputId="da95db30-a6ad-4043-86c6-ed587a44277b"
from dragonn.generators import * 
spi1_test_regression_gen=DataGenerator("SPI1.test.regression.hdf5",
                                       "hg19.genome.fa.gz",
                                         upsample=False,
                                         add_revcomp=False,
                                         batch_size=1000,
                                         tasks=['SPI1'])
spi1_test_regression_predictions=spi1_regression_model.predict_generator(spi1_test_regression_gen,
                                                               max_queue_size=5000, 
                                                               workers=40, 
                                                               use_multiprocessing=True, 
                                                               verbose=1)
spi1_test_regression_truth=spi1_test_regression_gen.data

# + colab={} colab_type="code" id="QMbEZo8XGd2m"
## find the indices of the non-zero coverage bins 
nonzero_bins=spi1_test_regression_truth.max(axis=1)>0

# + colab={} colab_type="code" id="-byNK4qMV41N" outputId="aa18b7e6-6ee2-439d-ca0b-d68d2a7aa8a1"
#Calculate spearman and pearson correlation between truth labels and predictions 
from scipy.stats import pearsonr, spearmanr
corr_pearson=pearsonr(spi1_test_regression_truth,spi1_test_regression_predictions)
corr_spearman=spearmanr(spi1_test_regression_truth,spi1_test_regression_predictions)
print("Pearson correlation on test set:"+str(corr_pearson))
print("Spearman correlation on test set:"+str(corr_spearman))

# + colab={} colab_type="code" id="YipLiXRFGd2t" outputId="f7e2d334-30b4-4661-d8bb-23832c8c5763"
# Calculate the spearman and pearson correlation, restricted to non-zero bins 
corr_pearson_nonzero_bins=pearsonr(spi1_test_regression_truth[nonzero_bins],spi1_test_regression_predictions[nonzero_bins])
corr_spearman_nonzero_bins=spearmanr(spi1_test_regression_truth[nonzero_bins],spi1_test_regression_predictions[nonzero_bins])
print("Pearson correlation on test set:"+str(corr_pearson_nonzero_bins))
print("Spearman correlation on test set:"+str(corr_spearman_nonzero_bins))

# + colab={} colab_type="code" id="kVEfZ1d0V41O"
#There is some overfitting, let's save this model and see if we can do better 
spi1_regression_model.save("SPI1.regression.model.hdf5")

# + colab={} colab_type="code" id="UGW2FQrgGd2y" outputId="e3069d0a-65d3-4c8c-be0d-eef77ed83581"
spi1_test_regression_truth.values[0:10].squeeze()

# + colab={} colab_type="code" id="SKa34qjbGd21"
test_df=pd.DataFrame({"Observed":list(spi1_test_regression_truth.values.squeeze()),
                     "Predicted":list(spi1_test_regression_predictions.squeeze())})

# + colab={} colab_type="code" id="L1tUvhU5Gd27"
test_df_nonzero=pd.DataFrame({"Observed":list(spi1_test_regression_truth[nonzero_bins].values.squeeze()),
                     "Predicted":list(spi1_test_regression_predictions[nonzero_bins].squeeze())})

# + colab={} colab_type="code" id="Xqvik10lGd29" outputId="4d35b751-13ce-4555-f798-916765d4e981"
import plotnine 
from plotnine import * 
print((ggplot(test_df,aes(x="Observed",y="Predicted"))
 +geom_bin2d(bins=100)
 +theme_bw()
 +xlab("Observed asinh(mean coverage in FC bigWig")
 +ylab("Model prediction")
 +ggtitle("SPI1 regression model test set prediction")))

print((ggplot(test_df_nonzero,aes(x="Observed",y="Predicted"))
 +geom_bin2d(bins=100)
 +theme_bw()
 +xlab("Observed asinh(mean coverage in FC bigWig")
 +ylab("Model prediction")
 +ggtitle("SPI1 regression model test set prediction: bins with nonzero coverage")))

# + colab={} colab_type="code" id="TOBvin88Gd3A" outputId="5fc16ecd-16d2-47f0-a39f-e91244e4d9cb"
# Plot observed vs predicted regression values 
plt.scatter(spi1_test_regression_truth, spi1_test_regression_predictions, alpha=0.01)
plt.xlabel("Observed asinh(mean coverage in FC bigWig)")
plt.ylabel("Model prediction")
plt.title("SPI1 regression model test set prediction")
plt.show()

# + colab={} colab_type="code" id="WCgR4DmGGd3C" outputId="82246537-1073-48b1-81f7-1dc8f2a4b178"
# Plot observed vs predicted regression values for the nonzero bins 
plt.scatter(spi1_test_regression_truth, spi1_test_regression_predictions, alpha=0.01)
plt.xlabel("Observed asinh(mean coverage in FC bigWig) for bins ")
plt.ylabel("Model prediction")
plt.title("SPI1 regression model test set prediction: bins with nonzero coverage")
plt.show()

# + [markdown] colab_type="text" id="48jGNSBtV41R"
# ## Genome-wide interpretation of true positive predictions in SPI1, with DeepLIFT <a name='6'>
# <a href=#outline>Home</a>
#
# ### Classification Model

# + colab={} colab_type="code" id="RWumi0mQV41S"
#get the true positive predictions with a threshold of 0.9 (i.e. high confidence true positive predictions)
spi1_test_classification_truth_bool=spi1_test_classification_truth.values.astype(bool)
true_pos_spi1=spi1_test_classification_truth[spi1_test_classification_truth_bool*spi1_test_classification_predictions >0.9]
true_pos_spi1.head

# + colab={} colab_type="code" id="ecJa0a2HV41U" outputId="bb860e8e-f083-41e3-ed79-e66bb8aec2ac"
true_pos_spi1.shape

# + colab={} colab_type="code" id="UH5WvDmXV41W" outputId="c549e164-1087-4f2b-93bd-7f0677075090"
from dragonn.utils import one_hot_from_bed
deep_lift_input_spi1=one_hot_from_bed([i for i in true_pos_spi1.index],"hg19.genome.fa.gz")
deep_lift_input_spi1.shape

# + colab={} colab_type="code" id="41dInIg-V41Y"
from dragonn.tutorial_utils import deeplift 

# + colab={} colab_type="code" id="hViDtFpGV41a"
deep_lift_scores_spi1=deeplift(spi1_classification_model,deep_lift_input_spi1)

# + colab={} colab_type="code" id="mFkS3UwpV41c" outputId="953cdbc2-d7c0-4a17-d6d1-227acb47c325"
deep_lift_scores_spi1.shape

# + [markdown] colab_type="text" id="skCZQ8j0V41e"
# Let's plot a few of the DeepLIFT tracks and see if the model successfully learned SPI1:

# + colab={} colab_type="code" id="Srt-dPa_V41e"
from dragonn.tutorial_utils import  plot_seq_importance

# + colab={} colab_type="code" id="nhcSoq48V41h" outputId="fd5de43b-d65f-4abb-d998-7115c5ad1c1e"
plot_seq_importance(deep_lift_scores_spi1[0],deep_lift_input_spi1[0])

# + colab={} colab_type="code" id="eyejvCeWV41j" outputId="9985d83c-d05f-4979-f472-bb352cce2f80"
plot_seq_importance(deep_lift_scores_spi1[1],deep_lift_input_spi1[1])

# + colab={} colab_type="code" id="CyuioSW_V41o" outputId="75f83e18-bf46-475d-c43b-22ab8cdc56d9"
plot_seq_importance(deep_lift_scores_spi1[2],deep_lift_input_spi1[2])

# + [markdown] colab_type="text" id="AfkuGJTGV41r"
# Let's zoom in to the center of one sequence so that it is easier to distinguish the motif: 

# + colab={} colab_type="code" id="YCr3vtPeV41s" outputId="a04a10b2-ff18-44cb-f918-005872cb3b59"
plot_seq_importance(deep_lift_scores_spi1[2].squeeze()[550:650],deep_lift_input_spi1[2].squeeze()[550:650])

# + [markdown] colab_type="text" id="tiw5YPInV41v"
# If we query the sequence "CACTTCCCCT" in the [TomTom](http://meme-suite.org/tools/tomtom) software from the MEME suite, we find that the motif is a good match for SPIB: 
# <img src="https://github.com/kundajelab/dragonn/blob/master/paper_supplement/tutorial_images/SPI1.Tut4.png?raw=1" alt="SPI12TomTom" width="400"/>
#

# + [markdown] colab_type="text" id="kEUxg-r1V41x"
# ### Regression model 

# + colab={} colab_type="code" id="bEJudwWXV41y" outputId="f2395f38-be1f-4b76-95a3-91cc85efa471"
#Sanity-check that the model is learning the SPI1 motif by running DeepLIFT on True Positives with high confidence (>0.9)
#get the true positive predictions 
true_pos=spi1_test_regression_truth[(spi1_test_regression_truth.values*spi1_test_regression_predictions)>2]
true_pos.shape

# + colab={} colab_type="code" id="z5dpuvzLV413" outputId="c0db98aa-22bf-4e44-f15f-0a75df7c9239"
deep_lift_input=one_hot_from_bed([i for i in true_pos.index],"hg19.genome.fa.gz")
deep_lift_input.shape

# + colab={} colab_type="code" id="D-3kc90ZV416" outputId="514a9b2b-0f73-45ac-fe0c-b695bb182dbb"
help(deeplift)

# + colab={} colab_type="code" id="YeuZFMZqV41-"
deep_lift_scores_spi1=deeplift(spi1_regression_model,deep_lift_input_spi1,target_layer_idx=-1)

# + colab={} colab_type="code" id="8qigmzDOV41-" outputId="1ef6a324-0b04-4000-a658-87c75cb0c50d"
plot_seq_importance(deep_lift_scores_spi1[0],deep_lift_input_spi1[0])

# + colab={} colab_type="code" id="EPHe9I8VV42A" outputId="b7ee3839-2d4e-466b-e2f7-ede00db8a16b"
plot_seq_importance(deep_lift_scores_spi1[1],deep_lift_input_spi1[1])

# + colab={} colab_type="code" id="bogEKZN2V42C" outputId="6465c01d-b842-4217-9766-aca7a14797b2"
plot_seq_importance(deep_lift_scores_spi1[2],deep_lift_input_spi1[2])

# + colab={} colab_type="code" id="Ck63kAsAV42F" outputId="5353b1fb-b881-4515-f931-96f49f196045"
plot_seq_importance(deep_lift_scores_spi1[2].squeeze()[550:650],deep_lift_input_spi1[2].squeeze()[550:650])

# + [markdown] colab_type="text" id="PWqxtR6NV42I"
# The motif learned by the regression model matches the canonical SPI1 motif, though the deepLIFT tracks are noisier compared to those for the classification model. 
#

# + [markdown] colab_type="text" id="-PGA_k3RV42J"
# ## Recovering bQTL effect sizes: Classification vs Regression <a name='7'>
# <a href=#outline>Home</a>

# + colab={} colab_type="code" id="GROAoPZDV42J"
from dragonn.generators import * 
bqtl_ref_gen=BQTLGenerator("SPI1.bQTLs.txt.gz","hg19.genome.fa.gz","POSTallele")
bqtl_alt_gen=BQTLGenerator("SPI1.bQTLs.txt.gz","hg19.genome.fa.gz","ALTallele")

# + colab={} colab_type="code" id="GJq0Ic_8V42L" outputId="f27e9d3e-87a2-479f-90ca-320bf3066fc0"
bqtl_ref_classification_predictions=spi1_classification_model.predict_generator(bqtl_ref_gen,
                                                               max_queue_size=5000, 
                                                               workers=40, 
                                                               use_multiprocessing=True, 
                                                               verbose=1)

# + colab={} colab_type="code" id="xlPJpmmLV42L" outputId="bf8b8ff1-0268-49af-ce7d-9d8875e2ae3f"
bqtl_alt_classification_predictions=spi1_classification_model.predict_generator(bqtl_alt_gen,
                                                               max_queue_size=5000, 
                                                               workers=40, 
                                                               use_multiprocessing=True, 
                                                               verbose=1)
bqtl_ref_classification_truth=bqtl_ref_gen.data['pvalue']

# + colab={} colab_type="code" id="MGXJaSvPV42N" outputId="32d58541-15ab-4c82-90d5-b3118ed18520"
print(bqtl_ref_classification_predictions.shape)
print(bqtl_alt_classification_predictions.shape)
print(bqtl_ref_classification_truth.shape)

# + colab={} colab_type="code" id="1NpOM03tV42P" outputId="9e58fdff-c60a-463a-fa43-2a6979d96064"
bqtl_ref_regression_predictions=spi1_regression_model.predict_generator(bqtl_ref_gen,
                                                               max_queue_size=5000, 
                                                               workers=40, 
                                                               use_multiprocessing=True, 
                                                               verbose=1)
bqtl_alt_regression_predictions=spi1_regression_model.predict_generator(bqtl_alt_gen,
                                                               max_queue_size=5000, 
                                                               workers=40, 
                                                               use_multiprocessing=True, 
                                                               verbose=1)


# + colab={} colab_type="code" id="hi_Pr6jcV42Q" outputId="62710f23-734d-4ea4-b9a9-e1c28948d597"
plt.scatter(bqtl_ref_classification_predictions, bqtl_alt_classification_predictions, alpha=0.01)
plt.xlabel("Ref")
plt.ylabel("Alt")
plt.title("BQTL Classification Model Predictions")
plt.show()


# + colab={} colab_type="code" id="tg-ZRo1tV42R" outputId="9cb78eaf-ed5f-4fe4-ca02-6f927b07bcca"
plt.scatter(bqtl_ref_regression_predictions, bqtl_alt_regression_predictions, alpha=0.01)
plt.xlabel("Ref")
plt.ylabel("Alt")
plt.title("BQTL Regression Model Predictions")
plt.show()

# + [markdown] colab_type="text" id="yxG07_SzV42T"
# ## Model-predicted SNP effect sizes vs bQTL effect sizes <a name='8'>
# <a href=#outline>Home</a>

# + colab={} colab_type="code" id="YdFUgn60V42T"
logpval=np.log10(bqtl_ref_classification_truth.values)
delta=bqtl_alt_classification_predictions-bqtl_ref_classification_predictions


# + [markdown] colab_type="text" id="WdQsnAM7Gd4B"
# ## Kat's Model Architecture (Classification)<a name='a'>
# <a href=#outline>Home</a>

# + colab={} colab_type="code" id="7yIyOg-AGd4B"
from concise.metrics import tpr, tnr, fpr, fnr, precision, f1
from keras.constraints import max_norm

def initialize_kat_classification_model(ntasks=1):
    #Define the model architecture in keras (regularized, 3-layer convolution model followed by 1 dense layer)
    model=Sequential() 
    
    model.add(Conv2D(filters=50,kernel_size=(1,15),padding="same", kernel_constraint=max_norm(7.0,axis=-1),input_shape=(1,1000,4)))
    model.add(BatchNormalization(axis=-1))
    model.add(Activation('relu'))

    model.add(Conv2D(filters=50,kernel_size=(1,15),padding="same"))
    model.add(BatchNormalization(axis=-1))
    model.add(Activation('relu'))

    model.add(Conv2D(filters=50,kernel_size=(1,13),padding="same"))
    model.add(BatchNormalization(axis=-1))
    model.add(Activation('relu'))
    
    model.add(MaxPooling2D(pool_size=(1,40)))
    
    model.add(Flatten())
    model.add(Dense(50))
    model.add(BatchNormalization(axis=-1))
    model.add(Activation('relu'))
    model.add(Dropout(0.2))
    
    model.add(Dense(ntasks))
    model.add(Activation("sigmoid"))

    ##compile the model, specifying the Adam optimizer, and binary cross-entropy loss. 
    model.compile(optimizer='adam', loss='binary_crossentropy',
                  metrics=[tpr,
                           tnr,
                           fpr,
                           fnr,
                           precision,
                           f1])
    return model


# + colab={} colab_type="code" id="EG64R1eNGd4C"
#create the generators, upsample positives to ensure they constitute 30% of each batch 
from dragonn.generators import * 
spi1_train_classification_gen=DataGenerator("SPI1.train.classification.hdf5","hg19.genome.fa.gz",upsample_ratio=0.3, batch_size=256)
spi1_valid_classification_gen=DataGenerator("SPI1.valid.classification.hdf5","hg19.genome.fa.gz",upsample_ratio=0.3, batch_size=256)


# + colab={} colab_type="code" id="3SeWRHimGd4D" outputId="7516c239-52e8-4c19-dafc-ee32c264808c"
#Train the SPI1 classification model 
spi1_kat_classification_model=initialize_kat_classification_model()

## use the keras fit_generator function to train the model with early stopping after 3 epochs 
history_kat_classification=spi1_kat_classification_model.fit_generator(spi1_train_classification_gen,
                                                  validation_data=spi1_valid_classification_gen,
                                                  steps_per_epoch=10000,
                                                  validation_steps=5000,
                                                  epochs=150,
                                                  verbose=1,
                                                  use_multiprocessing=True,
                                                  workers=40,
                                                  max_queue_size=100,
                                                  callbacks=[EarlyStopping(patience=3,restore_best_weights=True),History()])

# + colab={} colab_type="code" id="a3D2uhXwGd4E" outputId="b5450d6e-447b-4a4a-e5c0-19890b4949bc"
## Plot the learning curves for SPI1  
from dragonn.tutorial_utils import plot_learning_curve
plot_learning_curve(history_kat_classification)

# + colab={} colab_type="code" id="Jq6tLmcTGd4G" outputId="9757df93-a9e1-45c9-85f9-8cc5c51fd926"
from dragonn.generators import * 
spi1_test_classification_gen=DataGenerator("SPI1.test.classification.hdf5",
                                       "hg19.genome.fa.gz",
                                         upsample=False,
                                         add_revcomp=False,
                                         batch_size=1000,
                                         tasks=['SPI1'])
spi1_test_classification_predictions=spi1_kat_classification_model.predict_generator(spi1_test_classification_gen,
                                                               max_queue_size=5000, 
                                                               workers=40, 
                                                               use_multiprocessing=True, 
                                                               verbose=1)
spi1_test_classification_truth=spi1_test_classification_gen.data

# + colab={} colab_type="code" id="-fWPY4DoGd4I" outputId="ab351e3b-7f55-465f-8f74-a5e8b7981d18"
## Generate a ClassificationResult object to print performance metrics on held-out test set 
from dragonn.metrics import ClassificationResult
print(ClassificationResult(spi1_test_classification_truth.values.astype(bool),spi1_test_classification_predictions))


# + [markdown] colab_type="text" id="dPWznGgeGd4K"
# ## Kat's Model Architecture (Regression)<a name='b'>
# <a href=#outline>Home</a>

# + colab={} colab_type="code" id="2Nt_b4_BGd4K"
def initialize_kat_regression_model(ntasks=1):
    #Define the model architecture in keras (regularized, 3-layer convolution model followed by 1 dense layer)
    model=Sequential() 
    
    model.add(Conv2D(filters=50,kernel_size=(1,15),padding="same", kernel_constraint=max_norm(7.0,axis=-1),input_shape=(1,1000,4)))
    model.add(BatchNormalization(axis=-1))
    model.add(Activation('relu'))

    model.add(Conv2D(filters=50,kernel_size=(1,15),padding="same"))
    model.add(BatchNormalization(axis=-1))
    model.add(Activation('relu'))

    model.add(Conv2D(filters=50,kernel_size=(1,13),padding="same"))
    model.add(BatchNormalization(axis=-1))
    model.add(Activation('relu'))
    
    model.add(MaxPooling2D(pool_size=(1,40)))
    
    model.add(Flatten())
    model.add(Dense(50))
    model.add(BatchNormalization(axis=-1))
    model.add(Activation('relu'))
    model.add(Dropout(0.2))
    
    model.add(Dense(ntasks))


    ##compile the model, specifying the Adam optimizer, and binary cross-entropy loss. 
    model.compile(optimizer='adam',loss='mse')
    return model

# + colab={} colab_type="code" id="Wm8P3ABcGd4L"
#create the generators, no upsampling of positives is used for regression. 
from dragonn.generators import * 
spi1_train_regression_gen=DataGenerator("SPI1.train.regression.hdf5","hg19.genome.fa.gz",upsample_ratio=0.3,upsample_thresh=0.01)
spi1_valid_regression_gen=DataGenerator("SPI1.valid.regression.hdf5","hg19.genome.fa.gz",upsample_ratio=0.3,upsample_thresh=0.01)


# + colab={} colab_type="code" id="O2xnIWHXGd4M" outputId="716c5c00-0b5b-417d-aa33-72cd483c7546"
#Train the SPI1 regression model 
spi1_kat_regression_model=initialize_kat_regression_model()

## use the keras fit_generator function to train the model with early stopping after 3 epochs 
history_kat_regression=spi1_kat_regression_model.fit_generator(spi1_train_regression_gen,
                                                  validation_data=spi1_valid_regression_gen,
                                                  steps_per_epoch=10000,
                                                  validation_steps=5000,
                                                  epochs=150,
                                                  verbose=1,
                                                  use_multiprocessing=True,
                                                  workers=40,
                                                  max_queue_size=100,
                                                  callbacks=[EarlyStopping(patience=3,restore_best_weights=True),History()])

# + colab={} colab_type="code" id="S2sQWdgtGd4O" outputId="3db4acbc-8de5-4d24-96da-3c4c41eac26a"
plot_learning_curve(history_kat_regression)

# + colab={} colab_type="code" id="JvSYpEFZGd4P" outputId="7fcf68b2-7f0e-4191-ae45-71948da1cd3f"
from dragonn.generators import * 
spi1_test_regression_gen=DataGenerator("SPI1.test.regression.hdf5",
                                       "hg19.genome.fa.gz",
                                         upsample=False,
                                         add_revcomp=False,
                                         batch_size=1000,
                                         tasks=['SPI1'])
spi1_test_regression_predictions=spi1_kat_regression_model.predict_generator(spi1_test_regression_gen,
                                                               max_queue_size=5000, 
                                                               workers=40, 
                                                               use_multiprocessing=True, 
                                                               verbose=1)
spi1_test_regression_truth=spi1_test_regression_gen.data

# + colab={} colab_type="code" id="2vIyYceoGd4R"
## find the indices of the non-zero coverage bins 
nonzero_bins=spi1_test_regression_truth.max(axis=1)>0

# + colab={} colab_type="code" id="pJ8NbtgWGd4T" outputId="8a0c8f6c-512f-419e-dd49-8002895cdeb7"
#Calculate spearman and pearson correlation between truth labels and predictions 
from scipy.stats import pearsonr, spearmanr
corr_pearson=pearsonr(spi1_test_regression_truth,spi1_test_regression_predictions)
corr_spearman=spearmanr(spi1_test_regression_truth,spi1_test_regression_predictions)
print("Pearson correlation on test set:"+str(corr_pearson))
print("Spearman correlation on test set:"+str(corr_spearman))

# + colab={} colab_type="code" id="QpnVA7VsGd4U" outputId="94a59775-862e-48e4-f036-d7d3fb4d7654"
# Calculate the spearman and pearson correlation, restricted to non-zero bins 
corr_pearson_nonzero_bins=pearsonr(spi1_test_regression_truth[nonzero_bins],spi1_test_regression_predictions[nonzero_bins])
corr_spearman_nonzero_bins=spearmanr(spi1_test_regression_truth[nonzero_bins],spi1_test_regression_predictions[nonzero_bins])
print("Pearson correlation on test set:"+str(corr_pearson_nonzero_bins))
print("Spearman correlation on test set:"+str(corr_spearman_nonzero_bins))

# + colab={} colab_type="code" id="KCaEkJICGd4V"
test_df=pd.DataFrame({"Observed":list(spi1_test_regression_truth.values.squeeze()),
                     "Predicted":list(spi1_test_regression_predictions.squeeze())})

# + colab={} colab_type="code" id="9kktZbRCGd4W"
test_df_nonzero=pd.DataFrame({"Observed":list(spi1_test_regression_truth[nonzero_bins].values.squeeze()),
                     "Predicted":list(spi1_test_regression_predictions[nonzero_bins].squeeze())})

# + colab={} colab_type="code" id="af-4fxu0Gd4X" outputId="bc675856-7680-4a00-d016-a4e0fe3e6d55"
import plotnine 
from plotnine import * 
print((ggplot(test_df,aes(x="Observed",y="Predicted"))
 +geom_bin2d(bins=100)
 +theme_bw()
 +xlab("Observed asinh(mean coverage in FC bigWig")
 +ylab("Model prediction")
 +ggtitle("SPI1 regression model test set prediction")))

print((ggplot(test_df_nonzero,aes(x="Observed",y="Predicted"))
 +geom_bin2d(bins=100)
 +theme_bw()
 +xlab("Observed asinh(mean coverage in FC bigWig")
 +ylab("Model prediction")
 +ggtitle("SPI1 regression model test set prediction: bins with nonzero coverage")))

# + [markdown] colab_type="text" id="cZnRq3SzGd4Z"
# ## Kat's Model DeepLIFT profiles (Classification)

# + colab={} colab_type="code" id="UdWiImsYGd4Z"
spi1_test_classification_truth_bool=spi1_test_classification_truth.values.astype(bool)
true_pos_spi1=spi1_test_classification_truth[spi1_test_classification_truth_bool*spi1_test_classification_predictions >0.9]

# + colab={} colab_type="code" id="6iSq-J-iGd4a" outputId="578fdbcd-8cd1-4d5d-c6e2-4d28234b6610"
from dragonn.utils import one_hot_from_bed
deep_lift_input_spi1=one_hot_from_bed([i for i in true_pos_spi1.index],"hg19.genome.fa.gz")
deep_lift_input_spi1.shape

# + colab={} colab_type="code" id="XRIkb8oFGd4c"
from dragonn.tutorial_utils import deeplift, plot_seq_importance
deep_lift_scores_spi1=deeplift(spi1_kat_classification_model,deep_lift_input_spi1)

# + colab={} colab_type="code" id="sCAY5G5JGd4g" outputId="580aca0a-9a35-4cc8-84bc-15e9bf399d56"
plot_seq_importance(deep_lift_scores_spi1[0],deep_lift_input_spi1[0])
plot_seq_importance(deep_lift_scores_spi1[1],deep_lift_input_spi1[1])
plot_seq_importance(deep_lift_scores_spi1[2],deep_lift_input_spi1[2])

# + colab={} colab_type="code" id="fefAyEp4Gd4h" outputId="dca31e18-7f19-48bd-e76d-86a98e6cb221"
plot_seq_importance(deep_lift_scores_spi1[2].squeeze()[400:500],deep_lift_input_spi1[2].squeeze()[400:500])

# + [markdown] colab_type="text" id="NaizmyYbGd4i"
# If we query the sequence "GTTTCACTTCTGCAAA" in the [TomTom](http://meme-suite.org/tools/tomtom) software from the MEME suite, we find that the motif is a good match (p=3.55e-03) for SPIB: 
# <img src="https://github.com/kundajelab/dragonn/blob/master/paper_supplement/tutorial_images/SPIB.Kat.png?raw=1" alt="SPI12TomTom" width="400"/>

# + [markdown] colab_type="text" id="RFlubWfKGd4i"
# ## Kat's Model DeepLIFT profiles (Regression)

# + colab={} colab_type="code" id="5DdIebdVGd4j" outputId="450f7622-0471-49af-fcb4-cd442ec04c57"
#Sanity-check that the model is learning the SPI1 motif by running DeepLIFT on True Positives with high confidence (>0.9)
#get the true positive predictions 
true_pos=spi1_test_regression_truth[(spi1_test_regression_truth.values*spi1_test_regression_predictions)>4]
true_pos.shape

# + colab={} colab_type="code" id="My8t2AvgGd4k" outputId="65d55c2a-48ca-4878-f554-5dcf00e33c9f"
deep_lift_input=one_hot_from_bed([i for i in true_pos.index],"hg19.genome.fa.gz")
deep_lift_input.shape

# + colab={} colab_type="code" id="PzURc6r9Gd4l"
deep_lift_scores_spi1=deeplift(spi1_regression_model,deep_lift_input_spi1,target_layer_idx=-1)

# + colab={} colab_type="code" id="tqrN9urUGd4n" outputId="709f3bef-fe1b-40d3-b775-d20740a84deb"
plot_seq_importance(deep_lift_scores_spi1[0],deep_lift_input_spi1[0])
plot_seq_importance(deep_lift_scores_spi1[1],deep_lift_input_spi1[1])
plot_seq_importance(deep_lift_scores_spi1[2],deep_lift_input_spi1[2])

# + colab={} colab_type="code" id="leSwvZWFGd4n" outputId="848dd89f-eaa3-4872-ae04-3301ac0fff7d"
plot_seq_importance(deep_lift_scores_spi1[2].squeeze()[400:500],deep_lift_input_spi1[2].squeeze()[400:500])

# + [markdown] colab_type="text" id="QMBWHvnZV42V"
# ## Conclusions <a name='9'>
# <a href=#outline>Home</a>

# + [markdown] colab_type="text" id="PFyKUGQnV42X"
# ## Save tutorial outputs <a name='10'>
# <a href=#outline>Home</a>
#
# We save the models and test set predictions generated in this tutorial to an hdf5 file so that they can be loaded more readily in the future. 

# + colab={} colab_type="code" id="WBNQNhg3V42X"
#save the models 
#spi1_kat_classification_model.save("SPI1.kat.classification.model.hdf5")
#spi1_kat_regression_model.save("SPI1.kat.regression.model.hdf5")

#spi1_classification_model.save("SPI1.classification.model.hdf5")
#spi1_regression_model.save("SPI1.regression.model.hdf5")
#save the test predictions 
import h5py 
test_set_predictions=h5py.File("SPI1.test.kat.predictions.hdf5",'w')
test_set_predictions.create_dataset("classification",data=spi1_test_classification_predictions)
test_set_predictions.create_dataset("regression",data=spi1_test_regression_predictions)
test_set_predictions.close() 


# + colab={} colab_type="code" id="LOPcHHY5V42a"


StarcoderdataJupyterScriptsDedupFiltered