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 [1]:
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 [2]:
model_dir = Path(f"{ddir}/processed//chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/")
In [3]:
modisco_dir = model_dir / "modisco/by_peak_tasks/weighted/Oct4"
In [4]:
mr = ModiscoResult(modisco_dir / "modisco.h5")
mr.open()
In [5]:
mr.plot_pssm("metacluster_0", "pattern_0", trim_frac=0.08);
mr.plot_pssm("metacluster_0", "pattern_2", trim_frac=0.08);
In [6]:
pattern_name = 'metacluster_0/pattern_0'
In [7]:
pattern = mr.get_pattern(pattern_name)
In [8]:
pattern.plot(rotate_y=0);
In [9]:
im = HDF5Reader(model_dir / "grad.valid.h5")
im.open()
In [10]:
d = im.load_all()
In [11]:
tasks = ['Oct4', 'Sox2', 'Klf4', 'Nanog']
In [12]:
all_hyp_contrib = (im.f['/grads/Oct4/weighted/0'][:] + im.f['/grads/Oct4/weighted/1'][:]) / 2
all_seq = im.f['/inputs/'][:]
all_contrib = all_hyp_contrib * all_seq
all_profile = im.f['/targets/profile/Oct4'][:]
In [13]:
idx = top_max_count(all_profile[:, 400:600], 80)
In [14]:
example_idx = 11128
In [15]:
gt = 'weighted'
hyp_contrib_scores = {f"{task}": mean(d['grads'][task][gt])
for task in tasks}
contrib_scores = {f"{task}": hyp_contrib_scores[f"{task}"] * all_seq
for task in tasks}
targets = {f"t/{task}": d['targets']['profile'][task] for task in tasks}
In [16]:
viz_dict = filter_tracks({**get_dataset_item(targets, example_idx), **get_dataset_item(contrib_scores, example_idx)}, xlim=[400, 600])
In [17]:
viz_dict.keys()
Out[17]:
In [18]:
plot_tracks(viz_dict);
Focus on one specific example¶
In [19]:
contrib = all_contrib[example_idx, 400:600]
hyp_contrib = all_hyp_contrib[example_idx, 400:600]
seq = all_seq[example_idx, 400:600]
t = dict(contrib=contrib,
hyp_contrib=hyp_contrib,
seq=seq)
In [20]:
# target
plot_tracks(dict(contrib=contrib, hyp_contrib=hyp_contrib));
In [21]:
from basepair.modisco.results import *
In [22]:
pattern.trim_seq_ic(0.08).plot(['seq', 'contrib/Oct4', 'hyp_contrib/Oct4'], rotate_y=0, height=2, letter_width=0.5);
In [23]:
pattern = pattern.trim_seq_ic(0.08)
Open questions¶
- get only ~10 motifs (no long-ones) when not using multiple tasks
In [24]:
pattern.scan_importance?
Oct4 match¶
In [25]:
pattern = mr.get_pattern("metacluster_0/pattern_0").trim_seq_ic(0.08)
In [26]:
pattern.plot("seq");
In [28]:
match, importance = pattern.scan_importance(contrib_scores, hyp_contrib_scores, tasks, n_jobs=20)
In [29]:
# optional
seq_match = pattern.scan_seq(all_seq, n_jobs=20)
In [38]:
dfm = pattern.get_instances(match, importance, seq_match, fdr=0.01, verbose=True)
In [39]:
dfm.head()
Out[39]:
Run for all major patterns¶
In [27]:
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 [33]:
dfl = []
for tf, pattern_name in pattern_names:
pattern = mr.get_pattern(pattern_name).trim_seq_ic(0.08)
match, importance = pattern.scan_importance(contrib_scores, hyp_contrib_scores, tasks,
n_jobs=20, verbose=True)
seq_match = pattern.scan_seq(all_seq, n_jobs=20, verbose=True)
dfm = pattern.get_instances(match, importance, seq_match, fdr=0.01, verbose=True)
dfm['tf'] = tf
dfl.append(dfm)
dfp = pd.concat(dfl)
In [41]:
!mkdir -p {ddir}/cache/dev/
In [42]:
!rmdir {ddir}/cache/dev/hit-scoring.feather
In [39]:
dfp.info()
In [46]:
dfp.to_hdf(f"{ddir}/cache/dev/hit-scoring.hdf5", "/scores")
In [49]:
!du -sh {ddir}/cache/dev/hit-scoring.hdf5
In [34]:
dfp.head()
Out[34]:
Eyeball some of the results¶
Dev¶
Convert to pd.DataFrame¶
In [378]:
seq_matchf = seq_match.reshape((-1, 2))
choose_max_strand = np.argmax(matchf[idx_list], axis=-1)
matchfs = matchf[idx_list, choose_max_strand]
importancefs = importancef[idx_list]
seq_matchfs = seq_matchf[idx_list, choose_max_strand]
Sox2 match¶
In [35]:
sox2_pattern = mr.get_pattern("metacluster_0/pattern_2").trim_seq_ic(0.08)
In [36]:
sox2_pattern.plot("seq");
In [37]:
sox2_match, sox2_importance = sox2_pattern.scan_importance(contrib_scores, hyp_contrib_scores, tasks, n_jobs=20)
In [38]:
sox2_seq_match = sox2_pattern.scan_seq(all_seq, n_jobs=20)
Plots¶
In [41]:
match.shape
Out[41]:
In [45]:
import random
In [46]:
idx_list = np.random.randint(0, 19137000, 100000)
In [86]:
plt.hist(match.flatten()[idx_list], 100);
In [138]:
match.shape
Out[138]:
Correlation among tasks¶
In [378]:
matchf = match.reshape((-1, 4, 2)).mean(axis=-2)
importancef = importance.reshape((-1, 4, 2)).mean(axis=-1)
seq_matchf = seq_match.reshape((-1, 2))
choose_max_strand = np.argmax(matchf[idx_list], axis=-1)
matchfs = matchf[idx_list, choose_max_strand]
importancefs = importancef[idx_list]
seq_matchfs = seq_matchf[idx_list, choose_max_strand]
In [356]:
np.sum(matchf.max(axis=-1) > 0.0673)
Out[356]:
In [367]:
thr = fdr_threshold_norm_right(matchf.max(axis=-1), 99, fdr=0.001)
thr
Out[367]:
In [380]:
keep = matchf.max(axis=-1) > thr
In [401]:
thr
Out[401]:
In [383]:
plt.hist(importancef[keep,0].max(axis=-1), 100);
In [394]:
plt.plot(matchf[keep].max(axis=-1),importancef[keep,0], "." , alpha=0.05);
plt.ylabel("Importance")
plt.xlabel("Match");
In [400]:
g = sns.jointplot(matchf[keep].max(axis=-1), importancef[keep,0], alpha=.01);
g.set_axis_labels("Match", "Importance");
TODO¶
- sample a few points from this plot to get a sense what you are picking
In [387]:
np.sum(importancef[keep,0] > 0.2)
Out[387]:
In [385]:
plt.hist(importancef[keep,0], 100);
In [384]:
plt.hist(importancef.max(axis=-1), 100);
In [368]:
np.sum(matchf.max(axis=-1) > thr)
Out[368]:
In [369]:
np.mean(matchf.max(axis=-1) > thr)
Out[369]:
In [84]:
sox2_matchf = sox2_match.reshape((-1, 4, 2)).mean(axis=-2)
sox2_importancef = sox2_importance.reshape((-1, 4, 2))
sox2_seq_matchf = sox2_seq_match.reshape((-1, 2))
sox2_choose_max_strand = np.argmax(sox2_matchf[idx_list], axis=-1)
sox2_matchfs = sox2_matchf[idx_list, sox2_choose_max_strand]
sox2_importancefs = sox2_importancef[idx_list, :, sox2_choose_max_strand]
sox2_seq_matchfs = sox2_seq_matchf[idx_list, sox2_choose_max_strand]
In [42]:
seq_matchf.shape
Out[42]:
In [43]:
matchf.shape
Out[43]:
Correlation between tasks¶
In [49]:
import seaborn as sns
In [58]:
sns.pairplot(pd.DataFrame(matchf[idx_list[:10000], :, 0], columns=tasks))
Out[58]:
In [87]:
sns.pairplot(pd.DataFrame(importancefs, columns=tasks))
Out[87]:
Task¶
- for now, let's just try to get solid motif instances
- importance scores sufficiently high
- match sufficiently high
In [90]:
plt.hist(matchfs, 100);
In [97]:
plt.hist(sox2_matchfs, 100);
In [99]:
plt.hist(importancefs.max(axis=-1), 100);
In [100]:
plt.hist(sox2_importancefs.max(axis=-1), 100);
In [334]:
from basepair.stats import fdr_threshold_norm_right
In [335]:
fdr_threshold_norm_right(matchfs)
Out[335]:
In [92]:
g = sns.jointplot(matchfs, importancefs[:, 0], alpha=.1);
g.set_axis_labels("Match", "Importance");
In [94]:
g = sns.jointplot(matchfs, importancefs.max(axis=-1), alpha=.1);
g.set_axis_labels("Match", "Importance");
In [322]:
fdr_threshold_norm_right(sox2_matchfs)
Out[322]:
In [95]:
g = sns.jointplot(sox2_matchfs, sox2_importancefs[:, 0], alpha=.1);
g.set_axis_labels("Match", "Importance");
In [96]:
g = sns.jointplot(sox2_matchfs, sox2_importancefs.max(axis=-1), alpha=.1);
g.set_axis_labels("Match", "Importance");
In [346]:
g = sns.jointplot(matchfs, seq_matchfs, alpha=0.05);
g.set_axis_labels("Match", "PWM-scan");
In [344]:
fdr_threshold_norm_right(seq_matchfs)
Out[344]:
In [517]:
stats.probplot(seq_matchfs, dist="norm", plot=plt);
plt.title("asd");
Out[517]:
Interstingly, the PWM is a very poor predictor
In [78]:
g = sns.jointplot(matchf[idx_list, :,choose_max].mean(axis=-1), seq_matchf[idx_list,choose_max], alpha=0.05);
g.set_axis_labels("Match", "PWM-scan");
Figure out the right thresholding strategy¶
QQ-plots¶
In [101]:
import scipy.stats as stats
stats.probplot(matchfs, dist="norm", plot=plt);
In [102]:
stats.probplot(sox2_matchfs, dist="norm", plot=plt);
- Experiemented with fitting gaussian mixture models - didn't work at all
In [310]:
np.sum(matchfs > threshold(matchfs) )
Out[310]:
In [336]:
from scipy.stats import norm
In [295]:
norm.ppf(0.8, loc=loc, scale=scale)
Out[295]:
In [265]:
loc, scale = norm.fit(matchfs)
loc, scale
Out[265]:
In [327]:
upper = np.percentile(matchfs,99)
loc, scale = norm.fit(matchfs[(matchfs < upper)])
p = norm.cdf(matchfs, loc, scale)
In [271]:
keep, padj = fdrcorrection(1-p, alpha=0.1)
In [272]:
keep.sum()
Out[272]:
In [273]:
keep.mean()
Out[273]:
In [274]:
import scipy.stats as stats
#measurements = np.random.normal(loc = 20, scale = 5, size=100)
stats.probplot(matchfs[~keep], dist="norm", plot=plt);
In [275]:
import scipy.stats as stats
#measurements = np.random.normal(loc = 20, scale = 5, size=100)
stats.probplot(matchfs[(matchfs < upper)], dist="norm", plot=plt);
In [164]:
skip.mean()
Out[164]:
In [165]:
plt.hist(p[skip], 100);
In [166]:
plt.hist(p[p> 0.9], 100);
In [138]:
plt.hist(padj, 100);
In [117]:
var
Out[117]:
Explore called instances¶
In [349]:
keep = matchfs > fdr_threshold_norm_right(matchfs)
In [350]:
plt.hist(matchfs[keep], 100);
In [351]:
g = sns.jointplot(matchfs[keep], importancefs[keep, 0], alpha=.1);
g.set_axis_labels("Match", "Importance");
In [351]:
g = sns.jointplot(matchfs[keep], importancefs[keep, 0], alpha=.1);
g.set_axis_labels("Match", "Importance");
In [ ]:
plt.
Instance calling¶
- sufficiently good match in one of the tasks
In [92]:
hyp_contrib_score_scan.shape
Out[92]:
In [115]:
np.all(contrib_score_scan==np.maximum(contrib_score_scan_fwd, contrib_score_scan_rev))
Out[115]:
In [116]:
np.all(hyp_contrib_score_scan==np.maximum(hyp_contrib_score_scan_fwd, hyp_contrib_score_scan_rev))
Out[116]:
In [122]:
plt.plot(np.ravel(pwm_scan_fwd)[:100000], np.ravel(contrib_score_scan_fwd)[:100000], ".")
Out[122]:
In [119]:
plt.plot(np.ravel(pwm_scan_fwd)[:100000], np.ravel(hyp_contrib_score_scan_fwd)[:100000], ".")
Out[119]:
In [123]:
plt.hist(np.ravel(pwm_scan_fwd), 100);
plt.xlabel("pwm scan score");
In [128]:
plt.plot(np.ravel(pwm_scan_fwd)[:100000], np.ravel(contrib_score_scan_fwd_scale)[:100000], ".")
Out[128]:
In [224]:
plt.plot(np.ravel(contrib_score_scan_fwd)[:100000], np.ravel(contrib_score_scan_fwd_scale)[:100000], ".");
plt.xlabel("Match")
plt.ylabel("scale");
In [226]:
plt.hist(np.ravel(contrib_score_scan_fwd), 100);
In [134]:
plt.plot(np.ravel(contrib_score_scan_fwd)[:100000], np.ravel(hyp_contrib_score_scan_fwd)[:100000], ".")
Out[134]:
In [135]:
plt.plot(np.ravel(contrib_score_scan_fwd_scale)[:100000], np.ravel(hyp_contrib_score_scan_fwd_scale)[:100000], ".")
Out[135]:
Plot scores for the example¶
In [90]:
viz_dict = dict(contrib=all_contrib[example_idx, 400:600],
hyp_contrib=all_hyp_contrib[example_idx, 400:600],
pwm_match=seq_match[example_idx, 400:600],
contrib_match=match[example_idx, 400:600, 0],
contrib_scale=importance[example_idx, 400:600, 0])
In [91]:
from basepair.modisco.results import Seqlet
In [121]:
pos_hits = np.argsort(-match[example_idx, 400:600, 0].max(axis=-1))[:2] # top 2 hits
In [106]:
match[example_idx, 400:600, 0].sum()
Out[106]:
In [109]:
importance[example_idx, 400:600, 0].sum()
Out[109]:
In [96]:
pos_hits
Out[96]:
In [137]:
from basepair.utils import halve
def idx2seqlet(pos, pattern):
# TODO handle the reverse-complementation
i, j = halve(len(pattern))
return Seqlet(None, pos-i, pos+j, name=pattern.name)
In [131]:
seqlets = [idx2seqlet(pos, pattern) for pos in pos_hits]
In [129]:
pattern.plot("seq");
In [132]:
plot_tracks(viz_dict, seqlets, rotate_y=0, legend=True);
Benchmarking the functions¶
In [373]:
%timeit score_region_cont_jaccard(qa, ta)
In [375]:
%timeit score_region_cont_jaccard(qa, ta_all_contrib[:100])
In [376]:
%timeit score_region_cont_jaccard(qa, ta_all_contrib[:1000])
In [377]:
%timeit score_region_cont_jaccard(qa, ta_all_contrib[0]) # 1kb
In [378]:
%timeit score_full_regionarr_with_perpos_continjaccard(ta_all_contrib[0], qa) # 1kb
In [379]:
%timeit score_region_cont_jaccard(qa, ta_all_contrib[0][:, :100]) # 100bp
In [380]:
%timeit score_full_regionarr_with_perpos_continjaccard(ta_all_contrib[0][:, :100], qa) # 100bp
In [315]:
a, a_norm = score_full_regionarr_with_perpos_continjaccard(ta_all_contrib[0][:, :1000], qa)
b, b_norm = score_region_cont_jaccard(qa, ta_all_contrib[0][:, :1000]) # 100bp
In [336]:
assert b_norm[0] == a_norm[0]
In [382]:
%timeit np.stack([score_region_cont_jaccard(qa, ta_all_contrib[i]) for i in range(1000)])
In [390]:
scanned = np.stack(Parallel(10)(delayed(score_region_cont_jaccard)(qa, ta_all_contrib[i]) for i in tqdm(range(len(ta_all_contrib)))))
In [391]:
ta_all_contrib.shape
Out[391]:
In [415]:
# TODO - add unit-tests
In [337]:
other_jacc = jaccard_sim_func(ta_all_contrib[0][:, :16][np.newaxis]*qa_L1_norm / a_norm[0], qa[np.newaxis])[0,0]
In [343]:
np.allclose(other_jacc, a[0])
Out[343]:
In [344]:
np.allclose(other_jacc, b[0])
Out[344]:
In [348]:
other_jacc = jaccard_sim_func(ta_all_contrib[0][:, 1:17][np.newaxis]*qa_L1_norm / b_norm[1], qa[np.newaxis])[0,0]
In [349]:
np.allclose(other_jacc, a[1])
Out[349]:
In [350]:
np.allclose(other_jacc, b[1])
Out[350]:
In [316]:
plt.scatter(a_norm, b_norm)
Out[316]:
In [318]:
plt.scatter(np.ravel(a), np.ravel(b));
In [ ]:
%timeit score_region_cont_jaccard(qa, ta)
In [159]:
%timeit score_full_regionarr_with_perpos_continjaccard(ta, qa)
Pattern scanning for the whole group¶
In [187]:
from basepair.modisco.utils import ic_scale
def scan_pattern(pattern_grp, contrib_scores, hypothetical_contribs, one_hot, task_names, trim_frac=0, verbose=True):
"""
final modisco results scan structure
- sequence
- fwd
- rev
- <task>
- <imp score>
- match
- fwd
- rev
- scale
"""
# TODO - add i and j to the output?
pwm = pattern_grp['sequence']['fwd'][:]
trim_ij = trim_pssm_idx(ic_scale(pwm), frac=trim_frac)
pattern_len = len(pwm)
def get_group_match(grp, name):
for n in grp.keys():
if name in n:
return grp[n]
raise ValueError("{name} doesn't match any group keys: {l}".format(name=name, l=list(grp.keys())))
def parallel_score_continousjaccard_restructure(qa, ta, **kwargs):
fwd, fwd_scale = parallel_score_continousjaccard(qa, ta, **kwargs)
rev, rev_scale = parallel_score_continousjaccard(qa[::-1, ::-1], ta, **kwargs)
return dict(match=dict(fwd=fwd, rev=rev),
scale=dict(fwd=fwd_scale, rev=rev_scale),
)
trim_ij_rc = (pattern_len - trim_ij[1], pattern_len - trim_ij[0])
return dict(sequence=dict(fwd=parallel_score_pssm(pwm[trim_ij[0]:trim_ij[1]], one_hot, verbose=verbose),
rev=parallel_score_pssm(pwm[trim_ij[0]:trim_ij[1]][::-1, ::-1], one_hot, verbose=verbose)),
**{task: {"contrib_scores": parallel_score_continousjaccard_restructure(get_group_match(pattern_grp[task], "contrib_scores")['fwd'][trim_ij[0]:trim_ij[1]], contrib_scores),
"hypothetical_contribs": parallel_score_continousjaccard_restructure(get_group_match(pattern_grp[task], "hypothetical_contribs")['fwd'][trim_ij[0]:trim_ij[1]], hypothetical_contribs),}
for task in tqdm(task_names, disable=not verbose)}
)
In [190]:
tasks
Out[190]:
In [191]:
mp1 = scan_pattern(mr.get_pattern_grp("metacluster_0", "pattern_0"), all_contrib, all_hyp_contrib, all_seq, tasks, trim_frac=0.08)
In [199]:
all_seq.size
Out[199]:
In [217]:
idx = pd.Series(np.arange(all_seq.shape[0])).sample(10000).values
In [218]:
plt.plot(np.ravel(mp1['Oct4']['contrib_scores']['match']['fwd'])[idx], np.ravel(mp1['Sox2']['contrib_scores']['match']['fwd'])[idx], ".")
Out[218]:
In [219]:
plt.plot(np.ravel(mp1['Oct4']['contrib_scores']['match']['fwd'])[idx], np.ravel(mp1['Nanog']['contrib_scores']['match']['fwd'])[idx], ".")
Out[219]:
In [220]:
plt.plot(np.ravel(mp1['Oct4']['contrib_scores']['match']['fwd'])[idx], np.ravel(mp1['Klf4']['contrib_scores']['match']['fwd'])[idx], ".")
Out[220]:
In [221]:
mp1['Oct4']['contrib_scores']['match']['fwd'].shape
Out[221]:
Ideas¶
- how much of the blob can be explained by other motifs?
- what are the important motifs
- first try to assign the
- IDR: check Nathan's code - take from kundajelab/idr
- plot rank-rank
- affinity propagation
- run on the predictions as well as real data
- footprint -> DNA binding domain
- crystal structure?
- oct-sox
- crystal structure?
orthogonal validation?
- expression - MPRA
- Nanog SELEX dataset
- re-analyze the data?
- MPRA - HK27ac
- try to predict the enhancer activity
- stratify the peaks into acetylated vs non-acetylated
higher fraction of the
- allele specific binding?
- train in human on ips lines
transposable elements are not conserved in Human
simulations -> profile gives only a sense of what can be going on
- affinity with SELEX
- MPRA
- acetylation
- Human <-> mouse allelic effects
overlap with the histone data
question:
- mixture model also predicting things that are not bound
- repressors
- question: are there any repressors?
- idea:
- add more tasks to the output