Goal¶
- compute the cutoff within each modisco cluster
Conclusions¶
- The 'high' cutoff was at the lowest 5% percentile of the distribution
- need to modify the code
TODO¶
- [x] append the percentile scores to the table
- [x] determine the cuttofs: low, medium, high based on [0-33, 33-66, 66-100] ranges
- (discard seq_match < 0 for the motifs)
- [ ] add profile scoring to the funcitons
- figure out which metric to use?
- [ ] add profile match as well as total counts (also add the quantiles comparing both to the seqlet distribution)
- [x] for each modisco run, pre-compute these thresholds (make a CLI, embed into the original code)
In [247]:
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
from basepair.modisco.table import ModiscoData
from plotnine import *
import plotnine
In [256]:
model_dir = Path(f"{ddir}/processed//chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/")
In [257]:
modisco_dir = model_dir / "modisco/by_peak_tasks/weighted/Oct4"
In [251]:
mr = ModiscoResult(modisco_dir / "modisco.h5")
mr.open()
In [7]:
# load all of modisco data
md = ModiscoData.load(modisco_dir, model_dir / "grad.all.h5")
In [254]:
# load all of modisco data
mr = ModiscoResult(model_dir / "modisco/valid/modisco.h5")
mr.open()
In [258]:
model_dir = Path(f"{ddir}/processed//chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/")
In [259]:
modisco_dir = model_dir / "modisco/by_peak_tasks/weighted/Nanog"
In [260]:
dfm = md.get_centroid_seqlet_matches
In [ ]:
dfm =
In [262]:
dfm.head()
Out[262]:
In [ ]:
# TODO - get the whole table
In [265]:
dfm = md.get_centroid_seqlet_matches()
In [269]:
dfm.to_csv(modisco_dir / "centroid_seqlet_matches.csv.bak")
In [161]:
pattern = 'metacluster_0/pattern_0'
In [162]:
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"),
]
Load instances data¶
In [163]:
cached = True
In [164]:
cache_dir = Path(Path(ddir) / "cache/dev/")
cache_dir.mkdir(parents=True, exist_ok=True)
cache_file = cache_dir / "hit-scoring-contrib-only.hdf5"
In [165]:
dfp = pd.read_hdf(cache_file, "/scores")
In [166]:
# add some additional features to dfp
dfp['idx'] = np.arange(len(dfp))
dfp = dfp[dfp.score_seq_match > 0]
# -----------------------------------------------------------------
dfp['log_imp_weighted'] = np.log(1+dfp['imp_weighted'])
dfp['log_match_weighted'] = np.log(1+dfp['match_weighted'])
# 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()
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'})
dfp_subset = dfp.sample(10000) # focus on the subset of the data
Steps¶
- trim the motif
get the matrix of:
- sequences
- profiles
- contrib scores
scan these sequences using your scores
- Plot the histogram of scores
- Compare these to the obtained scores
In [167]:
pattern_name = 'metacluster_0/pattern_0'
1. Get an trim the motif¶
In [168]:
pattern = md.mr.get_pattern(pattern_name).trim_seq_ic(0.08)
2. Get the matrix of¶
- sequences
- profiles
- contrib scores
In [169]:
task = 'Oct4'
In [170]:
i, j = md.get_trim_idx(pattern_name)
In [171]:
seq = md.get_seq(pattern_name)[:, i:j]
profile = {task: md.get_profile_wide(pattern_name, task) for task in tasks}
contrib = {task: md.get_imp(pattern_name, task, 'profile')[:, i:j] for task in tasks}
In [174]:
plot_tracks(dict(recomputed=contrib[task].mean(axis=0), orig=pattern.contrib[task]));
In [175]:
from basepair.plot.heatmaps import heatmap_stranded_profile, multiple_heatmap_stranded_profile
from basepair.plot.profiles import plot_stranded_profile
In [177]:
plot_stranded_profile(profile[task].mean(axis=0))
In [178]:
heatmap_stranded_profile(profile[task][np.argsort(-profile[task].sum((1,2)))]);
3. Scan for patterns¶
In [179]:
tasks = md.tasks
In [182]:
match, importance = pattern.scan_importance(contrib, hyp_contrib=None, tasks=tasks,
n_jobs=1, verbose=True, pad_mode=None)
seq_match = pattern.scan_seq(seq, n_jobs=1, verbose=True, pad_mode=None)
In [183]:
dfm = pattern.get_instances(tasks, match, importance, seq_match, fdr=1, verbose=True)
In [184]:
print(f"Discarding {np.mean(dfm.score_seq_match < 0):.2%} of points with poor sequence match")
In [185]:
dfm = dfm[dfm.score_seq_match > 0]
In [188]:
dfm.head()
Out[188]:
4. Plot the histogram¶
In [186]:
fig ,axes = plt.subplots(1, 3, figsize=(12, 3))
dfm.match_weighted.plot.hist(100, ax=axes[0]);
axes[0].set_xlabel("Match");
dfm.imp_weighted.plot.hist(100, ax=axes[1]);
axes[1].set_xlabel("Importance");
dfm.score_seq_match.plot.hist(100, ax=axes[2]);
axes[2].set_xlabel("Sequence match");
plt.tight_layout()
5. Compare the previous¶
In [187]:
fig ,axes = plt.subplots(1, 3, figsize=(12, 3))
dfmp = dfp_subset[(dfp_subset.pattern == pattern.name) & (dfp_subset.match_weighted > 0.4)]
dfmp = dfp[(dfp.pattern == pattern.name) & (dfp.match_weighted > 0.2)]
dfmp.match_weighted.plot.hist(100, ax=axes[0]);
axes[0].set_xlabel("Match");
dfmp.imp_weighted.plot.hist(100, ax=axes[1]);
axes[1].set_xlabel("Importance");
dfmp.score_seq_match.plot.hist(100, ax=axes[2]);
axes[2].set_xlabel("Sequence match");
plt.tight_layout()
5. Wrap into a function. Run for all patterns of interest¶
In [192]:
from tqdm import tqdm
In [195]:
def get_all_dfm_seqlets(md, trim_frac=0.08, n_jobs=1):
return pd.concat([get_dfm_seqlets(md, pattern, trim_frac, n_jobs)
for pattern in tqdm(md.mr.patterns())])
def get_dfm_seqlets(md, pattern_name, trim_frac=0.08, n_jobs=1, verbose=False):
"""
Args:
md: ModiscoData
pattern_name
"""
tasks = md.tasks
pattern = md.mr.get_pattern(pattern_name).trim_seq_ic(trim_frac)
i, j = md.get_trim_idx(pattern_name)
seq = md.get_seq(pattern_name)[:, i:j]
profile = {task: md.get_profile_wide(pattern_name, task) for task in tasks}
contrib = {task: md.get_imp(pattern_name, task, 'profile')[:, i:j] for task in tasks}
match, importance = pattern.scan_importance(contrib, hyp_contrib=None, tasks=tasks,
n_jobs=n_jobs, verbose=False, pad_mode=None)
seq_match = pattern.scan_seq(seq, n_jobs=n_jobs, verbose=False, pad_mode=None)
dfm = pattern.get_instances(tasks, match, importance, seq_match, fdr=1, verbose=verbose, plot=verbose)
dfm = dfm[dfm.score_seq_match > 0]
return dfm
In [206]:
d = pd.Series(np.arange(10))
In [227]:
quantiles = [0, .33, .66, 1.0]
labels = ['low', 'medium', 'high']
In [234]:
percentile_steps = np.arange(0, 1.01, 0.01)
In [236]:
labels = percentile_steps[1:]
In [241]:
dfm = pd.read_csv(f"{modisco_dir}/../Nanog/centroid_seqlet_matches.csv")
In [237]:
q = d.quantile(quantiles)
In [238]:
pd.cut(d, q.values, labels=labels, include_lowest=True, right=True)
Out[238]:
In [ ]:
# TODO - get the quantile compared to another distribution - p_seqlet
In [ ]:
d.buck
In [201]:
dfm = get_all_dfm_seqlets(md)
In [205]:
md.mr.fpath
Out[205]:
In [203]:
dfm.tail()
Out[203]:
In [200]:
def plot(dfm_seqlets, dfm_orig, title):
dfm = dfm_seqlets
fig ,axes = plt.subplots(1, 3, figsize=(12, 3))
dfm.match_weighted.plot.hist(100, ax=axes[0]);
axes[0].set_xlabel("Match");
axes[0].set_title(title);
dfm.imp_weighted.plot.hist(100, ax=axes[1]);
axes[1].set_xlabel("Importance");
dfm.score_seq_match.plot.hist(100, ax=axes[2]);
axes[2].set_xlabel("Sequence match");
plt.tight_layout()
fig ,axes = plt.subplots(1, 3, figsize=(12, 3))
dfmp = dfm_orig
if len(dfmp) == 0:
import pdb
pdb.set_trace()
threshold = dfmp[dfmp['match_cat']=='high'].match_weighted.min()
dfmp.match_weighted.plot.hist(100, ax=axes[0]);
axes[0].set_xlabel("Match");
axes[0].axvline(threshold, color='red')
dfmp.imp_weighted.plot.hist(100, ax=axes[1]);
threshold = dfmp[dfmp['imp_cat']=='high'].match_weighted.min()
axes[1].set_xlabel("Importance");
axes[1].axvline(threshold, color='red')
dfmp.score_seq_match.plot.hist(100, ax=axes[2]);
axes[2].set_xlabel("Sequence match");
plt.tight_layout()
In [153]:
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 [154]:
for motif, pattern_name in pattern_names:
print(motif)
plot(get_dfm_seqlets(md, pattern_name, verbose=True), dfp_subset[(dfp_subset.pattern == pattern_name)], motif)
In [155]:
thresholds = {"Klf4": 0.55,
"Nanog": 0.45,
"Sox2": 0.55,
"Errb": 0.4,
"Oct4-Sox2": 0.4
}