Figure 5

Three things to show for the long motifs

  • % Match to the repbase
    • do the same plot as before
  • merge significantly similar motifs into one
    • merge also the seqlets to make the meta-plots?
  • show the importance scores within each motif
In [445]:
# Imports
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
from basepair.imports import *
from basepair.plot.profiles import multiple_plot_stranded_profile
from basepair.exp.paper.config import tf_colors
from basepair.plot.profiles import multiple_plot_stranded_profile
from basepair.exp.paper.config import tf_colors, tasks
from basepair.plot.table import render_mpl_table
from basepair.exp.paper.fig4 import *
from basepair.plot.tracks import tidy_motif_plot
In [498]:
paper_config()
In [446]:
sdir = Path('/oak/stanford/groups/akundaje/avsec/basepair/data/processed/comparison/output/')

model_dir = sdir / 'nexus,peaks,OSNK,0,10,1,FALSE,same,0.5,64,25,0.004,9,FALSE,[1,50],TRUE'
In [409]:
mr_dict = [(task, ModiscoResult(model_dir / f"deeplift/{task}/out/profile/wn/modisco.h5"))
           for task in tasks]
for k,v in mr_dict: v.open()
patterns_dict = [(task, read_pkl(model_dir / f"deeplift/{task}/out/profile/wn/patterns.pkl"))
                  for task in tasks]
footprints = {t: read_pkl(model_dir / f"deeplift/{t}/out/profile/wn/footprints.pkl")
              for t in tasks}
In [502]:
mr = ModiscoResult(model_dir / f"deeplift/all/out/profile/wn/modisco.h5")
mr.open()
In [503]:
patterns = [mr.get_pattern(pname) for pname in mr.patterns()]
In [447]:
patterns = read_pkl(model_dir / f"deeplift/all/out/profile/wn/patterns.pkl")
In [504]:
def get_seqlets_bedtool(p):
    return BedTool(str(model_dir / f"deeplift/all/out/profile/wn/seqlets/{p.name}.bed.gz"))
In [450]:
patterns2 = [p for p in out if p.attrs['n_seqlets'] > 100]
In [470]:
patterns_long_name = [p for p in patterns if p.attrs['features']['n seqlets'] > 100]

Overlapping seqlet regions with RepeatMasker regions

Load the repeat masker file

In [471]:
from basepair.exp.paper.repeat_masker import download_repeat_masker, read_repeat_masker, intersect_repeat_masker
from pybedtools import BedTool
In [453]:
dfrm = read_repeat_masker(download_repeat_masker())

Using downloaded and verified file: /users/avsec/workspace/basepair/data/raw/annotation/mm10/RepeatMasker/mm10.fa.out.gz

Intersect seqlet locations with repeat masker

In [454]:
from basepair.modisco.core import patterns_to_df
In [455]:
bt_te = BedTool.from_dataframe(dfrm)
In [456]:
dfi_raw = pd.concat(Parallel(n_jobs=20)(delayed(intersect_repeat_masker)(p.name, get_seqlets_bedtool(p), bt_te, f=0.9) for p in tqdm(patterns)))

100%|██████████| 151/151 [00:40<00:00,  3.33it/s]
In [457]:
dfi = dfi_raw.copy()
In [458]:
# Restrict only to LTR/
# dfi = dfi[dfi.repeat_family.str.startswith("LTR/")]
In [459]:
# Append some stats
unique_elements = dfi.groupby(['pattern_name']).interval.nunique()
dfiu = dfi[['pattern_name', 'n_pattern']].drop_duplicates().set_index('pattern_name').join(unique_elements)
dfiu['LTR_overlap_frac'] = dfiu.interval / dfiu.n_pattern
dfi = pd.merge(dfi, dfiu.reset_index()[['pattern_name', 'LTR_overlap_frac']], on='pattern_name', how='left')
dfi = dfi.drop_duplicates()
In [472]:
# Append properties
df = patterns_to_df(patterns_long_name, properties=['name', 'short_name', 'seq_info_content'])
df = df.set_index('name').join(dfiu).reset_index()

Scatter-plot: Seq IC ~ LTR_overlap_frac

In [473]:
TE_min_seq_IC = 30
TE_min_LTR_Overlap_frac = 0.7
In [474]:
height = get_figsize(0.25)[0]
fig = sns.jointplot("LTR_overlap_frac", "seq_info_content", marginal_kws=dict(bins=15), s=5, data=df, height=height);
fig.ax_joint.axvline(TE_min_LTR_Overlap_frac, linestyle='--', color='grey', alpha=.4);
fig.ax_joint.axhline(TE_min_seq_IC, linestyle='--', color='grey', alpha=.4);
# fig.savefig(f"{figures}/TE-match,seq-IC.scatter.pdf")
# fig.savefig(f"{figures}/TE-match,seq-IC.scatter.png")
In [475]:
df['group'] = (df['seq_info_content'] > 30).map({False: 'short', True: 'long'})

