Goal¶
- build a model to predict the exonuclease cuts given the sequence (fully seq->seq)
Design¶
Schema¶
- input: 200x4 (one-hot encoded DNA sequence)
- output: 200x2 (chip-exo counts, normalized/softmax across each channel)
Loss function¶
- Use multinomial log-likelihood
- this will control for the amount of reads at each position
Expected dificulties¶
- the netowork might be kind of correct but mis the signal for a single base.
- Surprizingly the network is able to predict the peaks very accurately so this is not a real problem
Architecture¶
- U-net like architecture with Transposed-Conv/upscale layers
Data split¶
- 60
- 20
- 20
TODO¶
- [x] get the data
- [x] split the data
- write the data function
- [x] try out the simplest possible architecture
- one conv layer
- batch norm?
- one de-conv layer (with large enough window size)
- try motif initialization for the two motifs
- [x] visualize the weights
- de-conv layer -> ala averaged signal over the motifs?
- [x] evaluate the accuracy on the test-set
- [x] visualize the prediction results for:
- most correct
- most wrong
- most uncertain even though it should be correct
- [ ] evaluate the performance for predicting high-count peaks
- how many can we predict correctly - with base-pair precision?
In [1]:
%env CUDA_VISIBLE_DEVICES=3
In [2]:
import os
In [3]:
os.environ['CUDA_VISIBLE_DEVICES']
Out[3]:
In [4]:
import tensorflow as tf
In [5]:
#tf.enable_eager_execution()
import keras
get the data and split it into train,test¶
In [16]:
from basepair.config import get_data_dir
import pyBigWig
import matplotlib.pyplot as plt
from concise.preprocessing import encodeDNA
from concise.utils.fasta import read_fasta
from pybedtools import BedTool
from tqdm import tqdm
import pandas as pd
import numpy as np
from basepair.data import get_sox2_data, seq_inp_exo_out
from basepair.math import softmax
In [10]:
keras.backend.get_session().list_devices()
Out[10]:
In [12]:
ddir = get_data_dir()
In [13]:
dfc = get_sox2_data()
In [15]:
train, test = seq_inp_exo_out()
In [17]:
train[0].shape
Out[17]:
In [18]:
train[1].shape
Out[18]:
In [19]:
test[0].shape
Out[19]:
In [20]:
test[1].shape
Out[20]:
In [21]:
train[1][2][:,0]
Out[21]:
In [22]:
y_batch = train[1][:64]
In [23]:
from basepair.losses import twochannel_multinomial_nll
try out the simplest possible architecture¶
- one conv layer
- batch norm?
- one de-conv layer (with large enough window size)
- try motif initialization for the two motifs
In [35]:
import concise.layers as cl
import concise.initializers as ci
import keras.layers as kl
from keras.optimizers import Adam
from keras.callbacks import EarlyStopping, ModelCheckpoint
from keras.models import Model, load_model
In [25]:
motifs = {"motif1": f"{ddir}/processed/chipnexus/motifs/sox2/homer_200bp/de-novo/motif1.motif",
"motif2": f"{ddir}/processed/chipnexus/motifs/sox2/homer_200bp/de-novo/motif2.motif"}
In [26]:
from basepair.motif.homer import load_motif, read_motif_hits
pwm_list = [load_motif(fname) for k,fname in motifs.items()]
In [27]:
for i,pwm in enumerate(pwm_list):
pwm.plotPWMInfo((5,1.5))
plt.title(f"motif{i+1}")
In [32]:
# TODO - build the reverse complement symmetry into the model
inp = cl.Input(shape=(201, 4))
filters = 32
first_conv = cl.Conv1D(filters, kernel_size=21, padding='same',
#kernel_initializer = ci.PSSMKernelInitializer(pwm_list, stddev=0, add_noise_before_Pwm2Pssm=False),
#bias_initializer = 'zeros',
activation='relu')(inp)
second_conv = kl.Conv1D(filters, kernel_size=3, padding='same', activation='relu', dilation_rate=2)(first_conv)
c_second_conv = kl.add([first_conv, second_conv])
third_conv = kl.Conv1D(filters, kernel_size=3, padding='same', activation='relu', dilation_rate=4)(c_second_conv)
c_third_conv = kl.add([first_conv, second_conv, third_conv])
fourth_conv = kl.Conv1D(filters, kernel_size=3, padding='same', activation='relu', dilation_rate=8)(c_third_conv)
c_fourth_conv = kl.add([first_conv, second_conv, third_conv, fourth_conv])
fifth_conv = kl.Conv1D(filters, kernel_size=3, padding='same', activation='relu', dilation_rate=16)(c_fourth_conv)
c_fifth_conv = kl.add([first_conv, second_conv, third_conv, fourth_conv, fifth_conv])
sixt_conv = kl.Conv1D(filters, kernel_size=3, padding='same', activation='relu', dilation_rate=32)(c_fifth_conv)
c_sixt_conv = kl.add([first_conv, second_conv, third_conv, fourth_conv, fifth_conv, sixt_conv])
seventh_conv = kl.Conv1D(filters, kernel_size=3, padding='same', activation='relu', dilation_rate=64)(c_sixt_conv)
combined_conv = kl.add([first_conv, second_conv, third_conv, fourth_conv, fifth_conv, sixt_conv, seventh_conv])
# Bottleneck layer
combined_conv = kl.Conv1D(3, kernel_size=1, padding='same', activation='relu')(combined_conv)
x = kl.Reshape((-1, 1, 3))(combined_conv)
x = kl.Conv2DTranspose(2, kernel_size=(25, 1), padding='same')(x)
#kl.Conv2DTranspose(32, kernel_size=(7, 1), padding='same', activation='relu'),
#kl.Conv2DTranspose(2, kernel_size=(3, 1), padding='same'),
out = kl.Reshape((-1, 2))(x)
model = Model(inp, out)
model.compile(Adam(lr=0.004), loss=twochannel_multinomial_nll)
Notes¶
Experiments¶
- resresnet_nconv=4_filters=256=_lr=0.004.h5 - 108238.6558
- looks pretty good
- resnest_allconnect_nconv=4_filters=256=_lr=0.004.h5 - 108588.4661
- resnest_allconnect_nconv=4_filters=64_lr=0.004.h5 - 108546.4132 - Amazing results
- resnest_allconnect_nconv=4_filters=32_lr=0.004.h5 108670.3445 - great looking as well
- resnest_allconnect_nconv=5_filters=32_lr=0.004.h5 107394.5353 - even better
- resnest_allconnect_nconv=6_filters=32_lr=0.004.h5 107187.9883 - even better
- resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=False.h5 109961.7913 - no dillation
- resnest_allconnect_nconv=7_filters=64_lr=0.004_dilated=True.h5 - 107460.7238 - great results
- resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True-2.h5 - 107225.4570 - amazing results
- resnest_allconnect_nconv=7_filters=16_lr=0.004_dilated=True.h5 - 107225.4570 - not as good as previous, but still cool
start messing around with output kernel size (before it was 25)
- resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True,out=13.h5 - 107437.4955 - qualitatively better looking
- resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True,out=31.h5 - 107032.3784 - not so reproducible
- resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True,out=21.h5 - 107032.3784 - similarly ok
- resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True,out=41.h5 - 107032.3784 - similarly ok
- resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True,out=25.h5 - 106892.1230 - qualitatively ok
- resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True,out=25,bottleneck=3.h5 - 109486.1274 - performance not as good, but hopefully more interepretable
Dillation helps
In [535]:
!mkdir -p {ddir}/processed/chipnexus/exp/models
In [34]:
# best valid so far: 108238.6558
ckp_file = f"{ddir}/processed/chipnexus/exp/models/resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True,out=25,bottleneck=3.h5"
model.fit(train[0], train[1],
batch_size=256,
epochs=200,
validation_split=0.2,
callbacks=[EarlyStopping(patience=5),
ModelCheckpoint(ckp_file, save_best_only=True)]
)
# get the best model
model = load_model(ckp_file, custom_objects={"twochannel_multinomial_nll": twochannel_multinomial_nll})
evaluate the accuracy on the test-set¶
visualize the prediction results for¶
- most correct
- most wrong
- most uncertain even though it should be correct
In [37]:
import basepair
In [38]:
model.evaluate(test[0], test[1])
Out[38]:
In [39]:
y_pred = softmax(model.predict(train[0]))
y_true = train[1]
np.allclose(y_pred.sum(axis=1),1)
Out[39]:
In [40]:
idx_list = pd.Series(np.arange(len(test[0]))).sample(10)
for idx in idx_list:
plt.figure(figsize=(10,2))
plt.subplot(121)
plt.plot(y_pred[idx,:,0], label='pos,m={}'.format(np.argmax(y_pred[idx,:,0])))
plt.plot(y_pred[idx,:,1], label='neg,m={}'.format(np.argmax(y_pred[idx,:,1])))
plt.legend();
plt.subplot(122)
plt.plot(y_true[idx,:,0], label='pos,m={}'.format(np.argmax(y_true[idx,:,0])))
plt.plot(y_true[idx,:,1], label='neg,m={}'.format(np.argmax(y_true[idx,:,1])))
plt.legend();
In [41]:
print("Pos")
print(np.argmax(y_pred[idx_list,:,0], axis=1))
print(np.argmax(y_true[idx_list,:,0], axis=1))
print("Neg")
print(np.argmax(y_pred[idx_list,:,1], axis=1))
print(np.argmax(y_true[idx_list,:,1], axis=1))
In [42]:
y_pred = softmax(model.predict(test[0]))
y_true = test[1]
np.allclose(y_pred.sum(axis=1),1)
Out[42]:
In [43]:
idx_list = pd.Series(np.arange(len(test[0]))).sample(10)
for idx in idx_list:
plt.figure(figsize=(10,2))
plt.subplot(121)
plt.plot(y_pred[idx,:,0], label='pos,m={}'.format(np.argmax(y_pred[idx,:,0])))
plt.plot(y_pred[idx,:,1], label='neg,m={}'.format(np.argmax(y_pred[idx,:,1])))
plt.legend();
plt.subplot(122)
plt.plot(y_true[idx,:,0], label='pos,m={}'.format(np.argmax(y_true[idx,:,0])))
plt.plot(y_true[idx,:,1], label='neg,m={}'.format(np.argmax(y_true[idx,:,1])))
plt.legend();
In [ ]:
# Inspect the intermediate activations
visualize the weights¶
- de-conv layer -> ala averaged signal over the motifs?
In [159]:
from concise.preprocessing.sequence import one_hot2string, DNA
In [315]:
one_hot2string(train[0][0][np.newaxis], DNA)[0]
Out[315]:
In [316]:
# Add a softmax layer to the model
In [319]:
w=model.layers[-2].get_weights()[0]
In [320]:
w.shape
Out[320]:
In [348]:
# Plot weights
w=model.layers[-2].get_weights()[0]
plt.figure(figsize=(20,3))
filters = w.shape[-1]
for i in range(filters):
plt.subplot(1,filters, i+1)
plt.plot(w[:,0,0,i], label='pos')
plt.plot(w[:,0,1,i], label='neg')
plt.title("Motif{}".format(i))
plt.legend();
In [349]:
from concise.utils.plot import seqlogo, seqlogo_fig
In [350]:
from concise.utils.pwm import pssm_array2pwm_array, _pwm2pwm_info
In [351]:
pssm_array2pwm_array
Out[351]:
In [352]:
pwm_list[0].plotPSSM(figsize=(5,4));
In [353]:
seqlogo_fig(model.layers[0].get_weights()[0], figsize=(10,4), ncol=2);
In [ ]:
# GC gets strongly enhanced -> GC bias?
In [339]:
model.layers[0].plot_weights();
In [340]:
seqlogo_fig(model.layers[0].get_weights()[0], figsize=(10,2), ncol=2);
In [55]:
l.plot_weights?
In [ ]:
model.layers[0].plot_weights