Questions¶
- can we discover events when Sox2 and Oct4 are both bound?
- could we discover/distinguish those events when using only Sox2 or only Oct4?
- does the multi-task model do better than the single task one?
Goal¶
- make a multi-task model predicting Sox2 and Oct4 ChipNexus signal git issue
TODO¶
- [x] setup a new dataset function based on Sox2 peaks
- [x] setup a new model having two output branches
- [x] train the model
- add the training curve plot bellow training
- [x] modify the plotting code
- show also the output tracks for Oct4
- show also the sequence importance for Oct4
- [ ] try out a few models to figure out the optimal hyper-parameters
- [x] figure out how to best visualize the sequence importances
- max w.r.t. each individual track
- the problem here is that Oct4 might have multiple peaks present in the sequence
- average?
- top2 peaks?
- max w.r.t. each individual track
- [x] run modisco - how do the motifs look like?
- [x] plot the average sequence distribution for each motif
setup a new dataset function based on Sox2 peaks¶
In [1]:
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.datasets import *
from basepair.datasets import sox2_oct4_peaks_sox2
from basepair import samplers
ddir = get_data_dir()
In [2]:
create_tf_session(1)
Out[2]:
In [3]:
train, valid, test = sox2_oct4_peaks_sox2()
setup a new model having two output branches¶
In [4]:
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 [5]:
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 [6]:
i=1
In [7]:
def get_model(mfn, mkwargs, i):
"""Get the model"""
import datetime
mdir = f"{ddir}/processed/chipnexus/exp/models/p-multi-task"
name = mfn + "_" + \
",".join([f'{k}={v}' for k,v in mkwargs.items()]) + \
"." + str(i)
i+=1
!mkdir -p {mdir}
ckp_file = f"{mdir}/{name}.h5"
return eval(mfn)(**mkwargs), name, ckp_file
train the model¶
- add the training curve plot bellow training
In [22]:
# hyper-parameters
mfn = "seq_mutlitask"
use_profile = True
use_counts = True
mkwargs = dict(filters=32,
conv1_kernel_size=21,
tconv_kernel_size=25,
n_dil_layers=6,
use_profile=use_profile,
use_counts=use_counts,
c_task_weight=10,
lr=0.004)
In [23]:
data = sox2_oct4_peaks_sox2()
In [24]:
from basepair.preproc import transform_data
In [25]:
train, valid, test = transform_data(data, use_profile, use_counts)
In [59]:
i
Out[59]:
In [65]:
# best valid so far: 108238.6558
i += 1
model, name, ckp_file = get_model(mfn, mkwargs, i)
history = model.fit(train[0],
train[1],
batch_size=256,
epochs=100,
validation_data=valid[:2],
callbacks=[EarlyStopping(patience=5),
History(),
ModelCheckpoint(ckp_file, save_best_only=True)]
)
# get the best model
model = load_model(ckp_file, custom_objects={"twochannel_multinomial_nll": twochannel_multinomial_nll,
"SpatialLifetimeSparsity": SpatialLifetimeSparsity})
In [90]:
yv_pred = model.predict(test[0])
In [97]:
d
In [98]:
ckp_file
Out[98]:
In [95]:
regression_eval(test[1][1][:,0], yv_pred[1][:,0])
In [103]:
# x10
# val_p_sox2_loss: 258.7645 - val_p_oct4_loss: 363.4225 - val_c_sox2_loss: 0.8482 - val_c_oct4_loss: 1.2600
modify the plotting code¶
- show also the output tracks for Oct4
- show also the sequence importance for Oct4
In [8]:
model = load_model("/users/avsec/workspace/basepair/basepair/../data/processed/chipnexus/exp/models/p-multi-task/seq_mutlitask_filters=32,conv1_kernel_size=21,tconv_kernel_size=25,n_dil_layers=6,use_profile=True,use_counts=True,c_task_weight=10,lr=0.004.5.h5")
In [9]:
yt_pred = model.predict(test[0])
yv_pred = model.predict(valid[0])
In [10]:
from basepair.plots import regression_eval
In [11]:
print(ckp_file)
Valid¶
TODO - the performance of Sox2 could be better¶
In [ ]:
# Sox2
regression_eval(valid[1][2].sum(1), yv_pred[2].sum(1))
In [ ]:
# Oct4
regression_eval(valid[1][3].sum(1), yv_pred[3].sum(1))
Test¶
In [ ]:
# Sox2
regression_eval(test[1][2].sum(1), yt_pred[2].sum(1))
In [ ]:
# Oct4
regression_eval(test[1][3].sum(1), yt_pred[3].sum(1))
In [12]:
from basepair.plots import *
In [123]:
p = Seq2Sox2Oct4(test[0], test[1], model)
In [124]:
p.plot(sort='max_sox2',figsize=(20,2))
In [125]:
p.plot(sort='max_oct4',figsize=(20,2))
try out a few models to figure out the optimal hyper-parameters¶
figure out how to best visualize the sequence importances¶
- max w.r.t. each individual track
- the problem here is that Oct4 might have multiple peaks present in the sequence
- average?
- top2 peaks?
In [19]:
from basepair import samplers
import pandas as pd
from basepair.math import softmax
import numpy as np
import keras.backend as K
from keras.models import Model
from concise.utils.plot import seqlogo_fig, seqlogo
import matplotlib.pyplot as plt
In [127]:
seqlogo
Out[127]:
In [158]:
p = Seq2Nexus(test[0], test[1],test[2], model)
In [159]:
p.plot(sort='max_sox2', seq_grad='max', figsize=(20,12))
run modisco - how do the motifs look like?¶
In [27]:
# Functions to get the gradients
out = model.outputs[0]
inp = model.inputs[0]
pos_strand_ginp_max = K.function([inp], K.gradients(K.max(out[:,:, 0], axis=1), inp))
neg_strand_ginp_max = K.function([inp], K.gradients(K.max(out[:,:, 1], axis=1), inp))
Get predictions and gradients¶
In [28]:
from basepair.data import numpy_minibatch
# Pre-compute the predictions and bottlenecks
x = valid[0]
y_true = valid[1][0]
y_pred = model.predict(x)[0]
grads_pos = np.concatenate([pos_strand_ginp_max([batch])[0]
for batch in numpy_minibatch(x, 512)])
grads_neg = np.concatenate([neg_strand_ginp_max([batch])[0]
for batch in numpy_minibatch(x, 512)])
igrads_pos = grads_pos * x
igrads_neg = grads_neg * x
Correlate pos and neg-strand gradients¶
In [29]:
from scipy.spatial.distance import cosine, correlation
grads_pos_ext = grads_pos.reshape((grads_pos.shape[0], -1))
grads_neg_ext = grads_neg.reshape((grads_neg.shape[0], -1))
distances = np.array([correlation(grads_neg_ext[i], grads_pos_ext[i]) for i in range(len(grads_neg_ext))])
In [30]:
import numpy as np
import seaborn as sns
In [42]:
fig = plt.figure(figsize=(4,2))
plt.hist(distances, bins=50);
plt.xlabel("Importance score cosine distance\n(pos., neg. strand)")
plt.ylabel("Frequency");
plt.savefig('fig/icml18/sox2-similarity.png', dpi=600)
plt.savefig('fig/icml18/sox2-similarity.pdf', dpi=600)
plt.close(fig) # close the figure
fig
Out[42]:
In [31]:
plt.figure(figsize=(6,6))
plt.subplot(211)
plt.hist(distances, bins=50);
plt.subplot(212)
values, base = np.histogram(distances, bins=40)
plt.plot(base[:-1], np.cumsum(values));
plt.grid()
plt.xlabel("Correlation Distance");
plt.ylabel("Fraction of data points");
In [29]:
include = distances<0.2
In [230]:
# We get roughly 1/3 of the points with high-correlation
In [231]:
top10_idx = pd.Series(np.where(distances<0.5)[0]).sample(10)
In [232]:
top10_idx = pd.Series(distances).sort_values(ascending=True).index[:10]
In [233]:
distances
Out[233]:
In [40]:
(distances<0.2).dump(f"{ddir}/processed/chipnexus/motifs/oct4/modisco/valid_include.npy")
In [ ]:
f"{ddir}/processed/chipnexus/motifs/sox2-oct4/modisco/counts.hdf5"
Plot some¶
In [174]:
# Top maxcount indicies
top10_idx = pd.Series(y_true.max(1).sum(1)).sort_values(ascending=False).index[10:30]
In [155]:
# Top count indicies
top10_idx = pd.Series(y_true.sum(1).sum(1)).sort_values(ascending=False).index[:10]
In [156]:
# Random indicies
top10_idx = pd.Series(np.arange(len(y_true))).sample(10)
In [157]:
for i in top10_idx:
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, sharex=True, figsize=(20,6))
ax1.plot(np.arange(1,202), y_true[i,:,0], label="pos")#
ax1.plot(np.arange(1,202), y_true[i,:,1], label="neg")
ax1.set_ylabel("Observed\ncounts")
ax1.legend()
ax2.plot(np.arange(1,202), y_pred[i,:,0], label="pos")#
ax2.plot(np.arange(1,202), y_pred[i,:,1], label="neg")
ax2.set_ylabel("Predicted\n")
ax2.legend()
ax3.set_ylabel("Pos. strand")
seqlogo(igrads_pos[i], ax=ax3);
ax4.set_ylabel("Neg. strand")
seqlogo(igrads_neg[i], ax=ax4);
x_range = [1, 201]
ax4.set_xticks(list(range(0, 201, 5)));
Run modisco¶
In [266]:
top_distances = distances<0.2
hyp_scores = grads_pos + grads_neg
hyp_scores = hyp_scores[top_distances]
hyp_scores = hyp_scores - hyp_scores.mean(-1, keepdims=True)
onehot_data = valid[0][top_distances]
scores = hyp_scores * onehot_data
In [267]:
len(hyp_scores)
Out[267]:
In [268]:
# Visualize
import modisco.visualization
from modisco.visualization import viz_sequence
viz_sequence.plot_weights(scores[0])
viz_sequence.plot_weights(hyp_scores[0])
viz_sequence.plot_weights(onehot_data[0])
In [269]:
from imp import reload
In [270]:
%env MKL_THREADING_LAYER=GNU
In [271]:
import theano
In [272]:
onehot_data.shape
Out[272]:
In [273]:
import h5py
import numpy as np
%matplotlib inline
import modisco
reload(modisco)
import modisco.backend
reload(modisco.backend.theano_backend)
reload(modisco.backend)
import modisco.nearest_neighbors
reload(modisco.nearest_neighbors)
import modisco.affinitymat
reload(modisco.affinitymat.core)
reload(modisco.affinitymat.transformers)
import modisco.tfmodisco_workflow.seqlets_to_patterns
reload(modisco.tfmodisco_workflow.seqlets_to_patterns)
import modisco.tfmodisco_workflow.workflow
reload(modisco.tfmodisco_workflow.workflow)
import modisco.aggregator
reload(modisco.aggregator)
import modisco.cluster
reload(modisco.cluster.core)
reload(modisco.cluster.phenograph.core)
reload(modisco.cluster.phenograph.cluster)
import modisco.core
reload(modisco.core)
import modisco.coordproducers
reload(modisco.coordproducers)
import modisco.metaclusterers
reload(modisco.metaclusterers)
tfmodisco_results = modisco.tfmodisco_workflow.workflow.TfModiscoWorkflow(
sliding_window_size=21,
flank_size=10,
histogram_bins=100,
percentiles_in_bandwidth=10,
overlap_portion=0.5,
min_cluster_size=50,
threshold_for_counting_sign=1.0,
weak_threshold_for_counting_sign=0.7)(
task_names=["task0"],
contrib_scores={'task0': scores},
hypothetical_contribs={'task0': hyp_scores},
one_hot=onehot_data)
In [279]:
mkdir -p {ddir}/processed/chipnexus/motifs/sox2-oct4/modisco
In [275]:
!du -sh {ddir}/processed/chipnexus/motifs/sox2/modisco
In [280]:
modisco_file = f"{ddir}/processed/chipnexus/motifs/sox2-oct4/modisco/counts.hdf5"
modisco_file_pkl = f"{ddir}/processed/chipnexus/motifs/sox2-oct4/modisco/counts.hdf5"
In [281]:
from basepair.utils import write_pkl
In [192]:
write_pkl(tfmodisco_results, modisco_file_pkl)
In [249]:
!rm {modisco_file}
In [282]:
import h5py
import modisco.util
reload(modisco.util)
grp = h5py.File(modisco_file)
tfmodisco_results.save_hdf5(grp)
# TODO - write pickle
In [286]:
hdf5_results.close()
In [283]:
from collections import Counter
from modisco.visualization import viz_sequence
reload(viz_sequence)
from matplotlib import pyplot as plt
import modisco.affinitymat.core
reload(modisco.affinitymat.core)
import modisco.cluster.phenograph.core
reload(modisco.cluster.phenograph.core)
import modisco.cluster.phenograph.cluster
reload(modisco.cluster.phenograph.cluster)
import modisco.cluster.core
reload(modisco.cluster.core)
import modisco.aggregator
reload(modisco.aggregator)
import sklearn.decomposition
import sklearn.manifold
hdf5_results = h5py.File(modisco_file)
#patterns = (tfmodisco_results
# .metacluster_idx_to_submetacluster_results[0]
# .seqlets_to_patterns_result.patterns);
patterns = (list(hdf5_results
["metacluster_idx_to_submetacluster_results"]
["metacluster0"]
["seqlets_to_patterns_result"]
["patterns"]["all_pattern_names"]))
print(len(patterns))
pattern_grp = (hdf5_results
["metacluster_idx_to_submetacluster_results"]
["metacluster0"]
["seqlets_to_patterns_result"]
["patterns"])
for pattern_name in patterns:
pattern = pattern_grp[pattern_name]
print(pattern_name)
print("total seqlets:",len(pattern["seqlets_and_alnmts"]["seqlets"]))
#pattern.plot_counts(counts=aggregated_seqlet.get_per_position_seqlet_center_counts())
background = np.array([0.27, 0.23, 0.23, 0.27])
print("fwd:")
viz_sequence.plot_weights(pattern["task0_contrib_scores"]["fwd"])
viz_sequence.plot_weights(pattern["task0_hypothetical_contribs"]["fwd"])
viz_sequence.plot_weights(viz_sequence.ic_scale(np.array(pattern["sequence"]["fwd"]),
background=background))
print("reverse:")
viz_sequence.plot_weights(pattern["task0_contrib_scores"]["rev"])
viz_sequence.plot_weights(pattern["task0_hypothetical_contribs"]["rev"])
viz_sequence.plot_weights(viz_sequence.ic_scale(np.array(pattern["sequence"]["rev"]),
background=background))
In [285]:
hdf5_results
Out[285]:
In [270]:
sl = tfmodisco_results.multitask_seqlet_creation_results.final_seqlets[0]
In [278]:
cl0 = hdf5_results.metacluster_idx_to_submetacluster_results[0]
In [276]:
sl.coor.example_idx
Out[276]:
In [277]:
len(pd.Series([sl.coor.example_idx for sl in cl0.seqlets]).unique())
Out[277]:
In [255]:
cl0.seqlets[0]
Out[255]:
In [254]:
len(cl0.seqlets)
Out[254]:
In [252]:
len(cl0.seqlets)
Out[252]:
In [ ]:
tfmodisco_results["metacluster_idx_to_submetacluster_results"]
["metacluster0"]
["seqlets_to_patterns_result"]
["patterns"]["all_pattern_names"]))
In [193]:
sl.coor.start
Out[193]:
In [197]:
sl.exidx_start_end_string
Out[197]:
In [222]:
tfmodisco_results.metaclustering_results.metacluster_indices
Out[222]:
In [231]:
sl = tfmodisco_results.multitask_seqlet_creation_results.final_seqlets[0]
In [236]:
sl.coor.
Out[236]:
In [217]:
tfmodisco_results.metaclustering_results.metacluster_indices.max()
Out[217]:
In [219]:
len(tfmodisco_results.metaclustering_results.metacluster_indices)
Out[219]:
In [198]:
sl.
Out[198]: