In [1]:
exp = 'nexus,peaks,OSNK,0,10,1,FALSE,same,0.5,64,25,0.004,9,FALSE,[1,50],TRUE'
imp_score = 'profile/wn'
In [2]:
task = "Nanog"
pattern = "metacluster_0/pattern_6"
Goal¶
- get PFM instances genome wide
Tasks¶
- TODO
- [x] export all patterns as MEME / FIMO motifs
- check if you have some code ready
- use the trimmed patterns (same as used for dfi scanning)
- [X] Do the FIMO scanning -> generates the bed files
- do it from the notebook or write a Snakemake file?
- fimo.smk
- do it from the notebook or write a Snakemake file?
- [X] load the fimo files as pd.DataFrame
- use your scripts from motif comparisons
- [X] overlap instances (join by chr start end strand)
- [X] load dfi
- is there a significant overlap?
- which ones overlap?
- which ones don't overlap?
- [X] Give the following columns for each motif (pd.DataFrame):
- pattern
#seqlets- is_te
- IC
#PFM genome-wide#PFM peaks- % PFM peaks in CWM
#CWM peaks- % CWM peaks in PFM
#all peaks
- [ ] export as csv to gdrive. load with gsheet
- [ ] subset the gsheet to only TE's in Figure 3B
- [x] export all patterns as MEME / FIMO motifs
- analysis
- compare with CWM and PWM in peaks
- give spreadsheet and or histogram
In [3]:
# Imports
from basepair.imports import *
from basepair.exp.paper.config import *
import pyranges as pr
hv.extension('bokeh')
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
paper_config()
In [4]:
ddir = get_data_dir()
model_dir = models_dir / exp
Export MEME motifs¶
In [5]:
patterns_dict = [(task, read_pkl(model_dir / f"deeplift/{task}/out/{imp_score}/patterns.pkl"))
for task in tasks]
In [409]:
p_meme = None
for task, patterns in patterns_dict:
mr = ModiscoResult(model_dir / f'deeplift/{task}/out/{imp_score}/modisco.h5')
for pname in mr.patterns():
untrimmed_pattern = mr.get_pattern(pname)
pattern = untrimmed_pattern.trim_seq_ic(0.08)
if task == 'Nanog' and pattern.name == 'metacluster_0/pattern_6':
print(len(pattern))
p_meme = untrimmed_pattern
motif_dir = model_dir / f"deeplift/{task}/out/{imp_score}/PFM/{pattern.name}"
motif_dir.mkdir(parents=True, exist_ok=True)
pattern.write_meme_file([0.25, 0.25, 0.25, 0.25],
motif_dir /"meme.txt")
Run FIMO¶
In [ ]:
# RUN
# snakemake -s fimo.smk -j 10
TODO¶
- [X] test the load time of a single motif vs full table
- fastparquet doesn't work properly :/ -> that's probably due to '/' in the category names
- write a function to load dfi for each pattern
- read_instances('Klf4', 'metacluster_0/pattern_1')
- write a function to perform the overlaps for each pattern separately
- use snakemake to run in parallel
In [5]:
# Load stats about the motif
mr = ModiscoResult(model_dir / f'deeplift/{task}/out/{imp_score}/modisco.h5')
# get TE's
motifs_te = [pattern
for pattern in mr.patterns()
if mr.get_pattern(pattern).seq_info_content > 30 and mr.n_seqlets(*pattern.split('/')) > 100]
is_te = pattern in motifs_te
# motifs_te_d = OrderedDict([(shorten_pattern(x), shorten_pattern(x)) for x in motifs_te])
n_seqlets = mr.n_seqlets(*pattern.split("/"))
p = mr.get_pattern(pattern).trim_seq_ic(0.08)
What fraction of the seqlets contains a CWM hit (and vice versa)?¶
- [X] convert the cwm matches into intervals using seqname strand
- [X] convert seqlets to pyranges
- [ ] perform the overlaps for a single pattern
- What fraction of the seqlets contains a CWM hit?
- [ ] perform the overlaps for multiple patterns
- all patterns
- major patterns
load_instancesshould load the full thing- load the major pattern names
- use
load_instancesto filter the instances only to the major motifs - include also the TE's
- make sure you are only looking at 57 main motifs (e.g. 5 must be excluded...)
In [7]:
# TOOD - there is a discrepancy between the intervals
TODO - use the absolute coordinates¶
In [6]:
from basepair.modisco.pattern_instances import load_instances, multiple_load_instances
def abs_coord(seqlet, ranges):
"""Create absolute seqlet coordinates
"""
row = ranges.loc[seqlet.seqname]
return {"Chromosome": row.chrom,
"Start": row.start + seqlet.start - 15, # make it 70 bp long
"End": row.start + seqlet.end + 15,
}
def dfi2pr(dfi):
dfi['Chromosome'] = dfi.example_chrom
dfi['Start'] = dfi.pattern_start_abs
dfi['End'] = dfi.pattern_end_abs
dfi['Strand'] = dfi.strand
dfi_pr = pr.PyRanges(dfi)
return dfi_pr
df_subset = pd.read_csv("https://docs.google.com/spreadsheets/d/1MVBUbMs_F0UPgZxR8-v9eMM6DCFhVYc7kdH-C6WJnYU/export?gid=0&format=csv")
df_subset.name = df_subset.name.str.replace('"', '')
df_subset.task = df_subset.task.str.replace('"', '')
In [8]:
df_info = pd.read_csv("https://docs.google.com/spreadsheets/d/1uAtJZJR_dmlWwer5Afcd2q72Uuj7A4Ln/export?gid=826411817&format=csv")
In [16]:
Out[16]:
In [110]:
len(dfi_pr_all.intersect(dfs, how='containment', strand=None))
Out[110]:
In [111]:
len(dfi_pr_all)
Out[111]:
In [127]:
len(dfi_pr_all.drop_duplicate_positions(strand=False))
Out[127]:
In [125]:
len(dfi_pr_main.intersect(dfs, how='containment', strand=None).drop_duplicate_positions(strand=False))
Out[125]:
In [112]:
len(dfi_pr_main)
Out[112]:
In [122]:
len(dfs.drop_duplicate_positions(strand=False))
Out[122]:
In [123]:
len(dfs.overlap(dfi_pr_all).drop_duplicate_positions(strand=False))
Out[123]:
In [129]:
len(dfs.overlap(dfi_pr_main).drop_duplicate_positions(strand=False))
Out[129]:
In [ ]:
# TODO - remove motif with low ic
In [17]:
for task in tasks:
print("Loading ranges")
# Load seqlets
mr = ModiscoResult(model_dir / f'deeplift/{task}/out/{imp_score}/modisco.h5')
seqlets = mr.all_seqlets()
ranges = mr.load_ranges()
ranges.set_index('idx', inplace=True)
assert np.all(ranges.index == np.arange(len(ranges)))
dfs = pd.DataFrame([abs_coord(s, ranges) for s in seqlets])
dfs = pr.PyRanges(dfs).drop_duplicate_positions(strand=False)
print("loading dfi")
# Load CWM instances
print("load instances")
dfi = load_instances(model_dir / f'deeplift/{task}/out/{imp_score}/instances.parq')
print("filter")
dfi = dfi.query('match_weighted_p > .2').query('imp_weighted_p > .01')
print("filter")
# That's only subsetting the main motifs for a particular task
motifs = {k:k for k in df_subset[df_subset.task == task].pattern}
all_motifs = {k:k for k in df_info[(df_info.IC > 4) & (df_info.task == task)].pattern.map(shorten_pattern)}
dfi_main_motifs = load_instances(dfi, motifs)
dfi_pr_main = dfi2pr(dfi_main_motifs).drop_duplicate_positions(strand=False)
dfi_all_motifs = load_instances(dfi, all_motifs)
dfi_pr_all = dfi2pr(dfi_all_motifs).drop_duplicate_positions(strand=False)
print("-"*40)
print(f"Task: {task}")
s = len(dfs)
sm = len(dfs.overlap(dfi_pr_main).drop_duplicate_positions(strand=False))
sa = len(dfs.overlap(dfi_pr_all).drop_duplicate_positions(strand=False))
print(f"seqlets overlaps main motifs: {sm}/{s} ({sm/s:.2f})")
print(f"seqlets overlaps all motifs: {sa}/{s} ({sa/s:.2f})")
print("seqlets overlaps main motifs: {}".format(len(dfs.overlap(dfi_pr_main))))
print("")
cm = len(dfi_pr_main)
cmo = len(dfi_pr_main.intersect(dfs, how='containment', strand=None).drop_duplicate_positions(strand=False))
ca = len(dfi_pr_all)
cao = len(dfi_pr_all.intersect(dfs, how='containment', strand=None).drop_duplicate_positions(strand=False))
print(f"CWM (main motifs) overlap seqlets: {cmo}/{cm} ({cmo/cm:.2f})")
print(f"CWM (all motifs) overlap seqlets: {cao}/{ca} ({cao/ca:.2f})")
print("-"*40)
In [25]:
# total seqlets
s = 42995+17665+101304+103514
s
Out[25]:
In [27]:
#overlap main motifs
mo = 28767 + 6377 + 51361 + 59243
ao = 37159 + 13144 + 87802 + 99336
In [28]:
print(f"total seqlets: {mo}")
print(f"main motif overlaps: {mo} / {s} ({mo/s:.2f})")
print(f"all motif overlaps: {ao} / {s} ({ao/s:.2f})")
In [407]:
p.plot();
In [408]:
p_meme.plot();
In [356]:
# Load CWM instances
from basepair.modisco.pattern_instances import load_instances
dfi = load_instances(model_dir / f'deeplift/{task}/out/{imp_score}/instances.parq',
{shorten_pattern(pattern): shorten_pattern(pattern)})
# filter only to the relevant ones
dfi = dfi.query('match_weighted_p > .2').query('imp_weighted_p > .01')
dfi['Chromosome'] = dfi.example_chrom
dfi['Start'] = dfi.pattern_start_abs
dfi['End'] = dfi.pattern_end_abs
dfi['Strand'] = dfi.strand
dfi_pr = pr.PyRanges(dfi)
In [357]:
# Load PWM instances
from basepair.exp.chipnexus.comparison import read_fimo_dfi
dfi_fimo = read_fimo_dfi(model_dir / f'deeplift/{task}/out/{imp_score}/PFM/{pattern}/fimo.tsv')
dfi_fimo['pattern'] = f"{task}/{pattern}"
dfi_fimo['Chromosome'] = dfi_fimo.chrom
dfi_fimo['Start'] = dfi_fimo.start - 1 # intervals are 1-based
dfi_fimo['End'] = dfi_fimo.end
dfi_fimo['Strand'] = dfi_fimo.strand
dfi_fimo = pr.PyRanges(dfi_fimo)
In [358]:
# Load the peak intervals
from basepair.cli.modisco import load_ranges
ranges = load_ranges(model_dir / f'deeplift/{task}/out/{imp_score}')
ranges['Chromosome'] = ranges.chrom
ranges['Start'] = ranges.start
ranges['End'] = ranges.end
ranges = pr.PyRanges(ranges)
In [359]:
ranges
Out[359]:
Perform the overlaps¶
In [360]:
dfi_fimo_peaks = pr.PyRanges(dfi_fimo.intersect(ranges, how='containment').as_df().reset_index())
In [361]:
stats = {}
stats['pattern'] = pattern
stats['task'] = task
stats['imp_score'] = imp_score
stats['is_te'] = is_te
stats['n_seqlets'] = n_seqlets
stats['IC'] = p.seq_info_content
stats['length'] = len(p)
stats['n_pwm_gw'] = len(dfi_fimo)
stats['n_pwm_peaks'] = len(dfi_fimo_peaks)
stats['n_cwm_peaks'] = len(dfi_pr)
stats['n_cwm_and_pwm_peaks'] = len(dfi_pr.intersect(dfi_fimo_peaks, how='containment'))
stats['frac_cwm_in_pwm_peaks'] = stats['n_cwm_and_pwm_peaks'] / stats['n_cwm_peaks']
stats['frac_pwm_in_cwm_peaks'] = stats['n_cwm_and_pwm_peaks'] / stats['n_pwm_peaks']
stats
Out[361]:
In [371]:
len(dfi_pr.overlap(dfi_fimo_peaks))
Out[371]:
In [373]:
dfi_pr.overlap(dfi_fimo_peaks).df
Out[373]:
In [375]:
dfi_pr['chr1', '+', 45358193:45358263].df
Out[375]:
In [376]:
dfi_fimo_peaks['chr1', '+', 45358193:45358263].df
Out[376]:
In [378]:
45358263 - 45358197
Out[378]:
Why do they differ?
In [380]:
len(p)
Out[380]:
In [384]:
!cat {model_dir}/deeplift/{task}/out/{imp_score}/PFM/{pattern}/meme.txt
In [366]:
dfi_pr.sort().df.head()
Out[366]:
In [365]:
dfi_fimo_peaks.sort().df.head()
Out[365]:
In [ ]:
# why is n_cwm_and_pwm_peaks 0?
In [315]:
# write to the json file
In [354]:
from basepair.utils import write_json
# write_json(stats, model_dir / f'deeplift/{task}/out/{imp_score}/PFM/{pattern}/overlaps.json')
Load the results¶
In [331]:
dfr = pd.read_csv(model_dir / 'deeplift/overlap_cwm_pwm.csv')
In [332]:
dfr.head()
Out[332]:
In [336]:
dict(dfr.mean())
Out[336]:
In [334]:
dict(dfr[dfr.is_te].mean())
Out[334]:
In [341]:
!cp {model_dir}/deeplift/overlap_cwm_pwm.csv figures/modisco/{exp}/overlap_cwm_pwm.csv
In [343]:
dfr.to_excel(f'figures/modisco/{exp}/overlap_cwm_pwm.xlsx')
In [410]:
!ls -latrh figures/modisco/{exp}/overlap_cwm_pwm.xlsx