Goal¶
- cluster the motif sites using modisco
Resources¶
Modisco videos¶
- NIPS MLCB 2017 TF MoDISco - https://www.youtube.com/watch?v=fXPGVJg956E
- Modisco MVP - https://www.youtube.com/watch?v=-f8N6XNbtLU&feature=youtu.be
TODO¶
- get the importance profiles
- add-up the grads for two strands
In [2]:
import os
import numpy as np
from keras.models import load_model
from basepair.config import create_tf_session, get_data_dir
from basepair import models, datasets
import basepair as bp
from basepair.losses import twochannel_multinomial_nll
from basepair.layers import SpatialLifetimeSparsity
from basepair.math import softmax
from basepair.plots import *
from keras.models import load_model
from tqdm import tqdm
ddir = get_data_dir()
from basepair.data import numpy_minibatch
Setup the model and the data¶
In [3]:
create_tf_session(1)
Out[3]:
In [4]:
ckp_file = f"{ddir}/processed/chipnexus/exp/models/resnest_allconnect_nconv=7_filters=32_lr=0.004_dilated=True,out=25.h5"
In [5]:
model = load_model(ckp_file)
In [6]:
train, valid, test = datasets.seq_inp_exo_out()
Plot¶
In [33]:
# Functions to get the gradients
out = model.outputs[0]
inp = model.inputs[0]
pos_strand_ginp_avg = K.function([inp], K.gradients(K.mean(out[:,:,0], axis=-1), inp))
neg_strand_ginp_avg = K.function([inp], K.gradients(K.mean(out[:,:,1], axis=-1), inp))
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 [36]:
# Pre-compute the predictions and bottlenecks
x = valid[0]
y_true = valid[1]
y_pred = softmax(model.predict(x))
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 [13]:
from scipy.spatial.distance import cosine, correlation
In [38]:
grads_pos_ext = grads_pos.reshape((grads_pos.shape[0], -1))
grads_neg_ext = grads_neg.reshape((grads_neg.shape[0], -1))
In [39]:
distances = np.array([correlation(grads_neg_ext[i], grads_pos_ext[i]) for i in range(len(grads_neg_ext))])
In [40]:
import numpy as np
import seaborn as sns
In [41]:
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 [ ]:
# We get roughly 1/3 of the points with high-correlation
In [42]:
top10_idx = pd.Series(np.where(distances<0.5)[0]).sample(10)
In [119]:
top10_idx = pd.Series(distances).sort_values(ascending=True).index[:10]
In [68]:
distances
Out[68]:
Plot some¶
In [49]:
# Top maxcount indicies
top10_idx = pd.Series(y_true.max(1).sum(1)).sort_values(ascending=False).index[10:30]
In [46]:
# Top count indicies
top10_idx = pd.Series(y_true.sum(1).sum(1)).sort_values(ascending=False).index[:10]
In [36]:
# Random indicies
top10_idx = pd.Series(np.arange(len(y_true))).sample(10)
In [44]:
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 [53]:
top_distances = distances<0.2
In [54]:
hyp_scores = grads_pos + grads_neg
In [55]:
hyp_scores = hyp_scores[top_distances]
In [56]:
hyp_scores = hyp_scores - hyp_scores.mean(-1, keepdims=True)
In [61]:
onehot_data = valid[0][top_distances]
In [59]:
scores = hyp_scores * one_hot_data
In [58]:
len(hyp_scores)
Out[58]:
In [62]:
# 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 [63]:
from imp import reload
In [64]:
%env MKL_THREADING_LAYER=GNU
In [65]:
import theano
In [66]:
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()(
task_names=["task0"],
contrib_scores={'task0': scores},
hypothetical_contribs={'task0': hyp_scores},
sliding_window_size=21,
flank_size=10,
histogram_bins=100,
percentiles_in_bandwidth=10,
overlap_portion=0.5,
min_cluster_size=200,
threshold_for_counting_sign=1.0,
weak_threshold_for_counting_sign=0.7,
one_hot=onehot_data)
In [74]:
mkdir -p {ddir}/processed/chipnexus/motifs/sox2/modisco
In [76]:
!du -sh {ddir}/processed/chipnexus/motifs/sox2/modisco
In [77]:
modisco_file = f"{ddir}/processed/chipnexus/motifs/sox2/modisco/defaults.hdf5"
In [75]:
import h5py
import modisco.util
reload(modisco.util)
grp = h5py.File(modisco_file)
tfmodisco_results.save_hdf5(grp)
In [78]:
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))
