Goal¶
- export bigwigs for Melanie and generate the logo-plots
In [1]:
import basepair
from basepair.config import get_data_dir, create_tf_session
from keras.models import load_model
from basepair.datasets import *
from basepair.datasets import sox2_oct4_peaks_sox2
from basepair import datasets
from basepair.preproc import AppendTotalCounts, transform_data, resize_interval
from basepair.plots import regression_eval
from basepair.BPNet import BPNetPredictor
from basepair.utils import write_pkl
Setup BPNetPredictor¶
Get the model¶
In [2]:
create_tf_session(0)
Out[2]:
In [3]:
model = load_model("/users/avsec/workspace/basepair/basepair/../data/processed/chipnexus/exp/models/p-multi-task/seq_mutlitask_filters=32,conv1_kernel_size=21,tconv_kernel_size=25,n_dil_layers=6,use_profile=True,use_counts=True,c_task_weight=10,lr=0.004.5.h5")
Get the preprocessor¶
In [4]:
data = sox2_oct4_peaks_sox2()
In [5]:
y = data[0][1]
# rename y -> 'profile/Sox2', 'profile/Oct4
ynew = {"profile/Sox2": y["sox2"], "profile/Oct4": y["oct4"]}
In [6]:
data[2][2].iloc[1788]
Out[6]:
In [7]:
preproc = AppendTotalCounts()
preproc.fit(ynew)
In [8]:
# Write to a file
write_pkl(preproc,
"/users/avsec/workspace/basepair/basepair/../data/processed/chipnexus/exp/models/p-multi-task/preproc.pkl")
Visualize the regions¶
In [9]:
bp = BPNetPredictor(model,
"/mnt/data/pipeline_genome_data/mm10/mm10_no_alt_analysis_set_ENCODE.fasta",
tasks = ["Sox2", "Oct4"],
preproc=preproc)
In [21]:
# All required paths
#!cat /users/avsec/workspace/basepair-workflow/models/0/dataspec.yaml
bws = {task: {s: f"/users/avsec/workspace/basepair-workflow/data/{task}/5prime.counts.{s}.bw" for s in ['pos', 'neg']}
for task in ["Sox2", "Oct4"]}
genome_file = "/mnt/data/pipeline_genome_data/mm10/mm10.chrom.sizes"
klf4_bed = "/users/avsec/workspace/basepair-workflow/data/klf4_sites_mm10.bed"
oct4_bed = "/users/avsec/workspace/basepair-workflow/data/oct4_sites_mm10.bed"
klf4_oct4_windowed_bed = "/users/avsec/workspace/basepair-workflow/data/klf4_oct4_windowed_sites_mm10.bed"
In [22]:
cat {klf4_bed} {oct4_bed} | bedtools makewindows -w 201 -b stdin -i srcwinnum > {klf4_oct4_windowed_bed}
In [23]:
!head {klf4_oct4_windowed_bed}
In [4]:
cat /users/avsec/workspace/basepair-workflow/data/klf4_sites_mm10.bed
In [24]:
!ls /users/avsec/workspace/basepair-workflow/data/
In [25]:
bt = BedTool(klf4_bed).cat(BedTool(oct4_bed), postmerge=False)
intervals = list(bt)
bt_windowed = BedTool(klf4_oct4_windowed_bed)
intervals_windowed = [resize_interval(interval, 201) for interval in bt_windowed]
In [26]:
intervals
Out[26]:
In [27]:
# Interval widths
[interval.stop - interval.start for interval in intervals]
Out[27]:
In [28]:
intervals_resized = [resize_interval(interval, 201) for interval in intervals]
In [20]:
# To be still managable, I resized interval around the center to have width=201
bp.predict_plot([Interval("chr17", 35503946, 35503946+201),
Interval("chr17", 35503946+150, 35503946+201+150),
], bws=bws,
profile_grad='weighted',
figsize=(20, 10))
In [20]:
# To be still managable, I resized interval around the center to have width=201
bp.predict_plot([Interval("chr17", 35503946, 35503946+201),
Interval("chr17", 35503946+150, 35503946+201+150),
], bws=bws,
profile_grad='weighted',
figsize=(20, 10))
In [18]:
# To be still managable, I resized interval around the center to have width=201
bp.predict_plot([Interval("chr17", 35503946, 35503946+201),
Interval("chr17", 35503946+150, 35503946+201+150),
], bws=bws, figsize=(20, 10))
In [59]:
# To be still managable, I resized interval around the center to have width=201
bp.predict_plot(intervals_resized, bws=bws, figsize=(20, 10))
Export Klf4 and Oct4 sites to bigwig¶
- this will generate the following files in
output_dir:<Task>.preds.pos.bw- prediction of positive counts<Task>.preds.neg.bw- prediction of negative counts<Task>.importance.profile.bw- importance scores for the max profile<Task>.importance.counts.bw- importance scores for the total counts
- files: http://mitra.stanford.edu/kundaje/avsec/chipnexus/BPNet/for_melanie/klf4_oct4_sites/
- washu session: https://epigenomegateway.wustl.edu/browser/?genome=mm10&session=RkHGEeJL2F&statusId=874323549
- 2-for-melanie-KLF4-enhancer
- 3-for-melanie-Oct4-Enhancers
In [30]:
[interval.stop - interval.start for interval in bt_windowed][:5]
Out[30]:
In [47]:
output_dir=f"{ddir}/www/BPNet/for_melanie/klf4_oct4_sites/"
!mkdir -p {output_dir}
bp.export_bw(intervals=list(BedTool(intervals_windowed).sort()),
batch_size=64,
output_dir=output_dir)
In [186]:
# Save the bedtool as well
bt.saveas(f"{output_dir}/intervals.bed")
Out[186]:
In [191]:
# Convert to the hammock file format
!basepair bed2hammock {output_dir}/intervals.bed
In [55]:
# Intervals
[f"{i.chrom}:{i.start}-{i.stop} {i.name}" for i in BedTool(f"{output_dir}/intervals.bed")]
Out[55]:
In [3]:
!ls {ddir}/
Export the whole dataset to the memory¶
In [34]:
oct4_bed = "/users/avsec/workspace/basepair-workflow/data/Oct4/peaks.bed.gz"
In [35]:
sox2_bed = f"{BED_DIR}//Sox2_123b_1_ppr.IDR0.05.filt.summit_centered_200bp.narrowPeak"
In [36]:
!wc -l {sox2_bed}
In [37]:
!head {sox2_bed}
In [38]:
btm = BedTool(oct4_bed).cat(sox2_bed, postmerge=False).sort()
In [39]:
btm.fn
Out[39]:
In [40]:
btm_windowed = "/tmp/pybedtools.w_g5cabn_n.tmp"
In [41]:
!bedtools makewindows -w 201 -b {btm.fn} -i srcwinnum > {btm_windowed}
In [42]:
intervals = [resize_interval(i, 201) for i in BedTool(btm_windowed).sort()]
In [43]:
len(intervals)
Out[43]:
In [44]:
pd.Series([i.stop-i.start for i in intervals]).unique()
Out[44]:
In [47]:
intervals = list(BedTool(intervals + intervals_windowed).sort())
In [62]:
output_dir=f"{ddir}/www/BPNet/for_melanie/all_sox2_oct4_sites_weighted_scaled/"
!mkdir -p {output_dir}
bp.export_bw(intervals=intervals,
batch_size=32,
profile_grad='weighted',
scale_importance=True,
output_dir=output_dir)
In [ ]:
from basepair.cli.washu import write_washu_config
In [53]:
BedTool(intervals).saveas(f"{output_dir}/intervals.bed")
Out[53]:
In [48]:
!basepair bed2hammock {output_dir}/intervals.bed
In [55]:
ls {output_dir}/*.bw
In [56]:
tasks = ['Oct4', 'Sox2']
features = ['preds.pos', 'preds.neg', 'importance.counts', 'importance.profile']
In [57]:
filespec = [
{"name": feat,
"path": f"{task}.{feat}.bw",
"type": "bigwig"}
for task in tasks for feat in features
]
In [58]:
filespec += [{"name": "peaks", "path": "peaks.hammock.gz", "type": "hammock"}]
In [59]:
filespec
Out[59]:
In [60]:
!wc -l {output_dir}/intervals.bed
In [63]:
write_washu_config(filespec,
output_dir,
"http://mitra.stanford.edu/kundaje/avsec/chipnexus/BPNet/for_melanie/all_sox2_oct4_sites_weighted_scaled/",
"mm10")
Screenshots¶
KLF4-Downstream_Enhancer_E1_mm10¶
KLF4-Downstream_Enhancer_E2_mm10¶
Distal and Proximal Oct4 Enhancers, Combined¶
