Goal¶
- develop an activity model predicting:
- H3K27ac
- Pol II
- GRO-seq
Design descisions¶
- dataset = all the chipnexus peaks
- output = total number of counts in the +- 1kb of the region center
Tasks¶
- [x] write the dataloader
- [x] load all the data
- [x] plot the response variable distribution (histogram, QQ-plots)
- [x] get the predictions across all the regions from bpnet
- [x] save to the hdf5 file
- [x] train a simple model on top
- [x] evaluate using simple scatterplots
- [ ] exclude promoter-regions and re-train
Other tasks¶
- [ ] fine-tune the whole model
- [ ] get the importance scores
In [1]:
from basepair.imports import *
import pybedtools
from genomelake.extractors import BigwigExtractor
from basepair.extractors import StrandedBigWigExtractor, bw_extract
from kipoiseq.transforms import ResizeInterval
from basepair.modisco.table import ModiscoData
from pybedtools import BedTool
from basepair.cli.modisco import load_ranges
from basepair.plot.heatmaps import (heatmap_importance_profile, normalize, multiple_heatmap_stranded_profile,
heatmap_stranded_profile)
from basepair.plot.profiles import plot_stranded_profile
import hvplot.pandas
from basepair.plots import regression_eval
In [2]:
create_tf_session(0)
Out[2]:
Load the data¶
In [3]:
df = pd.read_csv(f"{ddir}/processed/chipnexus/external-data.tsv", sep='\t')
df = df.set_index('assay')
In [4]:
df
Out[4]:
Write the bed file
In [5]:
from basepair.cli.schemas import DataSpec
ds = DataSpec.load(f"{ddir}/processed/chipnexus/exp/models/oct-sox-nanog-klf/dataspec.yml")
!mkdir -p {ddir}/processed/activity/data/
def read_factor(factor, filename):
df = pd.read_table(filename, header=None, usecols=[0, 1, 2])
df[3] = factor
df.columns = ['chrom', 'start', 'end', 'name']
return df
dfc = pd.concat([read_factor(k, ds.task_specs[k].peaks) for k in ds.task_specs], axis=0)
dfc.to_csv(f"{ddir}/processed/activity/data/peaks.bed", sep='\t', index=False, header=False)
In [6]:
!head -n 2 {ddir}/processed/activity/data/peaks.bed
Prepare the bigwigs¶
In [7]:
assays = ['H3K27ac', 'PolII', 'Groseq']
In [8]:
def tolist(s):
if isinstance(s, str):
return [s]
else:
return list(s)
In [9]:
bigwigs = {a: tolist(df.loc[a].path) for a in assays}
In [10]:
bigwigs
Out[10]:
Write the dataloader¶
In [13]:
from basepair.config import valid_chr, test_chr
from basepair.datasets import *
In [14]:
dl = ActivityDataset(f"{ddir}/processed/activity/data/peaks.bed", ds.fasta_file, bigwigs,
excl_chromosomes=valid_chr + test_chr)
In [15]:
# load all
train = dl.load_all(num_workers=10)
In [37]:
dfy = pd.DataFrame(train['targets'])
In [92]:
valid = ActivityDataset(f"{ddir}/processed/activity/data/peaks.bed", ds.fasta_file, bigwigs,
incl_chromosomes=valid_chr).load_all(num_workers=10)
In [94]:
dfy_valid = pd.DataFrame(valid['targets'])
In [93]:
test = ActivityDataset(f"{ddir}/processed/activity/data/peaks.bed", ds.fasta_file, bigwigs,
incl_chromosomes=test_chr).load_all(num_workers=10)
In [95]:
dfy_test = pd.DataFrame(test['targets'])
In [38]:
dfy.head()
Out[38]:
In [16]:
len(dl)
Out[16]:
In [17]:
dl[0]
Out[17]:
Output distribution¶
Histograms¶
In [58]:
np.log10(1 + dfy).hvplot.hist(bins=100, subplots=True, height=200, shared_axes=True).cols(1)
Out[58]:
QQ-plots¶
In [59]:
import scipy.stats as stats
In [72]:
fig, axes = plt.subplots(1, len(assays), figsize=(9, 3), sharex=True, sharey=True)
for a, ax in zip(assays, axes):
stats.probplot(np.log10(1 + dfy[a]), dist="norm", plot=ax);
ax.set_title(a)
plt.tight_layout()
Get bottleneck predictions¶
In [34]:
from basepair.BPNet import BPNetPredictor
from keras.models import Model, Sequential
import keras.layers as kl
In [19]:
mdir = f"{ddir}/processed/chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/"
In [20]:
bpnet = BPNetPredictor.from_mdir(mdir)
In [21]:
bpnet.tasks
Out[21]:
In [26]:
bpnet.model.get_layer("add_9").output
Out[26]:
In [27]:
bottleneck_model = Model(bpnet.model.inputs, bpnet.model.get_layer("add_9").output)
In [32]:
bottleneck_predictions = bottleneck_model.predict(train['inputs'], batch_size=32, verbose=1)
In [97]:
bottleneck_predictions_valid = bottleneck_model.predict(valid['inputs'], batch_size=32, verbose=1)
bottleneck_predictions_test = bottleneck_model.predict(test['inputs'], batch_size=32, verbose=1)
In [33]:
bottleneck_predictions.shape
Out[33]:
In [35]:
top_model = Sequential([
kl.GlobalAvgPool1D(input_shape=(1000, 64)),
kl.Dense(3)
])
In [111]:
top_model = Sequential([
kl.MaxPool1D(pool_size=50, input_shape=(1000, 64)),
kl.Flatten(),
kl.Dense(64, activation='relu'),
# kl.Dropout(.5)
kl.Dense(3)
])
In [112]:
from concise.metrics import var_explained
from keras.callbacks import EarlyStopping
from sklearn.preprocessing import StandardScaler
In [113]:
top_model.compile("adam", 'mse', metrics=[var_explained])
In [114]:
preproc = StandardScaler()
In [116]:
y = preproc.fit_transform(np.log10(1 + dfy).values)
y_valid = preproc.transform(np.log10(1 + dfy_valid).values)
y_test = preproc.transform(np.log10(1 + dfy_test).values)
In [117]:
top_model.fit(bottleneck_predictions, y, batch_size=512,
epochs=100,
validation_data=(bottleneck_predictions_valid, y_valid),
callbacks=[EarlyStopping(patience=5)]
)
Out[117]:
Evaluate¶
In [118]:
ypred_valid = top_model.predict(bottleneck_predictions_valid)
In [100]:
ypred_valid.shape
Out[100]:
In [120]:
fig, axes = plt.subplots(len(assays), 1, figsize=(5, 11), sharex=True, sharey=True)
for i, (a, ax) in enumerate(zip(assays, axes)):
regression_eval(ypred_valid[:,i], y_valid[:,i], alpha=0.05, task=a, ax=ax);
plt.tight_layout()
Fine-tune¶
In [122]:
whole_model = Sequential([bottleneck_model, top_model])
In [123]:
whole_model.compile("adam", 'mse', metrics=[var_explained])
In [127]:
whole_model.fit(train['inputs']['seq'], y, batch_size=512,
epochs=100,
validation_data=(valid['inputs']['seq'], y_valid),
callbacks=[EarlyStopping(patience=5)])
Out[127]:
In [128]:
ypred_valid = whole_model.predict(valid['inputs']['seq'])
In [129]:
ypred_valid.shape
Out[129]:
In [130]:
fig, axes = plt.subplots(len(assays), 1, figsize=(5, 11), sharex=True, sharey=True)
for i, (a, ax) in enumerate(zip(assays, axes)):
regression_eval(ypred_valid[:,i], y_valid[:,i], alpha=0.05, task=a, ax=ax);
plt.tight_layout()
Train from scratch¶
In [141]:
import keras.backend as K
from keras.models import load_model
def reset_weights(model):
session = K.get_session()
for layer in model.layers:
if hasattr(layer, 'kernel_initializer'):
layer.kernel.initializer.run(session=session)
In [137]:
whole_model.save(f"{ddir}/processed/activity/models/dense/fine-tuned.h5")
In [145]:
reinitialized_model = load_model(f"{ddir}/processed/activity/models/dense/fine-tuned.h5")
In [163]:
reset_weights(reinitialized_model.layers[0])
reset_weights(reinitialized_model.layers[1])
In [165]:
reinitialized_model.compile("adam", 'mse', metrics=[var_explained])
In [166]:
reinitialized_model.fit(train['inputs']['seq'], y, batch_size=512,
epochs=100,
validation_data=(valid['inputs']['seq'], y_valid),
callbacks=[EarlyStopping(patience=5)])
Out[166]:
In [167]:
ypred_valid = reinitialized_model.predict(valid['inputs']['seq'])
In [168]:
ypred_valid.shape
Out[168]:
In [169]:
fig, axes = plt.subplots(len(assays), 1, figsize=(5, 11), sharex=True, sharey=True)
for i, (a, ax) in enumerate(zip(assays, axes)):
regression_eval(ypred_valid[:,i], y_valid[:,i], alpha=0.05, task=a, ax=ax);
plt.tight_layout()
