Goal¶
- analyze instance locations. Focus on Oct4 peaks
TODO¶
- sample a few points from the scatterplot: importance vs match to see what you are getting
- try the same analysis for a few factors
- try to find the right cutoff based on the number of motifs you expect
- eyeball a few regions (higlight the found seqlets in there)
- before doing the search, merge all teh metaclusters (ignore the importance score magnitude)
Conclusions¶
- ignore hypothetical scores when scoring seqlet instances
In [212]:
from basepair.imports import *
from basepair.samplers import top_sum_count, top_max_count
from basepair.plot.tracks import plot_tracks, filter_tracks
from kipoi.data_utils import get_dataset_item
In [213]:
model_dir = Path(f"{ddir}/processed//chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/")
In [214]:
modisco_dir = model_dir / "modisco/by_peak_tasks/weighted/Oct4"
In [215]:
mr = ModiscoResult(modisco_dir / "modisco.h5")
mr.open()
In [216]:
tasks = ['Oct4', 'Sox2', 'Nanog', 'Klf4']
In [217]:
mr.plot_pssm("metacluster_0", "pattern_0", trim_frac=0.08);
mr.plot_pssm("metacluster_0", "pattern_2", trim_frac=0.08);
In [218]:
pattern_name = 'metacluster_0/pattern_0'
In [219]:
pattern = mr.get_pattern(pattern_name)
In [220]:
pattern.plot(rotate_y=0);
In [221]:
def load_scores(h5_file, importance='weighted'):
gt = importance
im = HDF5Reader(h5_file)
im.open()
d = im.load_all()
seq = d['inputs']
tasks = list(d['grads'])
hyp_contrib = {f"{task}": mean(d['grads'][task][gt])
for task in tasks}
contrib = {f"{task}": hyp_contrib[f"{task}"] * seq
for task in tasks}
profile = {f"t/{task}": d['targets']['profile'][task] for task in tasks}
im.close()
return seq, contrib, hyp_contrib, profile
# subset
def get_scores(im, idx, tasks):
hyp_contrib = {task: (im.f[f'/grads/{task}/weighted/0'][idx] +
im.f[f'/grads/{task}/weighted/1'][idx]) / 2
for task in tasks}
seq = im.f['/inputs/'][idx]
contrib = {t: hyp_contrib[t] * seq for t in tasks}
profile = {t: im.f[f'/targets/profile/{t}'][idx] for t in tasks}
return seq, contrib, hyp_contrib, profile
In [11]:
seq, contrib, hyp_contrib, profile = load_scores(model_dir / "grad.valid.h5")
In [12]:
idx = top_max_count(profile['t/Oct4'][:, 400:600], 80)
In [13]:
seq.shape
Out[13]:
In [14]:
example_idx = 11128
In [15]:
viz_dict = filter_tracks({**get_dataset_item(profile, example_idx),
**get_dataset_item(contrib, example_idx)}, xlim=[400, 600])
In [16]:
plot_tracks(viz_dict);
Focus on one specific example¶
In [17]:
t = dict(contrib=contrib['Oct4'][example_idx, 400:600],
hyp_contrib=hyp_contrib['Oct4'][example_idx, 400:600],
seq=seq[example_idx, 400:600])
In [18]:
# target
plot_tracks(dict(contrib=t['contrib'], hyp_contrib=t['hyp_contrib']));
In [19]:
from basepair.modisco.results import *
In [20]:
pattern.trim_seq_ic(0.08).plot(['seq', 'contrib/Oct4', 'hyp_contrib/Oct4'],
rotate_y=0, height=1, letter_width=0.2);
Get the instances¶
In [128]:
pattern_names = [
("Oct4-Sox2", "metacluster_0/pattern_0"),
("Errb", "metacluster_0/pattern_1"),
("Sox2", "metacluster_0/pattern_2"),
("Nanog", "metacluster_0/pattern_3"),
("Klf4", "metacluster_2/pattern_0"),
]
In [223]:
cached = True
In [130]:
cache_dir = Path(Path(ddir) / "cache/dev/")
cache_dir.mkdir(parents=True, exist_ok=True)
cache_file = cache_dir / "hit-scoring.hdf5"
In [131]:
dfp = pd.read_hdf(cache_file, "/scores")
In [132]:
!du -sh {cache_file}
In [133]:
#dfp.to_hdf(cache_file, "/scores")
In [224]:
if not cached:
dfl = []
match_l, importance_l, seq_match_l = [], [], []
for tf, pattern_name in pattern_names:
pattern = mr.get_pattern(pattern_name).trim_seq_ic(0.08)
match, importance = pattern.scan_importance(contrib, hyp_contrib, tasks,
n_jobs=20, verbose=True)
seq_match = pattern.scan_seq(seq, n_jobs=20, verbose=True)
dfm = pattern.get_instances(tasks, match, importance, seq_match, fdr=0.01, verbose=True)
dfm['tf'] = tf
dfl.append(dfm)
match_l.append(match)
importance_l.append(importance)
seq_match_l.append(seq_match)
dfp = pd.concat(dfl)
# cache
cached = True
dfp.to_hdf(cache_file, "/scores")
dfp['idx'] = np.arange(len(dfp))
dfp = dfp[dfp.score_seq_match > 0]
# get also the match scores
match = np.stack(match_l, axis=-1)
importance = np.stack(importance_l, axis=-1)
seq_match = np.stack(seq_match_l, axis=-1)
In [29]:
dfp.head()
Out[29]:
In [247]:
# convert a wide format to a long format
def melt_prefix(dfp, var='imp'):
prefix = var + "/"
dfp_imp = dfp[['idx'] + [prefix + t for t in tasks]].melt(id_vars='idx', value_name=var)
dfp_imp['task'] = dfp_imp.variable.str.replace(prefix, "")
del dfp_imp['variable']
return dfp_imp
dfpm = melt_prefix(dfp, 'imp').merge(melt_prefix(dfp, 'match'), on=['idx', 'task'])
dfpl = dfp[[c for c in dfp.columns
if not c.startswith("imp/")
and not c.startswith("match/")]].merge(dfpm, on='idx')
dfpl['log_imp'] = np.log(1+dfpl.imp)
dfp['log_imp_weighted'] = np.log(1+dfp['imp_weighted'])
dfp['log_match_weighted'] = np.log(dfp['match_weighted'])
In [248]:
# partition the scores into 4 categories (according to the median)
dfp_medians = pd.DataFrame({'log_match_weighted': dfp.groupby("tf").log_match_weighted.median(),
'log_imp_weighted': dfp.groupby("tf").log_imp_weighted.median()}).reset_index()
dfp
log_imp_weighted_median = dfp.tf.map(dfp.groupby("tf").log_imp_weighted.median())
log_match_weighted_median = dfp.tf.map(dfp.groupby("tf").log_match_weighted.median())
dfp['imp_cat'] = (dfp.log_imp_weighted > log_imp_weighted_median).map({False: 'low', True: 'high'})
dfp['match_cat'] = (dfp.log_match_weighted > log_match_weighted_median).map({False: 'low', True: 'high'})
TODO¶
- add percentiles
In [136]:
# dfp['match_max_task'] = dfp.match_max__task
# dfpl['match_max_task'] = dfpl.match_max__task
In [137]:
dfpl.head()
Out[137]:
TODO¶
- sample a few points from the scatterplot: importance vs match to see what you are getting
- try the same analysis for a few factors
- try to find the right cutoff based on the number of motifs you expect
- eyeball a few regions (higlight the found seqlets in there)
- before doing the search, merge all teh metaclusters (ignore the importance score magnitude)
In [249]:
from plotnine import *
import plotnine
In [250]:
dfpl.shape
Out[250]:
In [251]:
dfpl_subset = dfpl.sample(10000) # focus on the subset of the data
dfp_subset = dfp.sample(10000) # focus on the subset of the data
In [252]:
dfpl_subset.head()
Out[252]:
In [231]:
import warnings
warnings.filterwarnings("ignore")
In [143]:
plotnine.options.figure_size = (12, 7)
ggplot(aes(x="match", y="log_imp"), dfpl_subset) + geom_point(alpha=0.1) + \
facet_grid("tf~task", labeller='label_both') + theme_bw()
Out[143]:
Match score with highest importance across tasks¶
In [145]:
plotnine.options.figure_size = (12, 7)
ggplot(aes(x="match_weighted", y="log_imp"), dfpl_subset) + geom_point(alpha=0.1) + \
facet_grid("tf~task", labeller='label_both') + theme_bw()
Out[145]:
Match vs imp¶
In [242]:
plotnine.options.figure_size = (12, 3)
ggplot(aes(x="log_match_weighted", y="log_imp_weighted"), dfp_subset) + geom_point(alpha=0.1) + \
facet_grid(".~tf", labeller='label_both') + theme_bw()
Out[242]:
Seq vs match¶
In [243]:
plotnine.options.figure_size = (12, 3)
ggplot(aes(x="score_seq_match", y="log_match_weighted"), dfp_subset) + geom_point(alpha=0.1) + \
facet_grid(".~tf", labeller='label_both') + theme_bw()
Out[243]:
Seq vs importance¶
In [245]:
plotnine.options.figure_size = (12, 3)
ggplot(aes(x="score_seq_match", y="log_imp_weighted"), dfp_subset) + geom_point(alpha=0.1) + \
facet_grid(".~tf", labeller='label_both') + theme_bw()
Out[245]:
Seq vs match¶
In [234]:
ggplot(aes(x='log_match_weighted'), dfp_subset) + geom_histogram() + facet_grid(".~tf", labeller='label_both') + theme_bw()
Out[234]:
Non-PWM filtered¶
In [240]:
plotnine.options.figure_size = (12, 3)
ggplot(aes(x="log_match_weighted", y="log_imp_weighted"), dfp_subset) + geom_point(alpha=0.1) + \
geom_vline(aes(xintercept='log_match_weighted'), data=dfp_medians, alpha=1, color='orange', linetype='dashed') + \
geom_hline(aes(yintercept='log_imp_weighted'), data=dfp_medians, alpha=1, color='orange', linetype='dashed') + \
facet_grid(".~tf", labeller='label_both') + theme_bw()
Out[240]:
PWM filtered¶
In [253]:
plotnine.options.figure_size = (12, 3)
ggplot(aes(x="log_match_weighted", y="log_imp_weighted"), dfp_subset) + geom_point(alpha=0.1) + \
geom_vline(aes(xintercept='log_match_weighted'), data=dfp_medians, alpha=1, color='orange', linetype='dashed') + \
geom_hline(aes(yintercept='log_imp_weighted'), data=dfp_medians, alpha=1, color='orange', linetype='dashed') + \
facet_grid(".~tf", labeller='label_both') + theme_bw()
Out[253]:
In [ ]:
dfp_subset.log_imp.weighted.median()
In [66]:
np.log(1+dfp.imp_max).plot.hist(200);
#plt.xlim([0, 20])
plt.xlabel("Maximum importance plot");
In [44]:
plotnine.options.figure_size = (7, 7)
ggplot(aes(x='tf', color='task', y='match'), dfpl) + geom_boxplot() + theme_bw()
Out[44]:
In [45]:
plotnine.options.figure_size = (7, 7)
ggplot(aes(x='tf', color='task', y='log_imp'), dfpl) + geom_boxplot() + theme_bw()
Out[45]:
In [46]:
dfpl.groupby(['tf', 'task']).imp.mean().unstack().plot.bar()
plt.ylabel("Average importance");
In [47]:
dfpl.groupby(['tf', 'task']).match.mean().unstack().plot.bar()
plt.title("Average match per task");
plt.ylabel("Average Match");
In [48]:
dfp.groupby("tf").imp_max_task.value_counts().unstack().plot.bar();
plt.ylabel("Frequency")
plt.title("Number of times a particular task has max importance ");
In [49]:
dfp.groupby("tf").match_max_task.value_counts().unstack().plot.bar();
plt.ylabel("Frequency")
plt.title("Number of times a particular task has max match ");
In [50]:
pattern_names = [
("Oct4-Sox2", "metacluster_0/pattern_0"),
("Errb", "metacluster_0/pattern_1"),
("Sox2", "metacluster_0/pattern_2"),
("Nanog", "metacluster_0/pattern_3"),
("Klf4", "metacluster_2/pattern_0"),
]
In [51]:
pattern.plot("contrib");
In [52]:
tasks
task_imp = np.array([np.abs(pattern.contrib[t]).mean() for t in tasks])
task_imp = task_imp / task_imp.sum()
In [53]:
pattern = mr.get_pattern(pattern.name).trim_seq_ic(0.08)
In [40]:
dfpl_subset[(dfpl_subset.task == 'Klf4') & (dfpl_subset.tf == 'Oct4-Sox2')& (dfpl_subset.imp > 40)]
Out[40]:
In [38]:
viz_dict = filter_tracks({**get_dataset_item(profile, example_idx),
**get_dataset_item(contrib, example_idx)}, xlim=[400, 600])
In [41]:
example_idx = 1288
In [67]:
dfp_subset_idx = dfp[(dfp.example_idx == example_idx) & (dfp.pattern_start > 400) & (dfp.pattern_end < 600)]
dfp_subset_idx
Out[67]:
In [69]:
np.log(1+dfp_subset_idx.imp_max)
Out[69]:
In [70]:
seqlets = [Seqlet(row.example_idx, row.pattern_start-400, row.pattern_end-400, row.tf, row.strand)
for i,row in dfp_subset_idx.iterrows()]
In [72]:
# TODO - is this thing normalized by importance?
In [71]:
plot_tracks(filter_tracks({**get_dataset_item(profile, example_idx),
**get_dataset_item(contrib, example_idx)}, xlim=[400, 600]),
seqlets, rotate_y=0, legend=True);
Number of motifs per example¶
In [102]:
dfp[(np.log(1+dfp.imp_max) > 2) & (dfp.pattern_start > 400) & (dfp.pattern_end < 600)].groupby('example_idx').size().plot.hist()
Out[102]:
In [256]:
def prefix_dict(d, prefix):
return {prefix + d: v for d,v in d.items()}
Example of strong reliance on hypothetical scores¶
In [282]:
example_idx = dfp.example_idx.sample(1).iloc[0]
print(example_idx)
In [283]:
dfp_subset_idx = dfp[#(np.log(1+dfp.imp_max) > 2) &
(dfp.imp_cat == 'high') &
(dfp.match_cat == 'high') &
#(dfp.tf == "Errb") &
(dfp.example_idx == example_idx) &
(dfp.pattern_start > 400) &
(dfp.pattern_end < 600)]
seqlets = [Seqlet(row.example_idx, row.pattern_start-400, row.pattern_end-400, row.tf, row.strand)
for i,row in dfp_subset_idx.iterrows()]
plot_tracks(filter_tracks({**get_dataset_item(contrib, example_idx),
**get_dataset_item(prefix_dict(hyp_contrib, "h/"), example_idx)}, xlim=[400, 600]),
seqlets, rotate_y=0, legend=True);
In [287]:
dfp_subset_idx[['tf', 'pattern_center', 'strand', 'match_weighted', 'imp_weighted', 'score_seq_match']].sort_values('pattern_center')
Out[287]:
In [285]:
for tf, pattern_name in pattern_names:
p = mr.get_pattern(pattern_name).trim_seq_ic(0.08)
p.name = "+" + tf
p.plot(['seq'] + [f'contrib/{t}' for t in tasks], rotate_y=0)
p.name = "-" + tf
p.rc().plot(['seq'] + [f'contrib/{t}' for t in tasks], rotate_y=0)
Example of strong reliance on hypothetical scores¶
In [254]:
example_idx = dfp.example_idx.sample(1).iloc[0]
print(example_idx)
In [272]:
dfp_subset_idx = dfp[#(np.log(1+dfp.imp_max) > 2) &
(dfp.imp_cat == 'high') &
(dfp.match_cat == 'high') &
#(dfp.tf == "Errb") &
(dfp.example_idx == example_idx) &
(dfp.pattern_start > 400) &
(dfp.pattern_end < 600)]
seqlets = [Seqlet(row.example_idx, row.pattern_start-400, row.pattern_end-400, row.tf, row.strand)
for i,row in dfp_subset_idx.iterrows()]
plot_tracks(filter_tracks({**get_dataset_item(contrib, example_idx),
**get_dataset_item(prefix_dict(hyp_contrib, "h/"), example_idx)}, xlim=[400, 600]),
seqlets, rotate_y=0, legend=True);
In [275]:
dfp_subset_idx[['tf', 'pattern_center', 'strand', 'match_weighted', 'imp_weighted', 'score_seq_match']].sort_values('pattern_center')
Out[275]:
In [268]:
dfp_subset_idx = dfp[#(np.log(1+dfp.imp_max) > 2) &
(dfp.imp_cat == 'high') &
(dfp.match_cat == 'high') &
(dfp.tf == "Sox2") &
(dfp.example_idx == example_idx) &
(dfp.pattern_start > 400) &
(dfp.pattern_end < 600)]
seqlets = [Seqlet(row.example_idx, row.pattern_start-400, row.pattern_end-400, row.tf, row.strand)
for i,row in dfp_subset_idx.iterrows()]
plot_tracks(filter_tracks({**get_dataset_item(contrib, example_idx),
**get_dataset_item(prefix_dict(hyp_contrib, "h/"), example_idx)}, xlim=[400, 600]),
seqlets, rotate_y=0, legend=True);
Conclusion¶
- hypothetical contribution scores can be confusing to inspect
In [269]:
dfp_subset_idx = dfp[#(np.log(1+dfp.imp_max) > 2) &
(dfp.imp_cat == 'high') &
(dfp.match_cat == 'high') &
(dfp.tf == "Oct4-Sox2") &
(dfp.example_idx == example_idx) &
(dfp.pattern_start > 400) &
(dfp.pattern_end < 600)]
seqlets = [Seqlet(row.example_idx, row.pattern_start-400, row.pattern_end-400, row.tf, row.strand)
for i,row in dfp_subset_idx.iterrows()]
plot_tracks(filter_tracks({**get_dataset_item(contrib, example_idx),
**get_dataset_item(prefix_dict(hyp_contrib, "h/"), example_idx)}, xlim=[400, 600]),
seqlets, rotate_y=0, legend=True);
In [270]:
dfp_subset_idx = dfp[#(np.log(1+dfp.imp_max) > 2) &
(dfp.imp_cat == 'high') &
(dfp.match_cat == 'high') &
(dfp.tf == "Nanog") &
(dfp.example_idx == example_idx) &
(dfp.pattern_start > 400) &
(dfp.pattern_end < 600)]
seqlets = [Seqlet(row.example_idx, row.pattern_start-400, row.pattern_end-400, row.tf, row.strand)
for i,row in dfp_subset_idx.iterrows()]
In [279]:
plot_tracks(filter_tracks({**get_dataset_item(contrib, example_idx),
**get_dataset_item(profile, example_idx)}, xlim=[400, 600]),
rotate_y=0, legend=True);
In [270]:
plot_tracks(filter_tracks({**get_dataset_item(contrib, example_idx),
**get_dataset_item(prefix_dict(hyp_contrib, "h/"), example_idx)}, xlim=[400, 600]),
seqlets, rotate_y=0, legend=True);
