In [1]:
from basepair.imports import *
ddir = get_data_dir()
In [2]:
from basepair.plot.profiles import extract_signal
from basepair.math import softmax
from basepair.plot.heatmaps import heatmap_stranded_profile, multiple_heatmap_stranded_profile
from basepair.plot.profiles import plot_stranded_profile, multiple_plot_stranded_profile
import plotnine
from plotnine import *
import statsmodels.api as sm
import statsmodels.formula.api as smf
In [338]:
from basepair.modisco.results import Seqlet, resize_seqlets
from basepair.modisco.core import dfi2seqlets, annotate_profile
from basepair.cli.modisco import load_profiles
In [337]:
model_dir = Path(f"{ddir}/processed/chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/")
modisco_pdir = model_dir / "modisco/by_peak_tasks/weighted/"
In [34]:
# Load the data
d = HDF5Reader(model_dir / "grad.all.h5")
d.open()
In [6]:
dfi = pd.read_csv(f"{modisco_pdir}/Nanog/instances.tsv.gz", sep='\t')
In [7]:
dfi.tail()
Out[7]:
In [9]:
seqlets_medium = dfi2seqlets(dfi[dfi.match_weighted_cat == 'medium'])
In [10]:
dfi.match_weighted.plot.hist(100);
In [11]:
dfi.match_weighted_p.plot.hist(10);
In [12]:
dfi.match_weighted_cat.value_counts().plot.bar();
In [13]:
dfi.imp_weighted_p.plot.hist(100);
In [14]:
dfi.imp_weighted_cat.value_counts().plot.bar();
In [15]:
dfi.match_weighted_p.plot.hist(10);
In [16]:
from basepair.stats import low_medium_high
Re-compute the motif instances¶
In [57]:
p = mr.get_pattern(pattern).trim_seq_ic(0.08)
In [76]:
dfi.head()
Out[76]:
Nanog¶
In [86]:
task = "Nanog"
pattern = "metacluster_0/pattern_0"
Dev - compute the profile similarity¶
In [324]:
# load profiles
profiles = load_profiles(modisco_pdir, model_dir/'grad.all.h5')
In [ ]:
tasks = list(profiles)
In [345]:
dfi_anno = annotate_profile(dfi, mr, profiles)
In [351]:
dfi_anno['Klf4/profile_match'].isnull().mean()
Out[351]:
In [363]:
dfi_anno['match_weighted_cat'] = pd.Categorical(dfi_anno['match_weighted_cat'])
In [364]:
dfi_anno['imp_weighted_cat'] = pd.Categorical(dfi_anno['imp_weighted_cat'])
In [369]:
len(dfi_anno)
Out[369]:
In [370]:
len(dfi)
Out[370]:
In [392]:
dfif = dfi_anno[~np.isinf(dfi_anno['Oct4/profile_match'])]
In [393]:
dfif.groupby(['match_weighted_cat', 'imp_weighted_cat'])['Oct4/profile_match'].mean().plot.bar()
Out[393]:
In [371]:
ggplot(aes(x='match_weighted_cat',
color='imp_weighted_cat',
y='Oct4/profile_match'), dfi_anno.sample(100000)) + geom_boxplot()
In [367]:
ggplot(aes(x='match_weighted_cat',
color='imp_weighted_cat',
y='Oct4/profile_match'), dfi_anno.sample(100000)) + geom_boxplot()
Out[367]:
In [354]:
dfi_anno.head()
Out[354]:
In [344]:
np.isinf(np.array([np.inf]))
Out[344]:
In [285]:
mr = ModiscoResult(modisco_pdir / f"{task}/modisco.h5")
mr.open()
seqlets = mr._get_seqlets(pattern, trim_frac=0.08)
seqlets = resize_seqlets(seqlets, 200, seqlen=profile_obs.shape[1])
# mr.close()
In [325]:
list(profiles)
Out[325]:
In [286]:
len(seqlets)
Out[286]:
In [180]:
len(seqlets)
Out[180]:
In [181]:
from basepair.cli.modisco import load_ranges, load_included_samples
In [182]:
ranges = load_ranges(modisco_pdir / task)
In [183]:
load_included_samples(modisco_pdir / task)
Out[183]:
In [184]:
include_samples = np.load(read_json(modisco_pdir / f"{task}/kwargs.json")["filter_npy"])
In [187]:
out_task = 'Oct4'
profile_obs = d.f[f'/targets/profile/{out_task}'][:][include_samples]
In [188]:
ds = DataSpec.load(model_dir / "dataspec.yaml")
In [273]:
seqlets = dfi2seqlets(dfi[(dfi.match_weighted_cat == 'medium')])
seqlets = resize_seqlets(seqlets, 200, seqlen=profile_obs.shape[1])
In [274]:
seqlets[:3]
Out[274]:
In [291]:
seqlet_profile_obs = extract_signal(profile_obs, seqlets)
In [292]:
total_counts = seqlet_profile_obs.sum(axis=-1).sum(axis=-1)
sort_idx = np.argsort(-total_counts)
In [309]:
avg_profile = seqlet_profile_obs.mean(axis=0)
In [312]:
# now, compute the KL distance
In [310]:
plot_stranded_profile(avg_profile)
In [311]:
plot_stranded_profile(p.mean(axis=0))
In [300]:
# probabilities
p = seqlet_profile_obs[sort_idx] / seqlet_profile_obs[sort_idx].sum(axis=1, keepdims=True)
# drop NA's
notnan = ~np.any(np.any(np.isnan(p), axis=-1), axis=-1)
total_counts = seqlet_profile_obs[sort_idx].sum(axis=1)[notnan]
p = p[notnan]
seqlet_idx = np.array([s.seqname for s in seqlets])[notnan]
In [302]:
total_counts.shape
Out[302]:
In [279]:
# dropped
print("Dropped", seqlet_profile_obs.shape[0] - p.shape[0], "profiles with only 0's")
In [280]:
p.shape
Out[280]:
In [281]:
old_total_counts = total_counts
In [303]:
plt.scatter(total_counts[:,0], total_counts[:,1])
Out[303]:
In [282]:
plt.plot(total_counts.sum(axis=-1)[:400])
plt.plot(old_total_counts.sum(axis=-1)[:400], label='old')
plt.legend()
Out[282]:
In [283]:
heatmap_stranded_profile(p[:5000], figsize=(20,20))
In [272]:
heatmap_stranded_profile(p[:5000], figsize=(20,20))
1. Quantify the profile effect -> entropy and total counts¶
In [201]:
from scipy.stats import entropy
from scipy.special import rel_entr, kl_div
In [202]:
# S(p_obs)
entropies = entropy(p.swapaxes(0,1)).sum(axis=-1)
# KL(p_obs, p_average)
kl = kl_div(p, p.mean(axis=0, keepdims=True)).mean(axis=-1).sum(axis=-1)
crossentropy = rel_entr(p, p.mean(axis=0, keepdims=True)).mean(axis=-1).sum(axis=-1)
In [204]:
fig = plt.figure(figsize=(13,5))
plt.subplot(131)
plt.plot(entropies);
plt.ylabel("entropy")
plt.xlabel("idx");
plt.subplot(132)
plt.plot(kl, entropies, "."); # kl divergence and the entropy between the other factor are almost the same
plt.xlabel("kl")
plt.ylabel("Entropy");
plt.subplot(133)
plt.plot(crossentropy, kl, "."); # kl divergence and the entropy between the other factor are almost the same
plt.xlabel("crossentropy")
plt.ylabel("Entropy");
In [205]:
fig=plt.figure(figsize=(10,4))
plt.subplot(121)
plt.plot(entropies, np.log(1+total_counts.mean(axis=-1)), ".", alpha=0.3);
plt.xlabel("Entropy")
plt.ylabel("log(1+ counts)");
plt.subplot(122)
plt.plot(entropies**2, np.log(1+total_counts.mean(axis=-1)), ".", alpha=0.3);
plt.xlabel(r"Entropy^2")
plt.ylabel("log(1+ counts)");
Conclusion¶
- both metrics: counts and entropy are good for characterizing the signal
In [206]:
df = mr.seqlet_df_instances()
In [23]:
dfp = df.pivot_table("center", "seqname", "pattern", aggfunc=len, fill_value=0)
In [24]:
count_features = dfp.loc[seqlet_idx]
In [26]:
count_features[pattern].value_counts().plot.bar();
plt.xlabel("Number of occurences in the sequence")
plt.ylabel("Frequency");
In [28]:
count_features.iloc[:,count_features.columns!=pattern].sum().plot.bar(figsize=(20,5));
plt.ylabel("Sum");
Boxplot for each factor the change in entropy¶
In [29]:
def rename_pattern(p):
if "metacluster" in p:
return p.replace("metacluster_", "m").replace("/", "_").replace("pattern_", "p")
else:
return p
In [30]:
count_features.columns = [rename_pattern(p) for p in count_features.columns]
In [31]:
dfm = count_features.assign(entropy=entropies, counts=np.log10(1 + total_counts.mean(axis=-1)),
example_idx=count_features.index).melt(id_vars=['entropy', 'counts', 'example_idx'], var_name="pattern")
In [32]:
dfmf = dfm.groupby("pattern").filter(lambda x: x.value.sum()> 10)
In [33]:
dfmf.value = pd.Categorical(dfmf.value)
dfm.value = pd.Categorical(dfm.value)
In [45]:
plotnine.options.figure_size = (20,10)
ggplot(aes(x='value', y='entropy'), dfmf) + geom_boxplot() + facet_wrap("~pattern", ncol = 10, scales='free_x') + theme_bw()
Out[45]:
Fit a model to determine the effects¶
In [35]:
def ols_formula(df, dependent_var, *excluded_cols):
'''
Generates the R style formula for statsmodels (patsy) given
the dataframe, dependent variable and optional excluded columns
as strings
'''
df_columns = list(df.columns.values)
df_columns.remove(dependent_var)
for col in excluded_cols:
df_columns.remove(col)
return dependent_var + ' ~ ' + ' + '.join(df_columns)
In [36]:
dm = count_features.assign(counts=np.log10(1+total_counts.sum(axis=-1)))
In [37]:
ols_formula(dm, "counts")
Out[37]:
In [38]:
results = smf.ols(ols_formula(dm, "counts"), data=dm).fit()
In [39]:
def tidy_ols(ols_reults):
smry = results.summary()
coef = smry.tables[1]
return pd.DataFrame(coef.data[1:], columns=coef.data[0])
In [40]:
df_fit = tidy_ols(results)
In [41]:
# significant coefficients
df_fit_signif = df_fit[df_fit['P>|t|'].astype(float) < 0.05]
df_fit_signif = df_fit_signif[df_fit_signif[""] != "Intercept"] # don't show the intercept
df_fit_signif.iloc[df_fit_signif['coef'].astype(float).abs().argsort()].iloc[::-1] # sort by the effect size
Out[41]:
In [42]:
pattern_short = rename_pattern(pattern)
single_motif_idx = count_features.index[count_features[pattern_short] == 1]
df_center = df[df.seqname.isin(single_motif_idx)].query(f"pattern == '{pattern}'")[['seqname', 'center']]
df_counts = pd.DataFrame({"seqname": seqlet_idx,
"log_counts": np.log10(1+total_counts.mean(axis=-1))})
dfd = pd.merge(df[df.pattern != pattern], df_center, on='seqname', suffixes=("_other", "_core"))
dfd['rel'] = dfd.center_other - dfd.center_core
dfd = dfd.merge(df_counts, on="seqname")
In [50]:
plotnine.options.figure_size = (20,20)
ggplot(aes(x="rel", y='log_counts'), dfd) + geom_point(alpha=0.5) + facet_wrap("~pattern", ncol = 5) + theme_bw() #+ xlim([-400, 400])
Out[50]:
In [49]:
plotnine.options.figure_size = (20,20)
ggplot(aes(x="rel", y='log_counts'), dfd) + geom_point(alpha=0.5) + facet_wrap("~pattern", ncol = 5) + theme_bw() + xlim([-70, 70])
Out[49]:
Conclusions¶
- not much going on