In [156]:
from basepair.imports import *
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
from basepair.modisco.results import Seqlet, resize_seqlets
from basepair.modisco.pattern_instances import dfi2seqlets, annotate_profile
from basepair.cli.modisco import load_profiles
from basepair.BPNet import BPNetPredictor
from basepair.plot.tracks import plot_tracks, filter_tracks
from basepair.preproc import rc_seq
from copy import deepcopy
from basepair.exp.chipnexus.simulate import (insert_motif, generate_sim, plot_sim, generate_seq,
model2tasks, motif_coords, interactive_tracks, plot_motif_table,
plot_sim_motif_col)
from scipy.fftpack import fft, ifft
import warnings
warnings.filterwarnings("ignore")
hv.extension('bokeh', 'matplotlib')
# interval columns in dfi
interval_cols = ['example_chrom', 'pattern_start_abs', 'pattern_end_abs']
TODO¶
- [x] how were the quantiles computed? on restricted sets
- [x] annotate the in-vivo profiles using the same metrics as in vivo
- [X] use the average profile from modisco as the reference
- [X] sanity check - plot the motif
- [x] add also simulated data to the plot
- [ ] plot the frequency domain (spacing)
- [ ] mask TE's from the analysis
Load the required data¶
In [3]:
model_dir = Path(f"{ddir}/processed/chipnexus/exp/models/oct-sox-nanog-klf/models/n_dil_layers=9/")
modisco_dir = model_dir / "modisco/all/profile/"
mr = ModiscoResult(modisco_dir / "modisco.h5")
mr.open()
# Load the data
d = HDF5Reader(model_dir / "grad.all.h5")
d.open()
In [4]:
# specify motifs to use in the analysis
motifs = OrderedDict([
("Oct4-Sox2", "m0_p0"),
#("Oct4-Sox2-deg", "m6_p8"),
#("Oct4", "m0_p18"),
("Sox2", "m0_p1"),
# ("Essrb", "m0_p2"),
("Nanog", "m0_p3"),
#("Nanog-periodic", "m0_p9"),
("Klf4", "m2_p0"),
])
side_motifs = OrderedDict([
("Oct4-Sox2", ("m0_p0", "TTTGCATAACAA")),
("Sox2", ("m0_p1", "CCATTGTT")),
# ("Essrb", ("m0_p2", "TCAAGGTCA")),
("Nanog", ("m0_p3", "TTGATGGC")),
("Klf4", ("m2_p0", "GGGTGTGG")),
])
In [5]:
# plot them
for tf, p in motifs.items():
mr.get_pattern(longer_pattern(p)).trim_seq_ic(0.08).plot("seq");
plt.title(tf)
In [6]:
from basepair.modisco.pattern_instances import load_instances, filter_nonoverlapping_intervals, plot_coocurence_matrix
In [64]:
dfi_full = pd.read_parquet(f"{modisco_dir}/instances.parq", engine='fastparquet')
In [162]:
# TODO get transposable element locations
dfc_te = pd.read_csv("http://mitra.stanford.edu/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/motif_clustering/te_contrib.csv")
te_patterns = dfc_te[dfc_te['ic pwm mean'] > 0.8].pattern.map(longer_pattern).unique()
In [178]:
dfc_te.sort_values("ic pwm mean")
Out[178]:
In [163]:
dfi_te = load_instances(dfi_full[dfi_full.pattern.isin(te_patterns)], None)
dfi_te = dfi_te[(dfi_te.match_weighted_p > 0.1) & (dfi_te.seq_match > 20)]
In [7]:
dfi = load_instances(dfi_full, motifs)
dfi = filter_nonoverlapping_intervals(dfi)
total_examples = len(dfi.example_idx.unique())
total_examples
Out[7]:
In [8]:
# annotate with profiles
profiles = load_profiles(modisco_dir, model_dir/'grad.all.h5')
dfi_anno = annotate_profile(dfi, mr, profiles, profile_width=70, trim_frac=0.08)
In [9]:
dfi = dfi_anno
Visualize the distribution of different profile properties based on the different categories¶
Categories:
- imp_weighted_cat
- match_weighted_cat
In [10]:
dfi_subset = dfi.query('pattern_name=="Oct4-Sox2"')
In [11]:
profile_features = [f for f in dfi_subset if f.startswith("Oct4/profile") and f.endswith("_p")]
categories = ['imp_weighted_cat', 'match_weighted_cat','seq_match_cat']
dfm = dfi_subset.melt(id_vars=categories + ['id'], value_vars=profile_features)
dfm = dfm.dropna()
dfm = dfm[~dfm.variable.str.contains("match")]
dfm.value = dfm.value.astype(float)
In [12]:
dfi_subset[dfi_subset.match_weighted_p > 0.1].plot.scatter("imp_weighted_p", 'Oct4/profile_counts_p', alpha=0.05)
plt.title("Counts vs importance for match_weighted_p > 0.1");
In [13]:
ggplot(aes(x='match_weighted_cat', y='value', color='imp_weighted_cat'), dfm) + geom_boxplot() + \
facet_grid("variable~.", scales='free_y') + theme_classic() + ggtitle("High importance ~ higher profile counts")
Out[13]:
Plot the profiles - sanity check¶
In [14]:
dfi_subset = dfi[(dfi.pattern_center > 400) & (dfi.pattern_center < 600)].query('match_weighted_p > .5').\
query('imp_weighted_p > 0').query('pattern_name=="Oct4-Sox2"')
In [15]:
seqlets = dfi2seqlets(dfi_subset)
seqlets = resize_seqlets(seqlets, 70, seqlen=1000)
seqlet_profiles = {k: extract_signal(v, seqlets) for k,v in profiles.items()}
In [16]:
multiple_plot_stranded_profile({p:v for p,v in seqlet_profiles.items()}, figsize_tmpl=(2.55,2))
multiple_heatmap_stranded_profile(seqlet_profiles, sort_idx=np.arange(1000), figsize=(10,10));
In [17]:
# pattern lengths
dfi.groupby(['pattern_name', 'pattern_len']).size()
Out[17]:
In vivo plots¶
In [18]:
from basepair.modisco.pattern_instances import construct_motif_pairs
In [19]:
comp_strand_compbination = {
"++": "--",
"--": "++",
"-+": "-+",
"+-": "+-"
}
strand_combinations = ["++", "--", "+-", "-+"]
In [20]:
pairs = []
for i in range(len(motifs)):
for j in range(i, len(motifs)):
pairs.append([ list(motifs)[i], list(motifs)[j], ])
In [24]:
profile_mapping = {
"Oct4-Sox2": "Oct4",
"Sox2": "Sox2",
"Nanog": "Nanog",
"Klf4": "Klf4",
"Essrb": "Oct4"
}
In [25]:
from basepair.stats import smooth_window_agg, smooth_lowess, smooth_gam
In [26]:
# plt.style.use('default')
In [27]:
in_silico = read_pkl(f"{ddir}/processed/chipnexus/simulation/spacing.pkl")
In [28]:
sim_df_d, sim_res_dict_d = in_silico
In [29]:
sim_df_d.keys()
Out[29]:
In [30]:
def swap_orientation(o):
if o=='+':
return "-"
if o=="-":
return "+"
raise ValueError("")
# Very much a hacked function to map simulation to in vivo data
def get_xy_sim_single(sim_df_d, motif_pair, feature, orientation):
# For Nanog, always explicilty swap the orientation
orientation_pair = [orientation[0], orientation[1]]
# HACK! Nanog orientation didn't properly match the orientation
if motif_pair[0] == "Nanog":
orientation_pair[0] = swap_orientation(orientation_pair[0])
if motif_pair[1] == "Nanog":
orientation_pair[1] = swap_orientation(orientation_pair[1])
mp = list(deepcopy(motif_pair))
if orientation_pair[0] == "-":
mp[0] = mp[0] + "/rc"
if orientation_pair[1] == "-":
mp[1] = mp[1] + "/rc"
df = sim_df_d[mp[0]] # choose the central motif
df = df[df.motif == mp[1]] # choose the side motif
df = df[df.distance < 150]
# select the task
df = df[df.task == profile_mapping[motif_pair[0]]]
return df.distance.values, df[feature].values
def get_xy_sim(sim_df_d, motif_pair, feature, orientation):
x1,y1 = get_xy_sim_single(sim_df_d, motif_pair, feature, orientation)
x2,y2 = get_xy_sim_single(sim_df_d, list(reversed(motif_pair)), feature, comp_strand_compbination[orientation])
assert np.all(x1 == x2)
return x1, y1, y2
In [34]:
sim_df_d['Oct4-Sox2'].head()
Out[34]:
In [170]:
# motif_pair = ['Nanog', 'Nanog']
d_dfi_filtered=[]
for motif_pair in tqdm(pairs):
dfi_filtered = (dfi[(dfi.pattern_center > 400) & (dfi.pattern_center < 600)]
.query('match_weighted_p > 0.2')
.query('imp_weighted_p > 0'))
dftw_filt = construct_motif_pairs(dfi_filtered, motif_pair,
features=['match_weighted_p','imp_weighted_p', 'imp_weighted'] + \
[f for f in dfi if "profile" in f])
# assure the right strand combination
dftw_filt[dftw_filt.center_diff < 0]['strand_combination'] = dftw_filt[dftw_filt.center_diff < 0]['strand_combination'].map(comp_strand_compbination)
dftw_filt.center_diff = np.abs(dftw_filt.center_diff)
if motif_pair[0] == motif_pair[1]:
motif_pair = [f"{m}{i}" for i,m in enumerate(motif_pair)]
d_dfi_filtered.append((motif_pair, dftw_filt))
Profile¶
In [171]:
fig_profile, all_axes = plt.subplots(3*len(strand_combinations), len(d_dfi_filtered), figsize=(25,15), sharex=True, sharey='row')
for j, (motif_pair, dftw_filt) in enumerate(d_dfi_filtered):
axes = all_axes[:,j]
if motif_pair[0][:-1] == motif_pair[1][:-1]:
axes[0].set_title("{f}<>{f}".format(f=motif_pair[0][:-1]), fontsize=7)
motif_pair_c = [mp[:-1] for mp in motif_pair]
else:
motif_pair_c = motif_pair
axes[0].set_title("<>".join(motif_pair), fontsize=7)
dftw_filt = dftw_filt[(dftw_filt.center_diff < 150) & (dftw_filt.imp_weighted_p.max(1) > 0.3)]
ymax = max([np.log10(1+dftw_filt[profile_mapping[mp] + "/profile_counts"]).max().max() for mp in motif_pair_c])
ymin = min([np.log10(1+dftw_filt[profile_mapping[mp] + "/profile_counts"]).min().min() for mp in motif_pair_c])
#ymax = dftw_filt.imp_weighted.max().max()
#ymin = dftw_filt.imp_weighted.min().min()
for i, sc in enumerate(strand_combinations):
y1 = np.log10(1+ dftw_filt[dftw_filt.strand_combination==sc][profile_mapping[motif_pair_c[0]] + "/profile_counts"][motif_pair[0]])
y2 = np.log10(1+ dftw_filt[dftw_filt.strand_combination==sc][profile_mapping[motif_pair_c[1]] + "/profile_counts"][motif_pair[1]])
# y1 = dftw_filt[dftw_filt.strand_combination==sc]['imp_weighted'][motif_pair[0]]
# y2 = dftw_filt[dftw_filt.strand_combination==sc]['imp_weighted'][motif_pair[1]]
x = dftw_filt[dftw_filt.strand_combination==sc]['center_diff']
#dm,ym,confi = average_distance(x,y, window=5)
dm1,ym1,confi1 = smooth_lowess(x,y1, frac=0.15)
dm2,ym2,confi2 = smooth_lowess(x,y2, frac=0.15)
#dm,ym, confi = smooth_gam(x,y, 140, 20)
ax = axes[3*i]
ax.hist(dftw_filt[dftw_filt.strand_combination==sc]['center_diff'], np.arange(10, 150, 1));
if j == 0:
ax.set_ylabel(sc)
# second plot
ax = axes[3*i+1]
ax.scatter(x,y1, alpha=0.05, s=8)
if confi1 is not None:
ax.fill_between(dm1, confi1[:,0], confi1[:,1], alpha=0.2)
ax.plot(dm1, ym1, linewidth=2, alpha=0.8)
ax.scatter(x,y2, alpha=0.05, s=8)
if confi2 is not None:
ax.fill_between(dm2, confi2[:,0], confi2[:,1], alpha=0.2)
ax.plot(dm2, ym2, linewidth=2, alpha=0.8)
if j == 0:
ax.set_ylabel(sc)
# third plot, simulated
ax = axes[3*i+2]
sim_x, sim_y1, sim_y2 = get_xy_sim(sim_df_d, motif_pair_c, 'profile/counts_frac', sc)
ax.axhline(1, linestyle="--", color='grey', alpha=0.2)
ax.plot(sim_x, sim_y1, linewidth=1, alpha=0.8)
ax.plot(sim_x, sim_y2, linewidth=1, alpha=0.8)
ax.xaxis.set_minor_locator(plt.MultipleLocator(10))
if j == 0:
ax.set_ylabel(sc)
fig_profile.subplots_adjust(wspace=0, hspace=0)
Importance¶
In [172]:
fig_imp, all_axes = plt.subplots(3*len(strand_combinations), len(d_dfi_filtered), figsize=(25,15), sharex=True, sharey='row')
for j, (motif_pair, dftw_filt) in enumerate(d_dfi_filtered):
axes = all_axes[:,j]
if motif_pair[0][:-1] == motif_pair[1][:-1]:
axes[0].set_title("{f}<>{f}".format(f=motif_pair[0][:-1]), fontsize=7)
motif_pair_c = [mp[:-1] for mp in motif_pair]
else:
motif_pair_c = motif_pair
axes[0].set_title("<>".join(motif_pair), fontsize=7)
dftw_filt = dftw_filt[(dftw_filt.center_diff < 150) & (dftw_filt.imp_weighted_p.max(1) > 0.3)]
#ymax = max([np.log10(1+dftw_filt[profile_mapping[mp] + "/profile_counts"]).max().max() for mp in motif_pair_c])
#ymin = min([np.log10(1+dftw_filt[profile_mapping[mp] + "/profile_counts"]).min().min() for mp in motif_pair_c])
ymax = dftw_filt.imp_weighted.max().max()
ymin = dftw_filt.imp_weighted.min().min()
for i, sc in enumerate(strand_combinations):
#y1 = np.log10(1+ dftw_filt[dftw_filt.strand_combination==sc][profile_mapping[motif_pair_c[0]] + "/profile_counts"][motif_pair[0]])
#y2 = np.log10(1+ dftw_filt[dftw_filt.strand_combination==sc][profile_mapping[motif_pair_c[1]] + "/profile_counts"][motif_pair[1]])
y1 = dftw_filt[dftw_filt.strand_combination==sc]['imp_weighted'][motif_pair[0]]
y2 = dftw_filt[dftw_filt.strand_combination==sc]['imp_weighted'][motif_pair[1]]
x = dftw_filt[dftw_filt.strand_combination==sc]['center_diff']
#dm,ym,confi = average_distance(x,y, window=5)
dm1,ym1,confi1 = smooth_lowess(x,y1, frac=0.15)
dm2,ym2,confi2 = smooth_lowess(x,y2, frac=0.15)
#dm,ym, confi = smooth_gam(x,y, 140, 20)
ax = axes[3*i]
ax.hist(dftw_filt[dftw_filt.strand_combination==sc]['center_diff'], np.arange(10, 150, 1));
if j == 0:
ax.set_ylabel(sc)
# second plot
ax = axes[3*i+1]
ax.scatter(x,y1, alpha=0.05, s=8)
if confi1 is not None:
ax.fill_between(dm1, confi1[:,0], confi1[:,1], alpha=0.2)
ax.plot(dm1, ym1, linewidth=2, alpha=0.8)
ax.scatter(x,y2, alpha=0.05, s=8)
if confi2 is not None:
ax.fill_between(dm2, confi2[:,0], confi2[:,1], alpha=0.2)
ax.plot(dm2, ym2, linewidth=2, alpha=0.8)
if j == 0:
ax.set_ylabel(sc)
# third plot, simulated
ax = axes[3*i+2]
sim_x, sim_y1, sim_y2 = get_xy_sim(sim_df_d, motif_pair_c, 'imp/weighted_frac', sc)
ax.axhline(1, linestyle="--", color='grey', alpha=0.2)
ax.plot(sim_x, sim_y1, linewidth=1, alpha=0.8)
ax.plot(sim_x, sim_y2, linewidth=1, alpha=0.8)
ax.xaxis.set_minor_locator(plt.MultipleLocator(10))
if j == 0:
ax.set_ylabel(sc)
fig_imp.subplots_adjust(wspace=0, hspace=0)
In [51]:
fig_profile.savefig("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/spacing/profile_counts.pdf", dpi=400, bbox_inches='tight')
fig_profile.savefig("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/spacing/profile_counts.png", dpi=300, bbox_inches='tight')
fig_imp.savefig("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/spacing/importance_weighted.pdf", dpi=400, bbox_inches='tight')
fig_imp.savefig("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/spacing/importance_weighted.png", dpi=300, bbox_inches='tight')
Exclude TE's¶
In [173]:
from basepair.preproc import dfint_no_intersection
In [174]:
dfi_filtered = (dfi[(dfi.pattern_center > 400) & (dfi.pattern_center < 600)]
.query('match_weighted_p > 0.2')
.query('imp_weighted_p > 0'))
keep_nonte = dfint_no_intersection(dfi_filtered[interval_cols], dfi_te[interval_cols])
print("not overlapping TE", keep_nonte.mean()) # almost all were kept
dfi_filtered = dfi_filtered[keep_nonte]
d_dfi_filtered=[]
for motif_pair in tqdm(pairs):
dftw_filt = construct_motif_pairs(dfi_filtered, motif_pair,
features=['match_weighted_p','imp_weighted_p', 'imp_weighted'] + \
[f for f in dfi if "profile" in f])
# assure the right strand combination
dftw_filt[dftw_filt.center_diff < 0]['strand_combination'] = dftw_filt[dftw_filt.center_diff < 0]['strand_combination'].map(comp_strand_compbination)
dftw_filt.center_diff = np.abs(dftw_filt.center_diff)
if motif_pair[0] == motif_pair[1]:
motif_pair = [f"{m}{i}" for i,m in enumerate(motif_pair)]
d_dfi_filtered.append((motif_pair, dftw_filt))
Profile¶
In [175]:
fig_profile, all_axes = plt.subplots(3*len(strand_combinations), len(d_dfi_filtered), figsize=(25,15), sharex=True, sharey='row')
for j, (motif_pair, dftw_filt) in enumerate(d_dfi_filtered):
axes = all_axes[:,j]
if motif_pair[0][:-1] == motif_pair[1][:-1]:
axes[0].set_title("{f}<>{f}".format(f=motif_pair[0][:-1]), fontsize=7)
motif_pair_c = [mp[:-1] for mp in motif_pair]
else:
motif_pair_c = motif_pair
axes[0].set_title("<>".join(motif_pair), fontsize=7)
dftw_filt = dftw_filt[(dftw_filt.center_diff < 150) & (dftw_filt.imp_weighted_p.max(1) > 0.3)]
ymax = max([np.log10(1+dftw_filt[profile_mapping[mp] + "/profile_counts"]).max().max() for mp in motif_pair_c])
ymin = min([np.log10(1+dftw_filt[profile_mapping[mp] + "/profile_counts"]).min().min() for mp in motif_pair_c])
#ymax = dftw_filt.imp_weighted.max().max()
#ymin = dftw_filt.imp_weighted.min().min()
for i, sc in enumerate(strand_combinations):
y1 = np.log10(1+ dftw_filt[dftw_filt.strand_combination==sc][profile_mapping[motif_pair_c[0]] + "/profile_counts"][motif_pair[0]])
y2 = np.log10(1+ dftw_filt[dftw_filt.strand_combination==sc][profile_mapping[motif_pair_c[1]] + "/profile_counts"][motif_pair[1]])
# y1 = dftw_filt[dftw_filt.strand_combination==sc]['imp_weighted'][motif_pair[0]]
# y2 = dftw_filt[dftw_filt.strand_combination==sc]['imp_weighted'][motif_pair[1]]
x = dftw_filt[dftw_filt.strand_combination==sc]['center_diff']
#dm,ym,confi = average_distance(x,y, window=5)
dm1,ym1,confi1 = smooth_lowess(x,y1, frac=0.15)
dm2,ym2,confi2 = smooth_lowess(x,y2, frac=0.15)
#dm,ym, confi = smooth_gam(x,y, 140, 20)
ax = axes[3*i]
ax.hist(dftw_filt[dftw_filt.strand_combination==sc]['center_diff'], np.arange(10, 150, 1));
if j == 0:
ax.set_ylabel(sc)
# second plot
ax = axes[3*i+1]
ax.scatter(x,y1, alpha=0.05, s=8)
if confi1 is not None:
ax.fill_between(dm1, confi1[:,0], confi1[:,1], alpha=0.2)
ax.plot(dm1, ym1, linewidth=2, alpha=0.8)
ax.scatter(x,y2, alpha=0.05, s=8)
if confi2 is not None:
ax.fill_between(dm2, confi2[:,0], confi2[:,1], alpha=0.2)
ax.plot(dm2, ym2, linewidth=2, alpha=0.8)
if j == 0:
ax.set_ylabel(sc)
# third plot, simulated
ax = axes[3*i+2]
sim_x, sim_y1, sim_y2 = get_xy_sim(sim_df_d, motif_pair_c, 'profile/counts_frac', sc)
ax.axhline(1, linestyle="--", color='grey', alpha=0.2)
ax.plot(sim_x, sim_y1, linewidth=1, alpha=0.8)
ax.plot(sim_x, sim_y2, linewidth=1, alpha=0.8)
ax.xaxis.set_minor_locator(plt.MultipleLocator(10))
if j == 0:
ax.set_ylabel(sc)
fig_profile.subplots_adjust(wspace=0, hspace=0)
Importance¶
In [176]:
fig_imp, all_axes = plt.subplots(3*len(strand_combinations), len(d_dfi_filtered), figsize=(25,15), sharex=True, sharey='row')
for j, (motif_pair, dftw_filt) in enumerate(d_dfi_filtered):
axes = all_axes[:,j]
if motif_pair[0][:-1] == motif_pair[1][:-1]:
axes[0].set_title("{f}<>{f}".format(f=motif_pair[0][:-1]), fontsize=7)
motif_pair_c = [mp[:-1] for mp in motif_pair]
else:
motif_pair_c = motif_pair
axes[0].set_title("<>".join(motif_pair), fontsize=7)
dftw_filt = dftw_filt[(dftw_filt.center_diff < 150) & (dftw_filt.imp_weighted_p.max(1) > 0.3)]
#ymax = max([np.log10(1+dftw_filt[profile_mapping[mp] + "/profile_counts"]).max().max() for mp in motif_pair_c])
#ymin = min([np.log10(1+dftw_filt[profile_mapping[mp] + "/profile_counts"]).min().min() for mp in motif_pair_c])
ymax = dftw_filt.imp_weighted.max().max()
ymin = dftw_filt.imp_weighted.min().min()
for i, sc in enumerate(strand_combinations):
#y1 = np.log10(1+ dftw_filt[dftw_filt.strand_combination==sc][profile_mapping[motif_pair_c[0]] + "/profile_counts"][motif_pair[0]])
#y2 = np.log10(1+ dftw_filt[dftw_filt.strand_combination==sc][profile_mapping[motif_pair_c[1]] + "/profile_counts"][motif_pair[1]])
y1 = dftw_filt[dftw_filt.strand_combination==sc]['imp_weighted'][motif_pair[0]]
y2 = dftw_filt[dftw_filt.strand_combination==sc]['imp_weighted'][motif_pair[1]]
x = dftw_filt[dftw_filt.strand_combination==sc]['center_diff']
#dm,ym,confi = average_distance(x,y, window=5)
dm1,ym1,confi1 = smooth_lowess(x,y1, frac=0.15)
dm2,ym2,confi2 = smooth_lowess(x,y2, frac=0.15)
#dm,ym, confi = smooth_gam(x,y, 140, 20)
ax = axes[3*i]
ax.hist(dftw_filt[dftw_filt.strand_combination==sc]['center_diff'], np.arange(10, 150, 1));
if j == 0:
ax.set_ylabel(sc)
# second plot
ax = axes[3*i+1]
ax.scatter(x,y1, alpha=0.05, s=8)
if confi1 is not None:
ax.fill_between(dm1, confi1[:,0], confi1[:,1], alpha=0.2)
ax.plot(dm1, ym1, linewidth=2, alpha=0.8)
ax.scatter(x,y2, alpha=0.05, s=8)
if confi2 is not None:
ax.fill_between(dm2, confi2[:,0], confi2[:,1], alpha=0.2)
ax.plot(dm2, ym2, linewidth=2, alpha=0.8)
if j == 0:
ax.set_ylabel(sc)
# third plot, simulated
ax = axes[3*i+2]
sim_x, sim_y1, sim_y2 = get_xy_sim(sim_df_d, motif_pair_c, 'imp/weighted_frac', sc)
ax.axhline(1, linestyle="--", color='grey', alpha=0.2)
ax.plot(sim_x, sim_y1, linewidth=1, alpha=0.8)
ax.plot(sim_x, sim_y2, linewidth=1, alpha=0.8)
ax.xaxis.set_minor_locator(plt.MultipleLocator(10))
if j == 0:
ax.set_ylabel(sc)
fig_imp.subplots_adjust(wspace=0, hspace=0)
In [177]:
fig_profile.savefig("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/spacing/non-TE.profile_counts.pdf", dpi=400, bbox_inches='tight')
fig_profile.savefig("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/spacing/non-TE.profile_counts.png", dpi=300, bbox_inches='tight')
fig_imp.savefig("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/spacing/non-TE.importance_weighted.pdf", dpi=400, bbox_inches='tight')
fig_imp.savefig("/srv/www/kundaje/avsec/chipnexus/oct-sox-nanog-klf/models/n_dil_layers=9/modisco/all/profile/spacing/non-TE.importance_weighted.png", dpi=300, bbox_inches='tight')