Goal¶
- [x] analyze the data with BPNet
Notes¶
- to run BPNet, we should probably inject this sequence in the middle of a random sequence and average the predictions across multiple random sequences
Tasks¶
- [x] get the importance scores for all these sequences
- [x] visualize the importance scores as well as the footprints
- [x] correlate importance scores with activity
- [X] build a simple predictive model
- [x] compare to their RF model
- genomic sequences: gkSVM (top 25% vs bottom 75%) and got AUROC and AUPRC in the neighborhood of 0.75.
- See Figure 4
- genomic sequences: gkSVM (top 25% vs bottom 75%) and got AUROC and AUPRC in the neighborhood of 0.75.
Conclusions¶
- our model is not that good as theirs on synthetic constructs (however, their model greatly overfits the data (knows the location)
- model trained on genomics sequences achieves auc=0.83 instead of theirs 0.75
In [73]:
import warnings
warnings.filterwarnings("ignore")
In [1]:
from basepair.exp.chipnexus.data import load_BarakCohnen_mESC_MPRA, load_BarakCohnen_mESC_MPRA_raw
In [2]:
dfs_gen, dfs_syn = load_BarakCohnen_mESC_MPRA()
df_gen, df_syn = load_BarakCohnen_mESC_MPRA_raw()
In [3]:
dfs_gen.head()
Out[3]:
In [4]:
dfs_syn.head()
Out[4]:
In [5]:
df_gen.head()
Out[5]:
In [6]:
df_gen.groupby("CRE_id").size().value_counts()
Out[6]:
There are 8 different BC-seq per sequence
Create full sequences¶
Gen¶
Genomic sequences were represented in the pool by 150 bp oligos with the following sequences:
GACTTACATTAGGGCCCGT[SEQ]AAGCTT[FILL]GAATTCTCTAGAC[BC]TGAGCTCGGACTACGATACTG
Where [SEQ] is either a reference (gWT) or mutated (gMUT) genomic sequence of 81-82 bps
Syn¶
CTTCTACTACTAGGGCCCA[SEQ]AAGCTT[FILL]GAATTCTCTAGAC[BC]TGAGCTCTACATGCTAGTTCATG
where [SEQ] is a 40-80 bp synthetic element comprised of concatenated 20 bp building blocks
of pluripotency sites, as described previously, with the fifth position of the KLF4 site changed to
‘T’ to facilitate cloning (Fiore and Cohen 2016) . [FILL] is a random filler sequence of variable
length to bring the total length of each sequence to 150 bp, and [BC] is a random 9 bp barcode.
In [118]:
from basepair.exp.chipnexus.simulate import generate_seq, random_seq
def gen_padded_sequence(sequence, barcode, seqlen=1000):
"""Generate MPRA construct sequence for Genomic sequences
"""
fill = random_seq(150 - len(sequence) - len(barcode) - 59)
core_seq = f"GACTTACATTAGGGCCCGT{sequence}AAGCTT{fill}GAATTCTCTAGAC{barcode}TGAGCTCGGACTACGATACTG"
assert len(core_seq) == 150
return generate_seq(core_seq, seqlen=seqlen)
def syn_padded_sequence(sequence, barcode, seqlen=1000):
"""Generate MPRA construct sequence for Synthetic sequences
"""
fill = random_seq(150 - len(sequence) - len(barcode) - 61)
core_seq = f'CTTCTACTACTAGGGCCCA{sequence}AAGCTT{fill}GAATTCTCTAGAC{barcode}TGAGCTCTACATGCTAGTTCATG'
assert len(core_seq) == 150
return generate_seq(core_seq, seqlen=seqlen)
def get_core_seq_ij(seqlen, core_seq_len=150):
start = (seqlen - core_seq_len) // 2
return start, start + core_seq_len
In [8]:
a = gen_padded_sequence("ACGTA", "ACGST")
assert len(a) == 1000
In [9]:
b = syn_padded_sequence("ACGTA", "ACGST")
assert len(a) == 1000
Get the importance scores for all these sequences¶
In [10]:
from basepair.BPNet import BPNetPredictor
In [11]:
# TODO - change to {ddir}/...
mdir = '/s/project/avsec/basepair/models/oct-sox-nanog-klf'
In [12]:
ls {mdir}
In [13]:
bpnet = BPNetPredictor.from_mdir(mdir)
In [14]:
bpnet.tasks
Out[14]:
first, don't average across multiple flanking sequences
In [15]:
from concise.preprocessing import encodeDNA
In [16]:
from concise.preprocessing.sequence import one_hot2string, DNA
Predict¶
In [20]:
from basepair.utils import write_pkl, read_pkl
from pathlib import Path
In [21]:
cache_dir = Path("/s/project/avsec/basepair/BarakCohen-mESC-MPRA/cache")
In [22]:
cache=True
In [24]:
if cache and os.path.exists(cache_dir / "bpnet_preds.pkl"):
d = read_pkl(cache_dir / "bpnet_preds.pkl")
bpnet_seq_preds_gen, bpnet_seq_preds_syn = d['gen'], d['syn']
else:
bpnet_seq_syn = encodeDNA([syn_padded_sequence(s, "AAAGACGCG") for s in dfs_syn.Sequence.str.upper()])
bpnet_seq_gen = encodeDNA([gen_padded_sequence(s, "AAAGACGCG") for s in dfs_gen.Sequence.str.upper()])
bpnet_seq_preds_gen = bpnet.seq_predict(bpnet_seq_gen)
bpnet_seq_preds_syn = bpnet.seq_predict(bpnet_seq_syn)
write_pkl(dict(gen=bpnet_seq_preds_gen, syn=bpnet_seq_preds_syn), cache_dir / "bpnet_preds.pkl")
In [ ]:
!du -sh {cache_dir}/
Visualize the importance scores as well as the footprints¶
In [25]:
from basepair.plot.tracks import plot_tracks, filter_tracks
In [26]:
from collections import OrderedDict
def bpnet_format_output(obj, tasks, output_features=['p', 'i'], outlen=150):
out = dict(p=obj['preds']['scaled_profile'],
i={t: obj['seq'] * (obj['grads'][ti]['profile']['pos'] + obj['grads'][ti]['profile']['neg'])/2
for ti,t in enumerate(tasks)})
to_plot = OrderedDict([(f"{task}/{t}", out[t][task])
for t in output_features
for task in tasks
])
seqlen = obj['seq'].shape[0]
start = (seqlen - outlen) // 2
return filter_tracks(to_plot, xlim=[start, start + outlen])
In [27]:
max([np.abs(bpnet_format_output(bpnet_seq_preds_syn[i], bpnet.tasks, output_features=['i'])[f'{task}/i']).max() for task in bpnet.tasks for i in range(len(bpnet_seq_preds_syn))])
Out[27]:
In [28]:
max([np.abs(bpnet_format_output(bpnet_seq_preds_gen[i], bpnet.tasks, output_features=['i'])[f'{task}/i']).max() for task in bpnet.tasks for i in range(len(bpnet_seq_preds_gen))])
Out[28]:
In [29]:
# pre-compute the ranges for the importance scores
ylim_max = max(max([np.abs(bpnet_format_output(bpnet_seq_preds_syn[i], bpnet.tasks, output_features=['i'])[f'{task}/i']).max() for task in bpnet.tasks for i in range(len(bpnet_seq_preds_syn))]),
max([np.abs(bpnet_format_output(bpnet_seq_preds_gen[i], bpnet.tasks, output_features=['i'])[f'{task}/i']).max() for task in bpnet.tasks for i in range(len(bpnet_seq_preds_gen))]))
In [31]:
dfs_gen.head()
Out[31]:
In [244]:
mean_vals = OrderedDict([(k,v.mean()) for k,v in bpnet_format_output(bpnet_seq_preds_gen[i], bpnet.tasks, output_features=['p', 'i']).items()])
Genomic¶
In [238]:
# 10 most expressed sequences
for i in dfs_gen.expression_fc.sort_values().tail(10).index:
plot_tracks(bpnet_format_output(bpnet_seq_preds_gen[i], bpnet.tasks, output_features=['i']), # use output_features=['p', i'] to also show the profile
fig_width=14, fig_height_per_track=1.5,
ylim=(-ylim_max, ylim_max),
title=dfs_gen.iloc[i].CRE_id);
In [240]:
# 5 least expressed sequences
for i in dfs_gen.expression_fc.sort_values().head(5).index:
plot_tracks(bpnet_format_output(bpnet_seq_preds_gen[i], bpnet.tasks, output_features=['i']), # use output_features=['p', i'] to also show the profile
fig_width=14, fig_height_per_track=1.5,
ylim=(-ylim_max, ylim_max),
title=dfs_gen.iloc[i].CRE_id);
Synthetic¶
In [239]:
# for i in pd.Series(np.arange(len(bpnet_seq_preds_syn))).sample(10):
for i in dfs_syn.expression_fc.sort_values().tail(10).index:
plot_tracks(bpnet_format_output(bpnet_seq_preds_syn[i], bpnet.tasks, output_features=['i']), # use output_features=['p', i'] to also show the profile
fig_width=14, fig_height_per_track=1.5,
ylim=(-ylim_max, ylim_max),
title=dfs_syn.iloc[i].CRE_elements);
In [241]:
# 5 least expressed constructrs
for i in dfs_syn.expression_fc.sort_values().head(5).index:
plot_tracks(bpnet_format_output(bpnet_seq_preds_syn[i], bpnet.tasks, output_features=['i']), # use output_features=['p', i'] to also show the profile
fig_width=14, fig_height_per_track=1.5,
ylim=(-ylim_max, ylim_max),
title=dfs_syn.iloc[i].CRE_elements);
Correlate with importance scores¶
In [33]:
values_gen = pd.DataFrame([OrderedDict([(k + "/" + fn.__name__,fn(v)) for k,v in bpnet_format_output(bpnet_seq_preds_gen[i], bpnet.tasks, output_features=['p', 'i']).items() for fn in [np.max, np.mean]])
for i in range(len(bpnet_seq_preds_gen))])
In [34]:
values_syn = pd.DataFrame([OrderedDict([(k + "/" + fn.__name__,fn(v)) for k,v in bpnet_format_output(bpnet_seq_preds_syn[i], bpnet.tasks, output_features=['p', 'i']).items() for fn in [np.max, np.mean]])
for i in range(len(bpnet_seq_preds_syn))])
In [35]:
for f in ['p/amax', 'i/amax', 'p/mean', 'i/mean']:
values_gen[f] = values_gen[[t + "/" + f for t in bpnet.tasks]].sum(axis=1)
In [36]:
for f in ['p/amax', 'i/amax', 'p/mean', 'i/mean']:
values_syn[f] = values_syn[[t + "/" + f for t in bpnet.tasks]].sum(axis=1)
In [37]:
from basepair.modisco.motif_clustering import df_minmax_scale
In [38]:
dfs_gen_val = pd.concat([df_minmax_scale(values_gen), dfs_gen], 1)
In [39]:
dfs_syn_val = pd.concat([df_minmax_scale(values_syn), dfs_syn], 1)
In [40]:
dfs_gen_val.head()
Out[40]:
In [41]:
dfm_gen = dfs_gen_val.melt(['expression_fc_log2', 'CRE_type', 'interval'], list(values_gen))
dfm_gen.head()
Out[41]:
In [42]:
dfm_syn = dfs_syn_val.melt(['expression_fc_log2', 'CRE_elements'], list(values_syn))
dfm_syn.head()
Out[42]:
In [50]:
import matplotlib as mpl
In [55]:
mpl.style.use('default')
In [87]:
dfm_gen.groupby(['CRE_type', 'variable'])[['expression_fc_log2', 'value']].corr('spearman').iloc[0::2,-1].plot.barh(figsize=(5, 8))
Out[87]:
In [85]:
dfm_syn.groupby(['variable'])[['expression_fc_log2', 'value']].corr('spearman').iloc[0::2,-1].plot.barh(figsize=(3, 5))
Out[85]:
In [94]:
fig, ax = plt.subplots(figsize=(4,3))
dfs_gen_val.plot.scatter('i/amax', 'expression_fc_log2', alpha=0.2, ax=ax);
In [75]:
import plotnine
from plotnine import *
In [81]:
plotnine.options.figure_size = (25, 4)
ggplot(aes(x='value', y='expression_fc_log2'), dfm_gen[dfm_gen.CRE_type=='Genomic']) + geom_point(alpha=0.1) + geom_smooth(alpha=0.1, color='grey', linetype='dashed', method='lm') + facet_wrap("~variable", ncol=10) +theme_classic()
Out[81]:
In [80]:
plotnine.options.figure_size = (25, 4)
ggplot(aes(x='value', y='expression_fc_log2'), dfm_gen[dfm_gen.CRE_type=='All_Mutated']) + geom_point(alpha=0.1) + geom_smooth(alpha=0.1, color='grey', linetype='dashed', method='lm') + facet_wrap("~variable", ncol=10) +theme_classic()
Out[80]:
In [79]:
plotnine.options.figure_size = (25, 4)
ggplot(aes(x='value', y='expression_fc_log2'), dfm_syn) + geom_point(alpha=0.1) + geom_smooth(alpha=0.1, color='grey', linetype='dashed', method='lm') + facet_wrap("~variable", ncol=10) + theme_classic()
Out[79]:
Score motif instances¶
Build a simple predictive model based on motif instance locations¶
In [95]:
bpnet.model
Out[95]:
In [111]:
from keras.models import Model, Sequential
In [112]:
import keras.layers as kl
In [104]:
bpnet.model.get_layer("reshape_1").output
Out[104]:
In [105]:
bottleneck_model = Model(bpnet.model.inputs, bpnet.model.get_layer("reshape_1").output)
In [106]:
bpnet_seq_syn = encodeDNA([syn_padded_sequence(s, "AAAGACGCG") for s in dfs_syn.Sequence.str.upper()])
bpnet_seq_gen = encodeDNA([gen_padded_sequence(s, "AAAGACGCG") for s in dfs_gen.Sequence.str.upper()])
In [113]:
import skimage.measure
In [107]:
%time bpnet_bottleneck_syn = bottleneck_model.predict(bpnet_seq_syn)
In [108]:
%time bpnet_bottleneck_gen = bottleneck_model.predict(bpnet_seq_gen)
In [110]:
bpnet_bottleneck_syn.shape
Out[110]:
In [318]:
# get the features
def pool_bottleneck(bottleneck, pool_size=37, outlen=80, agg_fn=np.max):
i,j = get_core_seq_ij(bottleneck.shape[1])
i += 19
j = i + outlen
out = skimage.measure.block_reduce(bottleneck[:,i:j], (1, pool_size, 1, 1), agg_fn)
return out.reshape((out.shape[0], -1))
In [296]:
bpnet_bottleneck_syn_gen.shape
Out[296]:
In [297]:
## train a simple model
In [553]:
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.linear_model import LassoCV
from sklearn.model_selection import cross_val_predict
from basepair.plots import regression_eval
Synthetic¶
In the paper, they achieved R-squared of 0.87
In [609]:
bpnet_bottleneck_syn_feat = pool_bottleneck(bpnet_bottleneck_syn, pool_size=10, outlen=80, agg_fn=np.max)
# add also the individual features
bpnet_bottleneck_syn_feat = np.concatenate([values_syn.values, bpnet_bottleneck_syn_feat], 1)
preds = cross_val_predict(RandomForestRegressor(100), bpnet_bottleneck_syn_feat, dfs_syn.expression_fc_log2.values, cv=5, n_jobs=-1)
regression_eval(preds, dfs_syn.expression_fc_log2.values, alpha=0.4);
Genomic¶
In [612]:
bpnet_bottleneck_gen_feat = pool_bottleneck(bpnet_bottleneck_gen, 10, outlen=80)
bpnet_bottleneck_gen_feat = np.concatenate([values_gen.values, bpnet_bottleneck_gen_feat], axis=1)
preds = cross_val_predict(RandomForestRegressor(100), bpnet_bottleneck_gen_feat, dfs_gen.expression_fc_log2.values, cv=5, n_jobs=-1)
regression_eval(dfs_gen.expression_fc_log2, preds, alpha=0.4);
In [613]:
subset = dfs_gen.CRE_type=='Genomic'
regression_eval(dfs_gen.expression_fc_log2[subset], preds[subset], alpha=0.4);
Cross-data comparison¶
Trained on genetic data, eval on synthetic¶
In [614]:
rf_gen = RandomForestRegressor(100, n_jobs=-1)
rf_gen.fit(bpnet_bottleneck_gen_feat, dfs_gen.expression_fc_log2.values)
Out[614]:
In [615]:
regression_eval(dfs_syn.expression_fc_log2, rf_gen.predict(bpnet_bottleneck_syn_feat), alpha=0.4);
Trained on synthetic data, eval on genetic¶
In [616]:
rf_syn = RandomForestRegressor(100)
rf_syn.fit(bpnet_bottleneck_syn_feat, dfs_syn.expression_fc_log2.values)
regression_eval(dfs_gen.expression_fc_log2, rf_syn.predict(bpnet_bottleneck_gen_feat), alpha=0.4);
Compare to their RF model¶
- genomic sequences: gkSVM (top 25% vs bottom 75%) and got AUROC and AUPRC in the neighborhood of 0.75.
- See Figure 4
In [617]:
from basepair.stats import perc
from concise.eval_metrics import auc, auprc
In [618]:
genomic = dfs_gen.CRE_type=='Genomic'
In [619]:
neg = perc(dfs_gen.expression_fc_log2.values[genomic])< 25
pos = perc(dfs_gen.expression_fc_log2.values[genomic])> 75
In [620]:
subset = pos | neg
In [621]:
y = pos[subset]
In [622]:
y.mean()
Out[622]:
Bottleneck features¶
In [623]:
bpnet_bottleneck_gen_feat = pool_bottleneck(bpnet_bottleneck_gen, 5, outlen=80)
#bpnet_bottleneck_gen_feat = np.concatenate([values_gen.values, bpnet_bottleneck_gen_feat], axis=1)
preds = cross_val_predict(RandomForestClassifier(1000), bpnet_bottleneck_gen_feat[genomic][subset], y, cv=3, n_jobs=-1, method='predict_proba')
print("auc, auprc")
auc(y, preds[:,1]), auprc(y, preds[:,1])
Out[623]:
Both¶
In [624]:
bpnet_bottleneck_gen_feat = pool_bottleneck(bpnet_bottleneck_gen, 5, outlen=80)
bpnet_bottleneck_gen_feat = np.concatenate([values_gen.values, bpnet_bottleneck_gen_feat], axis=1)
preds = cross_val_predict(RandomForestClassifier(1000), bpnet_bottleneck_gen_feat[genomic][subset], y, cv=3, n_jobs=-1, method='predict_proba')
print("auc, auprc")
auc(y, preds[:,1]), auprc(y, preds[:,1])
Out[624]:
Only importance scores¶
In [625]:
bpnet_bottleneck_gen_feat = values_gen.values
preds = cross_val_predict(RandomForestClassifier(1000), bpnet_bottleneck_gen_feat[genomic][subset], y, cv=3, n_jobs=-1, method='predict_proba')
print("auc, auprc")
auc(y, preds[:,1]), auprc(y, preds[:,1])
Out[625]: