Goal¶
- make a model that will also predict the absolute counts of the regions
TODO¶
- [x] histogram of the counts
- what's the right distribution
- make a q-q plot with poisson, negative bionomial
- what's the right distribution
- [x] train a model predicting also the counts
- focus only on sox2
- make the correlation plot between these two
Modeling pooled counts¶
- [x] implement the data function which returns the pooled signal
- [x] allow the model function to also work on this dataset
Next steps¶
- [ ] run the same thing for the multi-task model
- [ ] setup the sacred experiment for it
- run the evaluation overnight
- run 3 runs for each model. Select the best valid loss
- report the test performance
- run 3 runs for each model. Select the best valid loss
- run the evaluation overnight
- [ ] model the counts using the negative binomial distribution
In [6]:
import pandas as pd
import numpy as np
from pybedtools import BedTool
from basepair.config import get_data_dir, create_tf_session
from tqdm import tqdm
from concise.preprocessing import encodeDNA
from basepair.datasets import get_sox2_data
from basepair.plots import plot_loss
import pyBigWig
from basepair.math import softmax
import matplotlib.pyplot as plt
from basepair import samplers
ddir = get_data_dir()
In [91]:
create_tf_session(1)
Out[91]:
In [101]:
from basepair.datasets import *
from basepair.models import *
from basepair.plots import *
In [5]:
train, valid, test = seq_inp_exo_out() # default sox2 dataste
In [97]:
import keras.layers as kl
from keras.optimizers import Adam
from keras.models import Model
import keras.backend as K
from concise.utils.helper import get_from_module
from basepair.losses import twochannel_multinomial_nll
from basepair.layers import SpatialLifetimeSparsity
from basepair.models import seq_mutlitask
In [98]:
import keras.layers as kl
from keras.optimizers import Adam
from keras.callbacks import EarlyStopping, ModelCheckpoint, History
from keras.models import Model, load_model
In [19]:
def get_model(mfn, mkwargs):
"""Get the model"""
import datetime
mdir = f"{ddir}/processed/chipnexus/exp/models/multi-task"
name = mfn + "_" + \
",".join([f'{k}={v}' for k,v in mkwargs.items()]) + \
"." + str(datetime.datetime.now()).replace(" ", "::")
!mkdir -p {mdir}
ckp_file = f"{mdir}/{name}.h5"
return eval(mfn)(**mkwargs), name, ckp_file
histogram of the counts¶
- what's the right distribution
- make a q-q plot with poisson, negative bionomial
In [9]:
train[1].shape
Out[9]:
In [13]:
plt.hist(train[1].sum(1).sum(1), bins=100)
plt.xlabel("Total counts");
In [16]:
counts = train[1].sum(1).sum(1)
In [14]:
plt.scatter(train[1].sum(1)[:,0], train[1].sum(1)[:,1])
plt.xlabel("Forward strand counts");
plt.ylabel("Reverse strand counts");
In [21]:
from statsmodels.graphics.gofplots import qqplot
from scipy.stats import poisson, nbinom, norm
In [30]:
import scipy.stats as stats
In [45]:
counts.mean()
Out[45]:
In [68]:
stats.probplot(counts, dist="norm", plot=plt);
In [67]:
# Log-normal seems to model the data pretty well
stats.probplot(np.log(counts+1), dist="norm", plot=plt);
In [63]:
stats.probplot(counts, sparams=(counts.mean()), dist="poisson", plot=plt);
In [54]:
stats.probplot(counts, sparams=(0.2, 0.001), dist="nbinom", plot=plt);
train a model predicting also the counts¶
- focus only on sox2
- make the correlation plot between these two
In [70]:
from keras.models import Sequential
In [128]:
import concise.metrics as ce
import concise.eval_metrics as cem
In [120]:
from concise.metrics import var_explained
In [302]:
filters = 32
def single_output(filters=32,
conv1_kernel_size=21,
dropout=0.1,
seq_len=201,
pool_type='avg',
lr=0.004):
def get_pool(pool_type):
if pool_type=='max':
return kl.MaxPool1D()
elif pool_type =='avg':
return kl.AveragePooling1D()
model = Sequential([
kl.Conv1D(filters, conv1_kernel_size, activation='relu',
input_shape=(seq_len,4), padding='same'),
kl.Conv1D(filters, 1, activation='relu', padding='same'),
get_pool(pool_type),
kl.Conv1D(2*filters, 7, activation='relu', padding='same'),
get_pool(pool_type),
kl.Conv1D(2*filters, 7, activation='relu', padding='same'),
#get_pool(pool_type),
kl.GlobalAveragePooling1D(),
#kl.Dense(8*filters, activation='relu'),
#kl.Dropout(dropout),
kl.Dense(1)
])
model.compile(Adam(lr=lr),
loss="mse",
metrics=[ce.var_explained]
)
return model
In [303]:
def get_model(mfn, mkwargs):
"""Get the model"""
import datetime
mdir = f"{ddir}/processed/chipnexus/exp/models/count-output"
name = mfn + "_" + \
",".join([f'{k}={v}' for k,v in mkwargs.items()]) + \
"." + str(datetime.datetime.now()).replace(" ", "::")
!mkdir -p {mdir}
ckp_file = f"{mdir}/{name}.h5"
return eval(mfn)(**mkwargs), name, ckp_file
In [312]:
# hyper-parameters
mfn = "single_output"
mkwargs = dict(filters=64,
conv1_kernel_size=21,
dropout=0.5,
pool_type='max',
lr=0.004)
In [313]:
valid_counts = (valid[0], np.log(valid[1].sum(1).sum(1)+1))
In [314]:
test_counts = (test[0], np.log(test[1].sum(1).sum(1)+1))
In [315]:
# best valid so far: 108238.6558
model, name, ckp_file = get_model(mfn, mkwargs)
history = model.fit(train[0],
np.log(train[1].sum(1).sum(1)+1),
batch_size=256,
epochs=100,
validation_data=(valid[0], np.log(valid[1].sum(1).sum(1)+1)),
callbacks=[EarlyStopping(patience=5),
History(),
ModelCheckpoint(ckp_file, save_best_only=True)]
)
# get the best model
model = load_model(ckp_file)
In [316]:
dfh = pd.DataFrame(history.history)
plt.figure(figsize=(8,2))
plt.subplot(121)
plt.plot(dfh.loss, label="loss")
plt.plot(dfh.val_loss, label="val_loss")
plt.legend()
plt.xlabel("epoch")
plt.subplot(122)
plt.plot(dfh.var_explained, label="var_explained")
plt.plot(dfh.val_var_explained, label="val_var_explained")
plt.legend()
plt.xlabel("epoch");
In [317]:
y_pred = model.predict(valid_counts[0])
In [318]:
regression_eval(valid_counts[1], y_pred[:,0])
In [311]:
cem.var_explained(valid_counts[1], y_pred)
Out[311]:
Train a joint model¶
In [627]:
def seq_dense_count(filters=21,
conv1_kernel_size=21,
tconv_kernel_size=25,
n_dil_layers=6,
use_profile=True,
seq_len=201,
profile_pool=None,
count_weight=100,
lr=0.004):
"""
Dense
"""
# TODO - build the reverse complement symmetry into the model
inp = kl.Input(shape=(seq_len, 4), name='seq')
first_conv = kl.Conv1D(filters, kernel_size=conv1_kernel_size, padding='same',
#kernel_initializer = ci.PSSMKernelInitializer(pwm_list, stddev=0, add_noise_before_Pwm2Pssm=False),
#bias_initializer = 'zeros',
activation='relu')(inp)
prev_layers = [first_conv]
for i in range(1, n_dil_layers+1):
if i == 1:
prev_sum = first_conv
else:
prev_sum = kl.add(prev_layers)
conv_output = kl.Conv1D(filters, kernel_size=3, padding='same', activation='relu', dilation_rate=2**i)(prev_sum)
prev_layers.append(conv_output)
# De-conv
combined_conv = kl.add(prev_layers)
if profile_pool is not None:
combined_conv = kl.AveragePooling1D(pool_size=profile_pool)(combined_conv)
# Bottleneck layer
#combined_conv = kl.Conv1D(3, kernel_size=1, padding='same', activation='relu')(combined_conv)
x = kl.Reshape((-1, 1, filters))(combined_conv)
x = kl.Conv2DTranspose(2, kernel_size=(tconv_kernel_size, 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), name='profile')(x)
pooled = kl.GlobalAvgPool1D()(combined_conv)
#hidden = kl.Dense(1, pooled)
count_out = kl.Dense(2, name='count')(pooled)
if use_profile:
model = Model(inp, [out, count_out])
model.compile(Adam(lr=lr),
loss=[twochannel_multinomial_nll, 'mse'],
loss_weights=[1, count_weight])
else:
model = Model(inp, count_out)
model.compile(Adam(lr=lr),
loss='mse',)
return model
In [628]:
# hyper-parameters
mfn = "seq_dense_count"
mkwargs = dict(filters=21,
conv1_kernel_size=21,
tconv_kernel_size=25,
n_dil_layers=6,
seq_len=201,
profile_pool=None,
use_profile=False,
lr=0.004)
In [629]:
valid_counts = (valid[0], np.log(valid[1].sum(1).sum(1)+1))
In [630]:
test_counts = (test[0], np.log(test[1].sum(1).sum(1)+1))
In [631]:
model, name, ckp_file = get_model(mfn, mkwargs)
history = model.fit(train[0],
[train[1], np.log(train[1].sum(1)+1)],
batch_size=256,
epochs=100,
validation_data=(valid[0], [np.log(valid[1].sum(1).sum(1)+1)]),
callbacks=[EarlyStopping(patience=5),
History(),
ModelCheckpoint(ckp_file, save_best_only=True)]
)
# get the best model
model = load_model(ckp_file)
In [242]:
y_pred = model.predict(valid[0])
yc_pred = y_pred[:,0]
yc_true = np.log(valid[1].sum(1).sum(1)+1)
In [243]:
cem.var_explained(yc_true, yc_pred)
Out[243]:
In [245]:
regression_eval(yc_true, yc_pred)
In [229]:
regression_eval(yc_true, yc_pred)
Make the plot¶
- train the profile model on a coarser and coarser profile
- report the prediction accuracy
In [320]:
x = train[1]
In [321]:
x.shape
Out[321]:
In [338]:
from basepair.preproc import bin_counts
In [324]:
coarsen = bin_counts(x,3)
In [325]:
coarsen.shape
Out[325]:
In [339]:
plt.plot(x[-1,:,0])
Out[339]:
In [340]:
plt.plot(coarsen[-1,:,0])
Out[340]:
Train the model on the coarsen signal¶
In [361]:
200/4
Out[361]:
In [369]:
200/8
Out[369]:
In [734]:
# hyper-parameters
mfn = "seq_dense_count"
bin_size=None
mkwargs = dict(filters=21,
conv1_kernel_size=21,
tconv_kernel_size=25,
n_dil_layers=6,
seq_len=200,
profile_pool=bin_size,
use_profile=True,
count_weight=100,
lr=0.004)
In [735]:
train, valid, test = seq_inp_exo_out(truncate_len=200, bin_size=bin_size)
In [736]:
valid_counts = (valid[0], np.log(valid[1].sum(1).sum(1)+1))
In [737]:
test_counts = (test[0], np.log(test[1].sum(1).sum(1)+1))
In [738]:
model, name, ckp_file = get_model(mfn, mkwargs)
history = model.fit(train[0],
[train[1], np.log(train[1].sum(1)+1)],
batch_size=256,
epochs=100,
validation_data=(valid[0], [valid[1], np.log(valid[1].sum(1)+1)]),
callbacks=[EarlyStopping(patience=5),
History(),
ModelCheckpoint(ckp_file, save_best_only=True)]
)
# get the best model
model = load_model(ckp_file)
In [739]:
y_pred = model.predict(valid[0])
yc_pred = y_pred[1][:,0]
yc_true = np.log(valid[1].sum(1)[:,0]+1)
In [740]:
cem.var_explained(yc_true, yc_pred)
Out[740]:
In [741]:
regression_eval(yc_true, yc_pred)
In [760]:
y_pred = model.predict(test[0])
yc_pred = y_pred[1][:,0]
yc_true = np.log(test[1].sum(1).sum(1)+1)
In [ ]:
cem.var_explained(yc_true, yc_pred)
In [745]:
ckp_file
Out[745]:
In [744]:
regression_eval(yc_true, yc_pred)
In [654]:
test[1].shape
Out[654]:
In [657]:
y_pred[0][:,:,0].shape
Out[657]:
In [662]:
y_pred[1].shape
Out[662]:
In [665]:
y_pred[0][:,:,0].shape
Out[665]:
In [757]:
binsize=1
regression_eval(np.log(1+np.ravel(bin_counts(test[1], binsize=binsize)[:,:,0])),
np.log(1+np.ravel((bin_counts(softmax(y_pred[0]), binsize=binsize)[:,:,0])*
y_pred[1][:,:1])))
In [755]:
from scipy.stats import pearsonr, spearmanr
In [802]:
binsizes = [1, 2, 4, 10, 25, 50, 100, 200]
perf = []
for binsize in binsizes:
yc_true = np.log(1+np.ravel(bin_counts(test[1], binsize=binsize)[:,:,0]))
yc_pred = np.log(1+np.ravel((bin_counts(softmax(y_pred[0]), binsize=binsize)[:,:,0])*
y_pred[1][:,:1]))
perf.append([pearsonr(yc_true, yc_pred)[0], spearmanr(yc_true, yc_pred)[0]])
In [803]:
df = pd.DataFrame(perf, columns=['pearson','spearman']).assign(binsize=binsizes)
In [811]:
plt.figure(figsize=(3,2))
plt.semilogx(df.binsize, df.spearman,'-o' )
plt.ylabel("R_spearman")
plt.xlabel("Bin size")
plt.xticks(binsizes, binsizes);
In [ ]:
# Show just the plot of different binning
In [710]:
binsize=200//25
regression_eval(np.log(1+np.ravel(bin_counts(test[1], binsize=binsize)[:,:,0])),
np.log(1+np.ravel((bin_counts(softmax(y_pred[0]), binsize=binsize)[:,:,0])*
y_pred[1][:,:1])))
The correlation performance is roughly the same if we train the model on pooled signal (of 25)
In [698]:
R = [0.45, 0.43, 0.39, 0.42, 0.38, 0.42, 0.42, 0.42, 0.38, 0.27]
pool_size = [0, 2, 4, 8, 10, 20, 25, 50, 100, 200]
In [579]:
plt.semilogx(np.array(pool_size)+1, R, "-o")
plt.xlabel("Pool Size")
plt.ylabel("R_spearman")
Out[579]:
In [580]:
plt.plot(np.array(pool_size), R, "-o")
plt.xlabel("Pool Size")
plt.ylabel("R_spearman")
Out[580]:
In [410]:
200/16
Out[410]: