Plot tracks abve the signal¶
In [ ]:
# TODO - list files to use
In [1]:
import pandas as pd
import numpy as np
from pybedtools import BedTool
from basepair.config import get_data_dir, create_tf_session
from tqdm import tqdm
from concise.preprocessing import encodeDNA
from basepair.datasets import get_sox2_data
from basepair.plots import plot_loss
import pyBigWig
from basepair.math import softmax
import matplotlib.pyplot as plt
from basepair.datasets import *
from basepair.datasets import sox2_oct4_peaks_sox2
from basepair import samplers
from concise.utils.plot import seqlogo_fig, seqlogo
from imp import reload
import h5py
from modisco.visualization import viz_sequence
from matplotlib import pyplot as plt
ddir = get_data_dir()
In [2]:
train, valid, test = sox2_oct4_peaks_sox2()
In [3]:
modisco_file = f"{ddir}/processed/chipnexus/motifs/sox2-oct4/modisco/counts.hdf5"
In [4]:
!ls {ddir}/processed/chipnexus/motifs/*/modisco/*
In [5]:
modisco_file = f"{ddir}/processed/chipnexus/motifs/oct4/modisco/counts.hdf5"
In [6]:
from kipoi.readers import HDF5Reader
In [7]:
f = HDF5Reader(modisco_file)
In [8]:
f.open()
In [9]:
f.f['/metacluster_idx_to_submetacluster_results/metacluster0/seqlets'][:][0]
Out[9]:
In [10]:
f.ls()
Out[10]:
In [11]:
import yaml
In [12]:
f.close()
In [36]:
from basepair.modisco import nan_like, parse_seqlet, ModiscoResult
In [14]:
from concise.utils.pwm import PWM
In [37]:
mr = ModiscoResult(modisco_file)
In [16]:
r = mr.extract_signal(valid[0])
In [17]:
def ic_scale(x):
background = np.array([0.27, 0.23, 0.23, 0.27])
return viz_sequence.ic_scale(x, background=background)
In [18]:
seqlogo_fig(ic_scale(r['pattern_1'].mean(axis=0)), figsize=(10,2));
In [19]:
# TODO - extract the signal
In [20]:
signal = mr.extract_signal(valid[1]['sox2'])['pattern_0'].mean(axis=0)
sequence = ic_scale(mr.extract_signal(valid[0])['pattern_0'].mean(axis=0))
In [21]:
# TODO - assemble figures into one?
In [22]:
def plot_profiles(x, y, mresult):
for pattern in mr.patterns():
signal = mr.extract_signal(y)[pattern].mean(axis=0)
sequence = ic_scale(mr.extract_signal(x)[pattern].mean(axis=0))
fig, (ax0, ax1) = plt.subplots(2, 1, sharex=True, figsize=(20,2))
ax0.set_title(f"{pattern}")
ax0.plot(signal[:,0], label='fwd')
ax0.plot(signal[:,1], label='rev')
ax0.legend()
seqlogo(sequence, ax=ax1)
In [23]:
mr.close()
In [ ]:
from basepair.modisco import ModiscoResult
Sox2 - signal¶
In [2]:
fpath = f"{ddir}/processed/chipnexus/motifs/sox2/modisco/multi-task.0.valid.h5"
mr = ModiscoResult(fpath)
incl = np.load(fpath+".valid.npy")
mr.plot_profiles(valid[0][incl], {"Sox2": valid[1]['sox2'][incl],
"Oct4": valid[1]['oct4'][incl]},
rc_vec=[False, True],
start_vec = [15, 15],
width=40,
legend=False,
ylim=[0, 7.45],
seq_height=1.5,
fpath_template='fig/icml18/motifs/sox2-signal.motif_{}',
figsize=(6,2.5))
mr.close()
Oct4 - signal¶
In [228]:
fpath = f"{ddir}/processed/chipnexus/motifs/oct4/modisco/multi-task.1.valid.h5"
mr = ModiscoResult(fpath)
df_oct4 = pd.DataFrame.from_dict(mr.seqlets()['pattern_0'])
df_oct4['id'] = np.arange(len(df_oct4))
df_oct4['chr'] = df_oct4.example.astype(str)
In [229]:
fpath = f"{ddir}/processed/chipnexus/motifs/sox2/modisco/multi-task.0.valid.h5"
mr = ModiscoResult(fpath)
df_sox2 = pd.DataFrame.from_dict(mr.seqlets()['pattern_1'])
df_sox2['id'] = np.arange(len(df_sox2))
df_sox2['chr'] = df_sox2.example.astype(str)
In [230]:
from pybedtools import BedTool
In [234]:
int_oct4 = BedTool.from_dataframe(df_oct4[['chr', 'start', 'end']])
In [236]:
int_sox2 = BedTool.from_dataframe(df_sox2[['chr', 'start', 'end']])
In [238]:
intersect =int_oct4.intersect(int_sox2, wa=True, u=True)
In [240]:
intersect_b =int_sox2.intersect(int_oct4, wa=True, u=True)
In [239]:
len(intersect)
Out[239]:
In [242]:
len(int_sox2)
Out[242]:
In [243]:
len(int_oct4)
Out[243]:
In [241]:
len(intersect_b)
Out[241]:
In [244]:
# Hence ~80/300 intervals overlap
In [246]:
80/300
Out[246]:
In [218]:
m = pd.merge(df_sox2, df_oct4, on='example', how='outer')
In [227]:
df_sox2['chr'] = df_sox2.example.astype(str)
Out[227]:
In [224]:
new_df = m.loc[(m["start_x"] >= m["start_y"]) & (m["end_x"] <= m["end_y"])]
In [226]:
new_df
Out[226]:
In [222]:
(m.start_y>m.start_x)*()
Out[222]:
In [207]:
df.head()
Out[207]:
Sox2 - signal¶
In [310]:
fpath = f"{ddir}/processed/chipnexus/motifs/sox2/modisco/multi-task.0.valid.h5"
mr = ModiscoResult(fpath)
incl = np.load(fpath+".valid.npy")
mr.plot_profiles(valid[0][incl], {"Sox2": valid[1]['sox2'][incl],
"Oct4": valid[1]['oct4'][incl]},
rc_vec=[False, True],
start_vec = [15, 15],
width=40,
legend=False,
ylim=[[0, 3.5],[0, 7.45]],
fpath_template='fig/icml18/motifs/sox2-signal.motif_{}',
n_bootstrap=100,
seq_height=1.5,
figsize=(6,2.5))
mr.close()
Oct4 - signal¶
In [311]:
fpath = f"{ddir}/processed/chipnexus/motifs/oct4/modisco/multi-task.1.valid.h5"
mr = ModiscoResult(fpath)
incl = np.load(fpath+".valid.npy")
mr.plot_profiles(valid[0][incl], {"Sox2": valid[1]['sox2'][incl],
"Oct4": valid[1]['oct4'][incl]},
rc_vec=[True, True],
start_vec = [18, 10],
width=40,
legend=False,
seq_height=1.5,
ylim=[[0, 3.5],[0, 7.45]],
fpath_template='fig/icml18/motifs/oct4-signal.motif_{}',
n_limit = 35,
n_bootstrap=100,
figsize=(6,2.5))
mr.close()
Sox2 - counts¶
In [312]:
fpath = f"{ddir}/processed/chipnexus/motifs/sox2/modisco/multi-task.2.valid.h5"
mr = ModiscoResult(fpath)
incl = np.load(fpath+".valid.npy")
mr.plot_profiles(valid[0][incl], {"Sox2": valid[1]['sox2'][incl],
"Oct4": valid[1]['oct4'][incl]},
rc_vec=[True, True, True, False],
start_vec = [23, 8, 8, 18],
width=40,
legend=False,
ylim=[[0, 3.5],[0, 7.45]],
fpath_template='fig/icml18/motifs/sox2-counts.motif_{}',
seq_height=1.5,
n_bootstrap=100,
figsize=(6,2.5))
mr.close()
In [302]:
a=1
Oct4 - counts¶
In [313]:
fpath = f"{ddir}/processed/chipnexus/motifs/oct4/modisco/multi-task.3.valid.h5"
mr = ModiscoResult(fpath)
incl = np.load(fpath+".valid.npy")
mr.plot_profiles(valid[0][incl], {"Sox2": valid[1]['sox2'][incl],
"Oct4": valid[1]['oct4'][incl]},
rc_vec=[False, True, False, False],
start_vec=[10, 18,30,20],
width=40,
legend=False,
ylim=[[0, 3.5],[0, 7.45]],
seq_height=1.5,
fpath_template='fig/icml18/motifs/oct4-counts.motif_{}',
n_bootstrap=100,
figsize=(6,2.5))
mr.close()
Test¶
In [134]:
fpath = f"{ddir}/processed/chipnexus/motifs/oct4/modisco/multi-task.0.test.h5"
mr = ModiscoResult(fpath)
incl = np.load(fpath+".test.npy")
mr.plot_profiles(test[0][incl], {"Sox2": test[1]['sox2'][incl],
"Oct4": test[1]['oct4'][incl]},
rc_vec=[True, True, False, False],
#start_vec=[10, 18,20,15],
#width=50,
figsize=(20,3))
mr.close()
In [138]:
fpath = f"{ddir}/processed/chipnexus/motifs/oct4/modisco/multi-task.1.test.h5"
mr = ModiscoResult(fpath)
incl = np.load(fpath+".test.npy")
mr.plot_profiles(test[0][incl], {"Sox2": test[1]['sox2'][incl],
"Oct4": test[1]['oct4'][incl]},
rc_vec=[False, True, False, False],
#start_vec=[10, 18,20,15],
#width=50,
figsize=(20,3))
mr.close()
In [137]:
fpath = f"{ddir}/processed/chipnexus/motifs/oct4/modisco/multi-task.1.valid.h5"
mr = ModiscoResult(fpath)
incl = np.load(fpath+".valid.npy")
mr.plot_profiles(valid[0][incl], {"Sox2": valid[1]['sox2'][incl],
"Oct4": valid[1]['oct4'][incl]},
rc_vec=[True, True, False, False],
#start_vec=[10, 18,20,15],
#width=50,
figsize=(20,3))
mr.close()
Sox2 - counts¶
In [94]:
modisco_file = f"{ddir}/processed/chipnexus/motifs/oct4/modisco/multi-task.3.valid.h5"
In [95]:
hdf5_results = h5py.File(modisco_file)
patterns = (list(hdf5_results
["metacluster_idx_to_submetacluster_results"]
["metacluster0"]
["seqlets_to_patterns_result"]
["patterns"]["all_pattern_names"]))
print(len(patterns))
pattern_grp = (hdf5_results
["metacluster_idx_to_submetacluster_results"]
["metacluster0"]
["seqlets_to_patterns_result"]
["patterns"])
In [100]:
for pattern_name in patterns:
pattern = pattern_grp[pattern_name]
print(pattern_name)
print("total seqlets:",len(pattern["seqlets_and_alnmts"]["seqlets"]))
background = np.array([0.27, 0.23, 0.23, 0.27])
print("fwd:")
viz_sequence.plot_weights(pattern["task0_contrib_scores"]["fwd"])
viz_sequence.plot_weights(pattern["task0_hypothetical_contribs"]["fwd"])
viz_sequence.plot_weights(viz_sequence.ic_scale(np.array(pattern["sequence"]["fwd"]),
background=background))
print("reverse:")
viz_sequence.plot_weights(pattern["task0_contrib_scores"]["rev"])
viz_sequence.plot_weights(pattern["task0_hypothetical_contribs"]["rev"])
viz_sequence.plot_weights(viz_sequence.ic_scale(np.array(pattern["sequence"]["rev"]),
background=background))
