Goal¶
- analyze instance locations
Links¶
Conclusions / questions¶
- does Klf4 compete with others for the site?
- there are many sites where the model says that if Klf4 binds, all other binding will be depleted?
- oct4-sox2 is realatively in-frequent in the dataset. Is this
- 3/2, 4/1 -> Nanog motif in 2 different context. One occuring in the middle while the other occuring on the side
- m1/p0 and m1/p1 are frequently together
- instance location assignment could be improved
- how to distinquish co-binding vs nearby binding?
- idea: have a look at the specific profile
Ideas¶
Shall we just run modisco for each peak-set separately?
- conditioned that we chipped the factor, what are other factors in the vicinity doing?
Improvements¶
- Aggregation plot
- plot negative strand as -
- always start at 0
- ylim = [0, max(abs(pos), abs(neg))]
calculate entropy for each region and the final aggregated distribution
- show a histogram
- or a scatterplot with the total counts in the particular region or in the total region
threshold all the computations w.r.t. the total importance score in the region
[ ] compare the PWM, and the aggregated profile across all pattern locations (not just those identified by clustering with modisco)
- can you improve the results by filtering either by seqlet_score or percnormscore?
[ ] scatterplot the spread or bi-modality of the distribution witin the sequence vs. seqlet_score
- are high seqlet scores mostly in the center?
In [5]:
from tqdm import tqdm
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
from kipoi.readers import HDF5Reader
Load the data¶
Modisco¶
In [6]:
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"
!ls {model_dir}
!ls {modisco_dir} # modisco dir
In [7]:
# Modisco
mr = ModiscoResult(modisco_dir / "modisco.h5")
mr.open()
tasks = mr.tasks()
In [9]:
ds = DataSpec.load(model_dir / "dataspec.yaml")
In [69]:
d = HDF5Reader(model_dir / "grad.all.h5")
d.open()
id_hash = pd.DataFrame({"peak_id": d.f['/metadata/interval_from_task'][:], "example_idx": np.arange(d.f['/metadata/interval_from_task'].shape[0])})
Labelled instances¶
In [70]:
# load the instances data frame
df = pd.read_csv(modisco_dir / "all.instances.tsv", sep='\t')
df['pattern_id'] = "metacluster_" + df.metacluster.astype(str) + "/pattern_" + df.pattern.astype(str)
df['seqlet_center'] = (df.seqlet_start + df.seqlet_end) / 2
patterns = df.pattern_id.unique().tolist()
pattern_pssms = {pattern: mr.get_pssm(*pattern.split("/")) for pattern in patterns}
append_pattern_loc(df, pattern_pssms, trim_frac=0.08)
df['pattern_center'] = (df.pattern_start + df.pattern_end) / 2
df = df.merge(id_hash, on="example_idx")
In [71]:
# Use thresholded values from now on
df_unfiltered = df
df = df[df.percnormed_score > .2]
In [178]:
ls {outdir}
In [17]:
seqlets = read_pkl(modisco_dir / "all.seqlets.pkl")
In [ ]:
assert len(df) == len(seqlets)
# total number of mapped seqlets
len(seqlets)
Grads and inputs¶
In [11]:
n_examples = d.f['/metadata/range/chr'].shape[0]
In [12]:
# total number of peaks / examples
n_examples
Out[12]:
In [92]:
# sanity check the raw peak file
for task, ts in ds.task_specs.items():
n_peaks = !zcat {ts.peaks} | wc -l
print(f"{task}: {n_peaks[0]}")
In [ ]:
df.head()
In [ ]:
# almost every peak has a seqlet
len(df.example_id.unique()) / n_examples
Write the interval_from_task vectors¶
In [175]:
d.f['/metadata/interval_from_task'][:4]
Out[175]:
In [177]:
outdir = model_dir / "clustering/peak-set/"
outdir.mkdir(parents=True, exist_ok=True)
metadata_task = d.f['/metadata/interval_from_task'][:]
for task in tasks:
mask = metadata_task == task
np.save(outdir / f"{task}.npy", mask)
Stats¶
Number of seqlets per example¶
In [21]:
from plotnine import *
import plotnine
In [22]:
plotnine.options.figure_size = (10, 1)
ggplot(aes(x="size"), df.groupby("example_id").pattern.agg(['size']).reset_index()) + \
geom_bar() + theme_bw() + \
xlab("Number of seqlets per example")
Out[22]:
Pattern frequency¶
In [23]:
# 20 most frequent patterns
fig = plt.figure(figsize=(3, 8))
df.pattern_id.value_counts().head(30).iloc[::-1].plot(kind='barh');
In [24]:
len(df.pattern_id.value_counts())
Out[24]:
Functions¶
In [88]:
from basepair.vdom_viz import fig2vdom
from vdom.helpers import *
def plot_n_per_example_ax(dfs, ax):
dfs.groupby("example_id").size().value_counts().plot("barh", ax=ax)
ax.set_title("# per example")
def plot_col_hist_ax(dfs, ax, col='score', xlim=None):
dfs[col].plot.hist(bins=100, ax=ax)
#ax.set_xlabel("# per example")
if xlim is not None:
ax.set_xlim(xlim)
ax.set_title(col)
def plot_all_dist(df, pattern_id, figsize=(12,2), add_title="", **kwargs):
"""Merge all plots into a single one
"""
dfs = df[df.pattern_id == pattern_id]
fig, ax = plt.subplots(1, 5, figsize=figsize)
plot_col_hist_ax(dfs, ax[0], "seqlet_center", xlim=[0, 1000])
ax[0].set_title(f"seqlet_center {add_title}")
plot_n_per_example_ax(dfs, ax[1])
plot_col_hist_ax(dfs, ax[2], "percnormed_score")
plot_col_hist_ax(dfs, ax[3], "score")
plot_col_hist_ax(dfs, ax[4], "seqlet_score")
plt.tight_layout()
return fig2vdom(fig, **kwargs)
import matplotlib.pyplot as plt
from concise.utils.plot import seqlogo_fig, seqlogo
from basepair.math import mean
from basepair.plots import draw_box
# helper functions for plotting the results
def get_scores(d, example_idx, task, grad_type="weighted"):
return mean([d.f[f'/grads/{task}/{grad_type}/0'][example_idx],
d.f[f'/grads/{task}/{grad_type}/1'][example_idx]])
def get_sequence(d, example_idx):
return d.f['/inputs'][example_idx]
def relevant_patterns(df, example_idx, xlim):
patterns = []
for i,row in df[df.example_idx == example_idx].sort_values("pattern_start").iterrows():
xmin = row.pattern_start - xlim[0] + 0.5
xmax = row.pattern_end - xlim[0]+ 0.5
if xmin < 0:
continue
if xmax > xlim[1] - xlim[0]:
continue
patterns.append((row.pattern_id, row.revcomp))
return patterns
def draw_patterns(df, example_idx, xlim, ax, add_label=True):
for i,row in df[df.example_idx == example_idx].iterrows():
xmin = row.pattern_start - xlim[0] + 0.5
xmax = row.pattern_end - xlim[0]+ 0.5
if xmin < 0:
continue
if xmax > xlim[1] - xlim[0]:
continue
draw_box(xmin, xmax, ax, 'g') # trimmed pattern location
if add_label:
y = ax.get_ylim()[1] + (ax.get_ylim()[1] - ax.get_ylim()[0]) * 0.1
if row.revcomp:
prefix = "-"
else:
prefix = ""
ax.text(row.pattern_start - xlim[0] + 0.5, y,
s=prefix + f"mc_{row.metacluster}/p_{row.pattern}")
def plot_region(df, d, tasks, example_idx, xlim=[400,600]):
fig, axes = plt.subplots(len(tasks), 1, figsize=(20, 2* len(tasks)), sharex=True)
seq = get_sequence(d, example_idx)
for i, task in enumerate(tasks):
ax = axes[i]
seqlogo((get_scores(d, example_idx, task) * seq)[xlim[0]:xlim[1]], ax=ax)
draw_patterns(df, example_idx, xlim, ax, add_label=i==0)
ax.set_ylabel(task)
ax.set_xticks(list(range(0, xlim[1] - xlim[0] + 1, 5)));
def rc_prefix(rc):
if rc:
return "-"
else:
return ""
def vdom_patterns(mr, patterns_rc_vec, trim_frac=0.08):
"""Get vdom patterns in a single row
Args:
mr: ModiscoResults object
patterns_rc_vec: list of tuples (pattern_id, rc). Example: ("metacluster_0/pattern_1", True)
trim_frac: how much to trim the pssm
"""
return p([fig2vdom(mr.plot_pssm(*pattern_idx.split("/"), rc=rc, trim_frac=trim_frac, title=rc_prefix(rc) + pattern_idx,
letter_width=0.35, height=1.7), height=75)
for pattern_idx, rc in patterns_rc_vec])
def dist_pairwise_hist(df, a, b, ax=None):
if a != b:
dfc_ab = df[df.pattern_id.isin([a,b])]
counts = dfc_ab.groupby("example_idx").pattern_id.value_counts().unstack(fill_value=0)
relevant_idx = counts.index[(counts==1).all(axis=1)]
distances = dfc_ab[dfc_ab.example_idx.isin(relevant_idx)].set_index(['example_idx', 'pattern_id']).pattern_center.unstack()
diff = distances[a] - distances[b]
title = f"{a} - {b}"
if a == b:
examples = df[df.pattern_id == a].groupby("example_idx").filter(lambda x: len(x) == 2)
diff = examples.groupby("example_idx").apply(lambda x: x.pattern_center.iloc[1] - x.pattern_center.iloc[0])
title = f"{a}"
if ax is None:
fig,ax = plt.subplots(figsize=(15,2))
ax.hist(diff, bins=100);
ax.set_title(title);
In [60]:
df.head()
Out[60]:
The plot bellow shows the following:
- positional distribution within the 1kb window
- number of seqlets per example
- PSSM
In [96]:
## fig_height = 100
div(*[div(span(str(i), ". ", row['index'] + f" (n={row.pattern_id})",
fig2vdom(mr.plot_pssm(*row['index'].split("/"), height=0.5, letter_width=0.15, trim_frac=0.08))),
br(),
*[plot_all_dist(df[df.peak_id == task], row['index'], add_title=task,
figsize=(16,1.8), height=100) for task in tasks],
)
for i, row in df.pattern_id.value_counts().reset_index().iterrows()])
Out[96]:
0. metacluster_1/pattern_0 (n=59068)
1. metacluster_2/pattern_0 (n=52170)
2. metacluster_0/pattern_0 (n=40564)
3. metacluster_0/pattern_1 (n=35094)
4. metacluster_2/pattern_3 (n=31215)
5. metacluster_0/pattern_2 (n=25087)
6. metacluster_1/pattern_1 (n=20328)
7. metacluster_3/pattern_2 (n=19276)
8. metacluster_4/pattern_1 (n=16044)
9. metacluster_1/pattern_2 (n=15511)
10. metacluster_3/pattern_6 (n=13137)
11. metacluster_4/pattern_0 (n=12320)
12. metacluster_4/pattern_2 (n=11480)
13. metacluster_3/pattern_1 (n=10655)
14. metacluster_4/pattern_3 (n=10154)
15. metacluster_7/pattern_0 (n=9145)
16. metacluster_6/pattern_4 (n=6858)
17. metacluster_1/pattern_4 (n=6775)
18. metacluster_10/pattern_0 (n=6037)
19. metacluster_1/pattern_5 (n=5770)
20. metacluster_3/pattern_0 (n=4978)
21. metacluster_3/pattern_5 (n=4455)
22. metacluster_28/pattern_0 (n=4423)
23. metacluster_0/pattern_4 (n=4301)
24. metacluster_0/pattern_3 (n=4243)
25. metacluster_15/pattern_0 (n=3882)
26. metacluster_0/pattern_8 (n=3825)
27. metacluster_21/pattern_0 (n=3779)
28. metacluster_6/pattern_2 (n=3603)
29. metacluster_2/pattern_1 (n=3526)
30. metacluster_2/pattern_4 (n=3384)
31. metacluster_3/pattern_7 (n=3371)
32. metacluster_22/pattern_0 (n=3371)
33. metacluster_31/pattern_0 (n=2979)
34. metacluster_6/pattern_1 (n=2965)
35. metacluster_7/pattern_2 (n=2949)
36. metacluster_1/pattern_3 (n=2770)
37. metacluster_7/pattern_1 (n=2767)
38. metacluster_0/pattern_6 (n=2755)
39. metacluster_0/pattern_5 (n=2597)
40. metacluster_9/pattern_1 (n=2592)
41. metacluster_7/pattern_3 (n=2567)
42. metacluster_3/pattern_4 (n=2468)
43. metacluster_9/pattern_0 (n=2197)
44. metacluster_9/pattern_2 (n=1940)
45. metacluster_25/pattern_0 (n=1922)
46. metacluster_6/pattern_3 (n=1833)
47. metacluster_0/pattern_9 (n=1778)
48. metacluster_7/pattern_4 (n=1673)
49. metacluster_18/pattern_0 (n=1631)
50. metacluster_13/pattern_0 (n=1629)
51. metacluster_10/pattern_1 (n=1527)
52. metacluster_0/pattern_11 (n=1519)
53. metacluster_6/pattern_6 (n=1511)
54. metacluster_6/pattern_0 (n=1477)
55. metacluster_0/pattern_10 (n=1369)
56. metacluster_3/pattern_3 (n=1108)
57. metacluster_2/pattern_2 (n=1062)
58. metacluster_1/pattern_6 (n=974)
59. metacluster_10/pattern_2 (n=886)
60. metacluster_23/pattern_1 (n=857)
61. metacluster_9/pattern_3 (n=791)
62. metacluster_18/pattern_2 (n=789)
63. metacluster_10/pattern_4 (n=688)
64. metacluster_18/pattern_1 (n=570)
65. metacluster_6/pattern_5 (n=438)
66. metacluster_23/pattern_0 (n=431)
67. metacluster_10/pattern_3 (n=423)
68. metacluster_9/pattern_4 (n=381)
69. metacluster_0/pattern_7 (n=239)
Analysis¶
- metacluster_2/pattern_0 = GGTGTGG = Klf4 motif http://meme-suite.org/opal-jobs/appTOMTOM_5.0.1_15351119239261901384692/tomtom.html
- metacluster_1/pattern_0 = Nanog
- metacluster_2/pattern_3 = Klf4
- metacluster_0/pattern_0 = Oct4-Sox2
- metacluster_0/pattern_2 = Nanog
Co-occurence stats¶
Compute the co-occurence stats.
Filter¶
In [28]:
from upsetplot import plot as upsetplt
In [14]:
all_patterns = mr.patterns()
In [15]:
# restrict to the central 500bp
df_center = df[(df.pattern_center > 250) & (df.pattern_center < 750)]
In [17]:
len(df)
Out[17]:
In [18]:
len(df_center)
Out[18]:
In [19]:
dfp = pd.DataFrame.from_dict([dict(pd.Series(all_patterns,
index=all_patterns).
isin(group.pattern_id.unique()))
for name, group in tqdm(df_center.groupby("example_id"))])
In [97]:
# make the upsetplot
data_dict = dfp.groupby(all_patterns).size()
dfc = data_dict.reset_index().rename(columns={0:"n"})
dfcp = dfc.sort_values("n", ascending=False).head(40) # take only 40 most frequent combinations
col_incl = dfcp.iloc[:, :-1].max()
dfcp_plot = dfcp.iloc[:, col_incl.tolist() + [True]]
dfcp_plot = dfcp_plot.set_index(list(dfcp_plot.columns[:-1]))
Plot¶
In [36]:
r = upsetplt(dfcp_plot.n)
In [98]:
# note interesting pairs to plot pairwise distances
pairs = [
["metacluster_2/pattern_0", "metacluster_1/pattern_0"],
["metacluster_2/pattern_0", "metacluster_0/pattern_0"],
["metacluster_2/pattern_0", "metacluster_28/pattern_0"],
["metacluster_1/pattern_0", "metacluster_2/pattern_3"],
["metacluster_0/pattern_0", "metacluster_0/pattern_1"],
["metacluster_0/pattern_0", "metacluster_0/pattern_2"],
["metacluster_0/pattern_0", "metacluster_3/pattern_2"],
["metacluster_1/pattern_0", "metacluster_1/pattern_1"],
["metacluster_1/pattern_2", "metacluster_0/pattern_0"],
["metacluster_1/pattern_2", "metacluster_0/pattern_1"],
]
In [37]:
# Visualie motifs for most frequent patterns from above
dfcp_plot_melt = pd.melt(dfcp_plot.reset_index().assign(idx=np.arange(len(dfcp_plot))), id_vars=["n","idx"])
patterns = dfcp_plot_melt.query('value').groupby("variable").n.sum().sort_values(ascending=False).index.tolist()
display(vdom_patterns(mr, zip(patterns, [False]*len(patterns)), trim_frac=0.08))
Pairwise distances with self¶
- restrict to sites with exactly 2 instances of the same pattern and plot a pairwise distance
In [54]:
nrow = int(np.ceil(len(patterns)/2))
fig, axes = plt.subplots(nrow,2, figsize=(20, nrow * 1.5), sharex=True)
for pattern, ax in tqdm(zip(patterns, axes.flat)):
dist_pairwise_hist(df[~df[['example_id', 'pattern_center']].duplicated()], pattern, pattern, ax)
ax.set_xlim([-1000, 1000])
plt.tight_layout()
Pairwise distances with a pair¶
- restrict to sites with exactly 1 instance for each motif and plot their pairwise distance
In [93]:
## TODO statify those plots by the peak_id
In [95]:
nrow = int(np.ceil(len(pairs)/2))
fig, axes = plt.subplots(nrow, 2, figsize=(20, nrow * 1.5), sharex=True)
for pair, ax in tqdm(zip(pairs, axes.flat)):
dist_pairwise_hist(df[~df[['example_id', 'pattern_center']].duplicated()], pair[0], pair[1], ax)
ax.set_xlim([-1000, 1000])
plt.tight_layout()
Eyeballing¶
- plots bellow this point are used to check how well the assignment of seqlets to patterns work
In [36]:
example_idx = 12
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, 12, xlim)
In [38]:
example_idx = 29934
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [39]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
In [40]:
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [41]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [601]:
# This example doesn't make any sense to me
In [66]:
df[df.example_idx == 5472]
Out[66]:
In [72]:
sl = seqlets[638799]
In [98]:
s = sl.track_name_to_snippet['Sox2_contrib_scores']
In [99]:
seqlogo_fig(s.fwd, figsize=(10,2));
In [42]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
example_idx = 5472
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [99]:
seqlogo_fig(s.fwd, figsize=(10,2));
In [101]:
mr.plot_pssm("metacluster_9","pattern_0");
In [43]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [44]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [45]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [46]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [47]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [49]:
example_idx = df.example_idx.sample(1).iloc[0]
print(example_idx)
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [30]:
df[df.example_idx == 36709].iloc[0][['example_chr', 'example_start', 'example_end']]
Out[30]:
Comparing models¶
Load old Oct-Sox BPNet¶
In [34]:
from basepair.BPNet import BPNetPredictor
from keras.models import load_model
from basepair.config import create_tf_session
from pybedtools import Interval
In [4]:
create_tf_session(0)
Out[4]:
In [5]:
model = load_model("/users/avsec/workspace/basepair/basepair/../data/processed/chipnexus/exp/models/p-multi-task/seq_mutlitask_filters=32,conv1_kernel_size=21,tconv_kernel_size=25,n_dil_layers=6,use_profile=True,use_counts=True,c_task_weight=10,lr=0.004.5.h5")
preproc = read_pkl("/users/avsec/workspace/basepair/basepair/../data/processed/chipnexus/exp/models/p-multi-task/preproc.pkl")
bp = BPNetPredictor(model,
"/mnt/data/pipeline_genome_data/mm10/mm10_no_alt_analysis_set_ENCODE.fasta",
tasks = ["Sox2", "Oct4"],
preproc=preproc)
bws = {task: {s: f"/users/avsec/workspace/basepair-workflow/data/{task}/5prime.counts.{s}.bw" for s in ['pos', 'neg']} for task in ["Sox2", "Oct4"]}
genome_file = "/mnt/data/pipeline_genome_data/mm10/mm10.chrom.sizes"
Make predictions¶
In [39]:
df[(df.example_chr == "chr17") & (df.example_start < 35503946) & (df.example_end >35503946+201 )]
Out[39]:
Core motif not recognized¶
In [47]:
# Filtered seqlets (percnormed_score > 0.2)
example_idx = 2
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df, example_idx, xlim), trim_frac=0.08))
plot_region(df, d, tasks, example_idx, xlim)
In [50]:
# Unfilterred seqlets
example_idx = 2
xlim=[400,600]
display(vdom_patterns(mr, relevant_patterns(df_unfiltered, example_idx, xlim), trim_frac=0.08))
plot_region(df_unfiltered, d, tasks, example_idx, xlim)
In [37]:
# old bpnet prediction
bp.plot_predict_grad([Interval("chr17", 35503946, 35503946+201)], bws=bws,
profile_grad='weighted',
figsize=(20, 10))