Density plot as an alternative to scatter plot

In [476]:
plotnine.options.figure_size = get_figsize(0.25)
(ggplot(aes(x='LTR_overlap_frac', color='group'), data=df) +
 geom_density() +
 scale_color_brewer(type='qual', palette=3) + 
 theme_classic(base_size=10, base_family='Arial')
)
Out[476]:
<ggplot: (-9223363279938752845)>
In [465]:
# TE's are in the top right
# pattern_te_names = df[(df.seq_info_content > TE_min_seq_IC) & (df.LTR_overlap_frac > TE_min_LTR_Overlap_frac)].sort_values('LTR_overlap_frac', ascending=False).name.unique()
In [477]:
# TE's are in the top right
pattern_te_names = df[(df.seq_info_content > TE_min_seq_IC)].sort_values('LTR_overlap_frac', ascending=False).name.unique()
In [478]:
len(pattern_te_names)
Out[478]:
57

We can see that motif with high IC are frequently in the region of known TE's

In [479]:
dfi.repeat_name.value_counts()
Out[479]:
RLTR9E          1618
ID_B1           1008
RLTR13D6         890
                ... 
(TGATGATGG)n       1
LTR81B             1
UCON9              1
Name: repeat_name, Length: 885, dtype: int64
In [480]:
dfi.repeat_family.value_counts()
Out[480]:
LTR/ERVK         11401
LTR/ERVL-MaLR     3824
SINE/Alu          2180
                 ...  
DNA/hAT              1
SINE/tRNA-RTE        1
DNA/TcMar            1
Name: repeat_family, Length: 51, dtype: int64

For each pattern, are all seqlets falling into regions of the same LTRs?

In [481]:
dfi_top = dfi.drop_duplicates()
dfi_top['n_repeat_name'] = dfi_top.groupby(['pattern_name', 'repeat_name']).repeat_name.transform('size')
dfi_top['n_repeat_frac'] = dfi_top['n_repeat_name'] / dfi_top['n_pattern']
dfi_top1 = dfi_top.groupby(['pattern_name']).apply(lambda x: x.loc[x.n_repeat_frac.idxmax()])
In [482]:
del dfi_top1['chrom']
del dfi_top1['start']
del dfi_top1['end']
del dfi_top1['interval']
del dfi_top1['pattern_name']
In [483]:
len(patterns)
Out[483]:
151
In [484]:
len(dfi_top1)
Out[484]:
151

Long patterns and their LTR overlap

In [485]:
long_patterns_clustered = cluster_align_patterns([p for p in patterns_long_name if p.seq_info_content >= 30], n_clusters=len(patterns_long_name))

100%|██████████| 57/57 [00:14<00:00,  4.25it/s]
100%|██████████| 57/57 [00:00<00:00, 105.03it/s]
In [486]:
dft = dfi_top1.loc[[p.name for p in long_patterns_clustered]]
In [487]:
main_patterns = dft.groupby("repeat_name").LTR_overlap_frac.idxmax().values
In [493]:
main_patterns
Out[493]:
array(['metacluster_4/pattern_18', 'metacluster_4/pattern_12', 'metacluster_1/pattern_8',
       'metacluster_2/pattern_2', 'metacluster_3/pattern_11', 'metacluster_3/pattern_12',
       'metacluster_3/pattern_20', 'metacluster_3/pattern_14', 'metacluster_4/pattern_8',
       'metacluster_2/pattern_13', 'metacluster_2/pattern_6', 'metacluster_4/pattern_4',
       'metacluster_3/pattern_9', 'metacluster_2/pattern_10', 'metacluster_2/pattern_14',
       'metacluster_4/pattern_34', 'metacluster_4/pattern_11', 'metacluster_3/pattern_3',
       'metacluster_4/pattern_26', 'metacluster_4/pattern_16', 'metacluster_4/pattern_25',
       'metacluster_6/pattern_11', 'metacluster_4/pattern_10', 'metacluster_1/pattern_16',
       'metacluster_4/pattern_17'], dtype=object)
In [506]:
for p in patterns:
    if p.name not in main_patterns:
        continue
    ymax = max([v.max() for v in p.contrib.values()])
    ymin = min([v.min() for v in p.contrib.values()])
    ltr_name = dft.loc[p.name].repeat_name
    ltr_frac = dft.loc[p.name].LTR_overlap_frac
    p.plot(['seq_ic', 'contrib'], 
           rotate_y=0,
           letter_width=0.1, 
           height=0.4,
           ylim=[(0, 2)] + [(ymin, ymax)]*4,
           title=f"{p.name},  LTR: {ltr_name}, any LTR frac: {ltr_frac:.2f}");
    sns.despine(top=True, right=True, left=False, bottom=True)