tornadoplot 绘制富集热图

Original JunJunLab 老俊俊的生信笔记 2022-08-15

收录于合集 #R包分享 93个

点击上方，关注老俊俊！

1引言

分享个小玩意, tornadoplot 用来绘制富集热图,输入 bigwig 文件和 peaks 文件即可。

使用的是 ggplot 绘图系统。

github 地址:

https://github.com/teunbrand/tornadoplot

2安装

# install.packages("devtools")
devtools::install_github("teunbrand/tornadoplot")

3使用

library(BiocFileCache)
#> Loading required package: dbplyr
bfc <- BiocFileCache()

# Declare files
# These bigwigs are derived from PMID: 28212747
sources <- c(
  "sox2" = "http://dc2.cistrome.org/genome_browser/bw/73468_treat.bw",
  "oct4" = "http://dc2.cistrome.org/genome_browser/bw/73466_treat.bw"
)

# 设置样本名
bigwigs <- setNames(bfcrpath(bfc, sources), c("Sox2 ChIP", "Oct4 ChIP"))

peaks 文件处理代码:

library(GenomicRanges)
library(rtracklayer)

# Files -------------------------------------------------------------------

# These files were downloaded from cistrome: http://cistrome.org/db/#/
# They are bed files from data deposited in GSE74112 from PMID: 28212747
# The data is chromatin immunoprecipitation followed by sequencing (ChIP-seq),
# where the proteins Sox2 and Oct4/Pou5f1 have been precipitated from
# mouse embryonic stem cells (mESC).

dir <- file.path("/DATA", "users", "t.vd.brand", "test_data")
files <- file.path(
  dir,
  c("73466_peaks_oct4.bed",
    "73468_peaks_sox2.bed")
)

# Data import -------------------------------------------------------------

si <- SeqinfoForUCSCGenome("mm10")
si <- keepStandardChromosomes(si, "Mus_musculus")

# We don't want any of the scaffold sequences

bed <- lapply(files, import)
bed <- lapply(bed, keepStandardChromosomes, "Mus_musculus", "coarse")
bed <- as(bed, "GRangesList")
bed <- setNames(bed, c("oct4", "sox2"))

# We want to find sites where Sox2 and Oct4 co-bind and where they bind
# uniquely.

red <- reduce(stack(bed))
red$sox2 <- overlapsAny(red, bed$sox2)
red$oct4 <- overlapsAny(red, bed$oct4)
red$cat <- ifelse(red$sox2 & red$oct4, "both",
                  ifelse(red$sox2, "sox2", "oct4"))
red <- GRanges(seqnames(red), IRanges(start(red), end(red)),
               seqinfo = si, cat = red$cat)

# Choosing sites ----------------------------------------------------------

set.seed(0)
octsox_peaks <- GRangesList(
  Sox2 = sample(red[red$cat == "sox2"], 500),
  Oct4 = sample(red[red$cat == "oct4"], 500),
  Both = sample(red[red$cat == "both"], 500)
)
octsox_peaks <- sort(octsox_peaks)

也可以拿网上的示例数据 Rdata,load 一下即可:

https://github.com/teunbrand/tornadoplot/tree/master/data

画图:

library(tornadoplot)

# These are a preprocessed set of ChIP-seq peaks
feats <- octsox_peaks

# Recall that the files are named, which the next function uses for labels
bigwigs <- BigWigFileList(bigwigs)

# Give the function locations and data
tornado_plot(features = feats, data = bigwigs)

看着 g 里 g 气的。

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,数据代码已上传至QQ群,欢迎加入下载。