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 [23]:
from basepair.config import valid_chr, test_chr
from basepair.datasets import *
In [53]:
ds.fasta_file
Out[53]:
In [30]:
dl = ActivityDataset(f"{ddir}/processed/activity/data/peaks.bed", ds.fasta_file, bigwigs,
excl_chromosomes=valid_chr + test_chr)
# load all
train = dl.load_all(num_workers=10)
dfy = pd.DataFrame(train['targets'])
valid = ActivityDataset(f"{ddir}/processed/activity/data/peaks.bed", ds.fasta_file, bigwigs,
incl_chromosomes=valid_chr).load_all(num_workers=10)
dfy_valid = pd.DataFrame(valid['targets'])
test = ActivityDataset(f"{ddir}/processed/activity/data/peaks.bed", ds.fasta_file, bigwigs,
incl_chromosomes=test_chr).load_all(num_workers=10)
dfy_test = pd.DataFrame(test['targets'])
In [31]:
dfy.head()
Out[31]:
In [20]:
len(dl)
Out[20]:
In [21]:
dl[0]
Out[21]:
Output distribution¶
Histograms¶
QQ-plots¶
In [23]:
import scipy.stats as stats
In [24]:
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 [11]:
from basepair.BPNet import BPNetPredictor
from keras.models import Model, Sequential
import keras.layers as kl
In [12]:
dfe = pd.read_csv("https://docs.google.com/spreadsheets/d/1n3l2HXKSNpmNUOifD41uRzDEAgmOqXMQDxquRaz6WLg/export?gid=0&format=csv")
In [13]:
dfe = dfe[dfe.train == True]
In [14]:
model_exps = {
#"nexus/binary.gw": 'nexus,gw,OSNK,1,0,0,FALSE,valid,0.5,64,25,0.001,9,FALSE',
"nexus/profile.peaks.bias-corrected": 'nexus,peaks,OSNK,0,10,1,FALSE,same,0.5,64,25,0.004,9,FALSE,[1,50],TRUE',
"nexus/profile.peaks.non-bias-corrected": 'nexus,peaks,OSNK,0,10,1,FALSE,same,0.5,64,25,0.004,9,FALSE-2',
#"seq/binary.gw": 'seq,gw,OSN,1,0,0,FALSE,valid,0.5,64,50,0.001,9,FALSE/',
"seq/profile.peaks.bias-corrected": 'seq,peaks,OSN,0,10,1,FALSE,same,0.5,64,50,0.004,9,FALSE,[1,50],TRUE',
"seq/profile.peaks.non-bias-corrected": 'seq,peaks,OSN,0,10,1,FALSE,same,0.5,64,50,0.004,9,FALSE',
'nexus/profile.peaks-union.bias-corrected': 'nexus,nexus-seq-union,OSN,0,10,1,FALSE,same,0.5,64,25,0.004,9,FALSE,[1,50],TRUE',
'seq/profile.peaks-union.bias-corrected': ' seq,nexus-seq-union,OSN,0,10,1,FALSE,same,0.5,64,50,0.004,9,FALSE,[1,50],TRUE',
}
In [15]:
ls ../../chipnexus/train/seqmodel
In [16]:
base_mdir = Path('/oak/stanford/groups/akundaje/avsec/basepair/data/processed/comparison/output')
In [17]:
mdir = base_mdir / model_exps['nexus/profile.peaks.bias-corrected']
In [18]:
from basepair.seqmodel import SeqModel
In [19]:
bpnet = SeqModel.from_mdir(mdir)
In [20]:
import tensorflow as tf
In [21]:
bottleneck_model = bpnet.bottleneck_model()
In [22]:
bottleneck_model.predict(np.ones((1, 1000, 4))).shape
Out[22]:
In [34]:
bottleneck_predictions = bottleneck_model.predict(train['inputs'], batch_size=32, verbose=1)
In [35]:
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 [35]:
top_model = Sequential([
kl.GlobalAvgPool1D(input_shape=(1000, 64)),
kl.Dense(3)
])
In [36]:
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 [37]:
from concise.metrics import var_explained
from keras.callbacks import EarlyStopping
from sklearn.preprocessing import StandardScaler
In [14]:
mdir = '../../../src/chipnexus/train/seqmodel/output/nexus,peaks,OSNK,0,10,1,FALSE,same,0.5,64,25,0.004,9,FALSE,[1,50],TRUE/'
In [15]:
# mdir = '../../../src/chipnexus/train/seqmodel/output/nexus,gw,OSNK,1,0,0,FALSE,same,0.5,64,25,0.001,9,FALSE/'
# mdir = '../../../src/chipnexus/train/seqmodel/output/seq,peaks,OSN,0,10,1,FALSE,same,0.5,64,50,0.004,9,FALSE,[1,50],TRUE/'
# mdir = '../../../src/chipnexus/train/seqmodel/output/nexus,peaks,OSNK,0,10,1,FALSE,same,0.5,64,25,0.004,9,FALSE-2/'
In [38]:
top_model.compile("adam", 'mse', metrics=[var_explained])
In [39]:
preproc = StandardScaler()
In [40]:
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 [41]:
top_model.fit(bottleneck_predictions, y, batch_size=512,
epochs=100,
validation_data=(bottleneck_predictions_valid, y_valid),
callbacks=[EarlyStopping(patience=5, restore_best_weights=True)]
)
Out[41]:
Evaluate¶
In [42]:
ypred_valid = top_model.predict(bottleneck_predictions_valid)
In [43]:
ypred_valid.shape
Out[43]:
In [44]:
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()
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 [45]:
whole_model = Sequential([bottleneck_model, top_model])
In [46]:
whole_model.compile("adam", 'mse', metrics=[var_explained])
In [47]:
whole_model.fit(train['inputs']['seq'], y, batch_size=512,
epochs=100,
validation_data=(valid['inputs']['seq'], y_valid),
callbacks=[EarlyStopping(patience=5, restore_best_weights=True)])
Out[47]:
In [48]:
ypred_valid = whole_model.predict(valid['inputs']['seq'])
In [ ]:
ypred_valid.shape
In [ ]:
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()
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()
