Goal¶
- subcluster the importance scores
- write the instance vectors of the importance scores to file
- use them when running modisco (e.g. subset)
In [1]:
import seaborn as sns
from basepair.config import get_data_dir
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from basepair.modisco import ModiscoResult, append_pattern_loc
from basepair.cli.schemas import DataSpec
from basepair.utils import read_pkl
In [2]:
ddir = Path(get_data_dir())
model_dir = ddir / "processed/chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/"
modisco_dir = model_dir / "modisco/valid"
In [3]:
ds = DataSpec.load(model_dir / "dataspec.yaml")
In [4]:
from kipoi.readers import HDF5Reader
In [5]:
d = HDF5Reader(model_dir / "grad.all.h5")
d.open()
In [6]:
tasks = list(ds.task_specs)
In [77]:
df = pd.DataFrame({task: d.f[f"/targets/profile/{task}"][:].sum(axis=-1).sum(axis=-1) for task in tasks})
In [78]:
df.head()
Out[78]:
In [79]:
dfl = np.log(1+df)
In [7]:
from sklearn.preprocessing import StandardScaler
In [11]:
x = StandardScaler().fit_transform(dfl)
In [12]:
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
principalComponents = pca.fit_transform(x)
principalDf = pd.DataFrame(data = principalComponents, columns = ['pc1', 'pc2'])
In [13]:
dfl.head()
Out[13]:
In [14]:
pca.components_
Out[14]:
In [15]:
principalDf.plot("pc1", "pc2", kind='scatter', alpha=0.1);
In [ ]:
sns.jointplot(x, y, kind="hex", color="#4CB391")
In [27]:
principalDf.plot("pc1", "pc2", gridsize=100, kind='hexbin');
In [29]:
sns.jointplot(principalDf.pc1, principalDf.pc2, kind="hex", color="#4CB391")
Out[29]:
In [16]:
import seaborn as sns
In [36]:
dfls = dfl.sample(frac=0.1)
In [37]:
g = sns.PairGrid(dfls)
g.map_diag(sns.kdeplot)
g.map_offdiag(sns.kdeplot, n_levels=6);
In [18]:
sns.pairplot(dfl,
diag_kind = 'kde',
markers='o',
plot_kws = {'alpha': 0.01, 'edgecolor': None})
Out[18]:
In [ ]:
def plot_emb(emb, legend=False):
"""Plot embedding. Highlights the held-out points on the plot
"""
plt.scatter(emb[:,0], emb[:,1], alpha=0.3, s=60, label="Training")
In [ ]:
plot_emb(tsnec.embedding_)
In [38]:
import umap
In [39]:
um = umap.UMAP(n_neighbors=5,
min_dist=0.3)
In [ ]:
plot_emb(mds.embedding_)
In [ ]:
from sklearn. import Pip
In [8]:
d.ls()
Out[8]:
Cluster activations¶
In [8]:
ls {model_dir}
In [9]:
from keras.models import load_model, Model
In [10]:
from basepair.config import create_tf_session
In [11]:
create_tf_session(0)
Out[11]:
In [12]:
m = load_model(model_dir / 'model.h5')
In [13]:
m.summary()
In [14]:
seq = d.f['/inputs'][:]
In [16]:
import keras.layers as kl
In [18]:
m_avg = Model(m.input, m.get_layer("global_average_pooling1d_1").output)
hid_act = m.get_layer("add_9").output
In [19]:
gmax_pool = kl.GlobalMaxPooling1D()(hid_act)
In [20]:
max_avg_vec = kl.concatenate([gmax_pool, m.get_layer("global_average_pooling1d_1").output])
In [21]:
m_avg_max = Model(m.input, max_avg_vec)
In [22]:
act = m_avg_max.predict(seq, verbose=1)
In [80]:
# append also the total counts
act_w_counts = np.concatenate([act, dfl.values], axis=1)
In [81]:
act_norm = StandardScaler().fit_transform(act_w_counts)
In [82]:
from sklearn.decomposition import PCA
pca = PCA(n_components=4)
principalComponents = pca.fit_transform(act_norm)
principalDf = pd.DataFrame(data = principalComponents, columns = ['pc1', 'pc2', 'pc3', 'pc4'])
In [26]:
idx = pd.Series(np.arange(len(act))).sample(frac=0.1, replace=True)
In [27]:
act_norm_small = act_norm[idx]
In [28]:
len(act_norm_small)
Out[28]:
In [ ]:
act_norm
In [48]:
def plot_emb(emb, legend=False, alpha=0.3, **kwargs):
"""Plot embedding. Highlights the held-out points on the plot
"""
plt.scatter(emb[:,0], emb[:,1], alpha=alpha, s=60, label="Training", **kwargs)
In [30]:
import umap
In [31]:
um = umap.UMAP(n_neighbors=5,
min_dist=0.3)
In [32]:
um.fit(act_norm_small)
Out[32]:
In [38]:
sns.jointplot(um.embedding_[:,0], um.embedding_[:,1], kind="hex", color="#4CB391")
Out[38]:
In [51]:
from sklearn.cluster import AgglomerativeClustering
In [86]:
cl = AgglomerativeClustering(10)
In [87]:
cl.fit(act_norm)
Out[87]:
In [89]:
pd.Series(cl.labels_).value_counts().plot(kind='bar')
Out[89]:
In [ ]:
# Save the clusters to different masked files
In [91]:
outdir = model_dir / "clustering/agglomerative.model,counts/10/"
In [93]:
outdir.mkdir(parents=True, exist_ok=True)
In [94]:
from basepair.utils import write_pkl
In [98]:
for i in range(10):
mask = cl.labels_ == i
np.save(outdir / f"{i}.npy", mask, )
write_pkl(cl, outdir / "clustering.pkl")
In [101]:
ls {outdir}
In [100]:
mask.sum()
Out[100]:
In [55]:
from sklearn import manifold
In [56]:
x = manifold.SpectralEmbedding(n_components=2).fit_transform(act_norm_small)
In [57]:
def plot_clustering(X_red, X, labels, title=None):
x_min, x_max = np.min(X_red, axis=0), np.max(X_red, axis=0)
X_red = (X_red - x_min) / (x_max - x_min)
plt.figure(figsize=(6, 4))
#for i in range(X_red.shape[0]):
plt.scatter(X_red[:, 0], X_red[:, 1], color=plt.cm.nipy_spectral(labels / 10.),
fontdict={'weight': 'bold', 'size': 9})
In [60]:
cl.labels_
Out[60]:
In [62]:
plt.cm.nipy_spectral(cl.labels_ / 10)
Out[62]:
In [ ]:
plt.scatter(x[:,0], x[:,1], color=plt.cm.nipy_spectral(cl.labels_ / 10), alpha=0.01)
In [67]:
cm = plt.get_cmap("tab10")
In [68]:
cm(np.array([10]))
Out[68]:
In [72]:
plt.scatter(um.embedding_[:,0], um.embedding_[:,1], color=cm(cl.labels_), alpha=0.05)
Out[72]:
In [74]:
pd.Series(cl.labels_).value_counts().plot(kind='bar')
Out[74]:
In [71]:
plt.scatter(x[:,0], x[:,1], color=cm(cl.labels_), alpha=0.05)
Out[71]:
In [50]:
plot_emb(um.embedding_, alpha=0.01)
In [42]:
principalDf.plot("pc1", "pc2", kind='scatter', alpha=0.1);
In [84]:
sns.jointplot(principalDf.pc1, principalDf.pc2, kind="hex", color="#4CB391")
Out[84]:
