Goal¶
- get the motif spacing information
Questions¶
- [x] How often do motifs co-occur in the sequences?
- [ ] What is their distance distribution?
- [ ] Are there signs of nearby weak binding sites in the presence of one major motif?
- [ ] Does modisco actually flag those weak binding sites?
TODO¶
- [x] make venn-like plot using upsetr
- [ ] get the motif spacing information
- [ ] for a particular combination of motifs, plot the distance distribution
TODO¶
- add number of sequences
Conclusions¶
-
Motif spacing¶
In [65]:
import pandas as pd
import numpy as np
from basepair.modisco import ModiscoResult
from basepair.config import get_data_dir
In [66]:
from basepair.datasets import sox2_oct4_peaks_sox2
In [67]:
train, valid, test = sox2_oct4_peaks_sox2()
In [68]:
ddir = get_data_dir()
In [84]:
incl.sum()
Out[84]:
How often do motifs co-occur in the sequences?¶
In [108]:
def plot_modisco_results(pattern, split):
fpath = pattern.format(split=split)
mr = ModiscoResult(fpath)
incl = np.load(fpath + f".{split}.npy")
print(f"# sequences: {incl.sum()}")
mr.stats()
mr.plot_profiles(valid[0][incl], {"Sox2": eval(split)[1]['sox2'][incl],
"Oct4": eval(split)[1]['oct4'][incl]},
#rc_vec=[False, True],
#start_vec = [15, 15],
#width=40,
legend=False,
#ylim=[0, 7.45],
#seq_height=1.5,
n_bootstrap=100,
fpath_template=None,
figsize=(6,2.5))
mr.plot_seqlet_upset()
mr.close()
In [98]:
mdir = f"{ddir}/processed/chipnexus/motifs"
Sox2 - signal¶
In [99]:
plot_modisco_results(f"{mdir}/sox2/modisco/multi-task.0.{{split}}.h5",
"valid")
Oct4 - signal¶
In [100]:
plot_modisco_results(f"{mdir}/oct4/modisco/multi-task.1.{{split}}.h5",
"valid")
Sox2 - counts¶
In [109]:
plot_modisco_results(f"{mdir}/sox2/modisco/multi-task.2.{{split}}.h5",
"valid")
In [106]:
plot_modisco_results(f"{mdir}/sox2/modisco/multi-task.2.{{split}}.h5",
"valid")
Oct4 - counts¶
In [107]:
plot_modisco_results(f"{mdir}/oct4/modisco/multi-task.3.{{split}}.h5",
"valid")
In [42]:
mr.f.ls()
Out[42]:
In [79]:
mr.f.f['/multitask_seqlet_creation_results/task_name_to_coord_producer_results/task0/vals_to_threshold'][:4]
Out[79]:
In [83]:
valid[0].shape
Out[83]:
In [50]:
d.shape
Out[50]:
In [38]:
sl = mr.seqlets()
In [15]:
## TODO - add no patterns to the plot
In [136]:
bdir = f"{ddir}/processed/chipnexus/motifs/sox2-oct4/modisco/bed"
In [137]:
mkdir -p {bdir}
In [138]:
mr.export_seqlets_bed(f"{bdir}/seqlet", True)
In [141]:
!head {bdir}/seqlet.pattern_0.bed
In [25]:
import pandas as pd
In [68]:
instances = []
for p, seqlets in sl.items():
for seqlet in seqlets:
seqlet["id"] = p
instances.append(seqlet)
In [85]:
dfi = pd.DataFrame(instances)
dfi['n'] = dfi.groupby('example').example.transform('size')
dfi.groupby('example').size().plot.hist()
Out[85]:
In [80]:
examples = {}
for p, seqlets in sl.items():
for seqlet in seqlets:
seqlet["id"] = p
if seqlet['example'] in examples:
examples[seqlet['example']].append(seqlet)
else:
examples[seqlet['example']] = [seqlet]
In [83]:
examples2 = {e: l for e,l in examples.items() if len(l) == 2}
In [89]:
examples20 = {e: sorted([x['id'] for x in l]) for e,l in examples.items() if len(l) == 2}
In [92]:
c = {}
for v in examples20.values():
if tuple(v) not in c:
c[tuple(v)] = 1
else:
c[tuple(v)] +=1
In [93]:
c
Out[93]:
In [96]:
ids = [k for k,v in examples20.items() if v == ['pattern_0', 'pattern_1']]
In [99]:
len(ids)
Out[99]:
In [109]:
examples2_01 = [(l[0]['end'] + l[0]['start'])/2 - (l[1]['end'] + l[1]['start'])/2 for e,l in examples2.items() if e in ids]
In [112]:
import matplotlib.pyplot as plt
In [124]:
plt.hist(examples2_01, bins=30);
In [122]:
plt.hist((dfi.end + dfi.start)/2, 30);
plt.xlabel("Seqlet position distribution")
Out[122]:
In [86]:
len(examples2)
Out[86]:
In [87]:
examples2
Out[87]:
In [70]:
dfi.head()
Out[70]:
In [ ]:
pd.pivot()
In [ ]:
df[dfi.n==2].pivot("example", "id", "n")
In [79]:
dfi[dfi.n==2]
Out[79]:
In [76]:
dfi[dfi.n==2].groupby('example')['id'].value_counts()
Out[76]:
In [49]:
n_examples.index[n_examples==2]
Out[49]:
In [52]:
dfi.set_index("example", inplace=True)
In [54]:
n_examples[n_examples==2]
Out[54]:
In [51]:
dfi.loc[n_examples.index[n_examples==2]]
Out[51]: