表观调控13张图之四，peaks区域注释分类比例
表观调控13张图之三转录组非标准分析之MA图，logFC散点图，韦恩图
我把表观调控数据分析，拆分成为了13张图，分别录制为13个视频，即将免费发布在B站，这个期间我们的视频编辑师还在兢兢业业的奋斗，希望这13张图能带领大家学会表观调控数据分析的一般流程， 然后应用到自己的课题哈!同时我也招募了4位视频审查员和细节知识补充员。接下来我们就连载第一位视频审查员的13个笔记：
第一位视频审查员大家也许并不陌生了，早在2018-08-29我发布给学徒的ATAC-seq数据实战（附上收费视频） 他就学完了全部课程内容，还写了笔记在简书，大家搜索二货潜，就可以看到。
本文目录：
批量对Peak进行注释
01：环境准备
02：Peak注释
注释Peak函数编写
批量进行注释
本章节我们的视频审查员（刘博-二货潜）将继续带领大家学习视频，并且复现附件中Figure S1b的一张图，如下：
这里我们是下载文章附件中的 Peak 文件。
"我们要学会变通，不要只会 copy 代码，要学会举一反三。不是背代码，而是去探寻规律。"因为文章当时做的时候是采用d3版的注释信息，而我们是下载了文章附件提供的Peak文件，所以我们这里分析要用d3版本的注释信息。
文章 ChIP-seq 分析方法：ChIP-Seq data were first processed using Trim Galore, then aligned to reference genome (dm3, BDGP Release 5) using Bowtie v1.1.2, 可以看到基因版本为d3.
由于国内镜像网速普遍不好，所以我们在进行相关R包安装之前需要指定镜像源。
# 运行
file.edit("~/.Rprofile")
# 然后再打开的命令中设置清华源信息
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror = "https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
# 然后点击 Rstudio 右上角的 source 即可
# 放在前面的话：一般要安装什么包直接搜索包名，对应的包的 manual 说怎么装就怎么装
BiocManager::install("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
BiocManager::install("org.Dm.eg.db")
BiocManager::install("ChIPseeker")
BiocManager::install("clusterProfiler")
# 加载包
require(ChIPseeker)
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
require(clusterProfiler)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
# 读取当前目录 oldBedFiles 文件夹下的 Peak 文件
bedPeaksFile = './oldBedFiles/Cg_WT.narrowPeak.bed';
bedPeaksFile
# [1] "./oldBedFiles/Cg_WT.narrowPeak.bed"
peak <- readPeakFile( bedPeaksFile)
# 查看 Peak 文件中染色体信息
seqlevels(peak)
# [1] "chr2L" "chr2LHet" "chr2R" "chr2RHet" "chr3L" "chr3LHet" "chr3R" "chr3RHet" "chr4" "chrM" "chrX" "chrXHet" "chrYHet"
# 过滤掉带有 Het 字眼的染色体
# 请留意 `grepl()` 函数
keepChr = !grepl('Het',seqlevels(peak))
seqlevels(peak, pruning.mode = "coarse") <- seqlevels(peak)[keepChr]
pruning.mode 参数有几种选项：c("error", "coarse", "fine", "tidy")
# 默认是 "error"
# 大致说一下每个参数对应的意思（实际上我就是翻译 manual 哈）
# "error"：没有修剪（pruning），并且引发错误
# "coarse": 删除 x 中正在使用的要删除的 seqlevel 的元素
# "fine": 支持 `list` 格式，删除要删除的序列上的范围，不影响原来文件的长度和排序，修剪后的对象与原始对象一样。
# "tidy":
# emmm, 解释的过为牵强，还是有需要的看原文吧。
# https://web.mit.edu/~r/current/arch/i386_linux26/lib/R/library/GenomeInfoDb/html/seqinfo.html
## ---------------------------------------------------------------------
## 重命名 TxDb 对象的染色体信息
## ---------------------------------------------------------------------
# 比如
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
seqlevels(txdb)
# "chr2L" "chr2R" "chr3L" "chr3R" "chr4" "chrX" "chrU" "chrM" "chr2LHet" "chr2RHet" "chr3LHet" "chr3RHet" "chrXHet" "chrYHet" "chrUextra"
# 使用 `sub()` 函数将染色体前缀 `chr` 替换为空，即 `chr1` 变成 `1`
seqlevels(txdb) <- sub("chr", "", seqlevels(txdb))
seqlevels(txdb)
# [1] "2L" "2R" "3L" "3R" "4" "X" "U" "M" "2LHet" "2RHet" "3LHet" "3RHet" "XHet" "YHet" "Uextra"
# 使用 `paste0` 函数体加前缀，即 `1`变成 `CH1`
seqlevels(txdb) <- paste0("CH", seqlevels(txdb))
seqlevels(txdb)
# [1] "CH2L" "CH2R" "CH3L" "CH3R" "CH4" "CHX" "CHU" "CHM" "CH2LHet" "CH2RHet" "CH3LHet" "CH3RHet" "CHXHet" "CHYHet" "CHUextra"
# 修改某一条染色体名称，比如这里将染色体 `CHM` 改为 `M`
seqlevels(txdb)[seqlevels(txdb) == "CHM"] <- "M"
seqlevels(txdb)
# [1] "CH2L" "CH2R" "CH3L" "CH3R" "CH4" "CHX" "CHU" "M" "CH2LHet" "CH2RHet" "CH3LHet" "CH3RHet" "CHXHet" "CHYHet" "CHUextra"
# 使用 `seqlevels0()` 函数恢复原来的染色体水平
seqlevels(txdb) <- seqlevels0(txdb)
seqlevels(txdb)
# [1] "chr2L" "chr2R" "chr3L" "chr3R" "chr4" "chrX" "chrU" "chrM" "chr2LHet" "chr2RHet" "chr3LHet" "chr3RHet" "chrXHet" "chrYHet" "chrUextra"
## ---------------------------------------------------------------------
##
## ---------------------------------------------------------------------
tx <- transcripts(txdb)
seqlevels(tx)
# [1] "chr2L" "chr2R" "chr3L" "chr3R" "chr4" "chrX" "chrU" "chrM" "chr2LHet" "chr2RHet" "chr3LHet" "chr3RHet" "chrXHet" "chrYHet" "chrUextra"
# 丢弃 `chrM` 染色体，保留其它的染色体
seqlevels(tx, pruning.mode="coarse") <- seqlevels(tx)[seqlevels(tx) != "chrM"]
seqlevels(tx)
# [1] "chr2L" "chr2R" "chr3L" "chr3R" "chr4" "chrX" "chrU" "chr2LHet" "chr2RHet" "chr3LHet" "chr3RHet" "chrXHet" "chrYHet" "chrUextra"
# 除了 `ch3L` 和 `ch3R` 染色体保留，其余的都丢掉
seqlevels(tx, pruning.mode = "coarse") <- c("ch3L", "ch3R")
seqlevels(tx)
# [1] "ch3L" "ch3R"
## 番外详情见开始链接，这个知识点必须掌握，特别是有的是支持`chr1`类的，有的是支持`1`类的，就涉及要转换，当然其它方法也可以。
接 Peak 注释
# 查看有多少个peak
cat(paste0('there are ',length(peak),' peaks for this data' ))
# 使用 annotatPeak 进行注释，
peakAnno <- annotatePeak(peak, tssRegion = c(-3000, 3000),
TxDb = txdb, annoDb = "org.Dm.eg.db")
# 转变成 data.frame 格式文件，方便查看与后续操作
peakAnno_df <- as.data.frame(peakAnno)
# 获取样本名
sampleName = basename(strsplit(bedPeaksFile,'\\.')[[1]][1])
print(sampleName)
# [1] "Cg_WT"
## ---------------------------------------------------------------------
## 这里我拆开给大家一步一步看，大神可以直接跳过这里
## ---------------------------------------------------------------------
bedPeaksFile
# [1] "./oldBedFiles/Cg_WT.narrowPeak.bed"
# 使用 strsplit() 函数对 bedPeakFile 中的内容用 '.' 来切割
strsplit(bedPeaksFile,'\\.')
#[[1]]
#[1] "" "/oldBedFiles/Cg_WT" "narrowPeak" "bed"
# 然后我们要提取 "/oldBedFiles/Cg_WT" 信息
strsplit(bedPeaksFile,'\\.')[[1]][2]
# [1] "/oldBedFiles/Cg_WT"
# 然后通过 `basename()` 函数去掉前面的路径信息，即获得样本名
basename(strsplit(bedPeaksFile,'\\.')[[1]][2])
# [1] "Cg_WT"
对 peak 进行可视化，可以看到大部分 peak 都位于启动子区。
plotAnnoPie(peakAnno)
plotAnnoBar(peakAnno)
但是我们可以看到文章中的图是有好几个 sample 在一起的，而且注释的分布类也没这么细致，所以我们需要进行批量的 Peak 注释，并将其绘制出来
anno_bed <- function(bedPeaksFile){
peak <- readPeakFile( bedPeaksFile ) # 第一步肯定是读取 peak 文件
keepChr = !grepl('Het',seqlevels(peak)) # 提取不包含 Het 的染色体
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr] # 第三步就是过滤包含 Het 的染色体
cat(paste0('there are ',length(peak),' peaks for this data' )) # 第四步查看此 peak 文件有多少行
peakAnno <- annotatePeak(peak, tssRegion = c(-3000, 3000),
TxDb = txdb, annoDb = "org.Dm.eg.db") # 第五步对 Peak 进行注释
peakAnno_df <- as.data.frame(peakAnno) # 第六步将 peakAnno 文件转换为 dataframe 格式文件，方便操作
sampleName = basename(strsplit(bedPeaksFile,'\\.')[[1]][1]) # 第七步提取此 peak 的样本信息
return(peakAnno_df)
}
# 可以看到就是替换掉了前面所导入的 Peak 文件变量，比较知识换一下就好了
require(ChIPseeker)
# BiocManager::install("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
# BiocManager::install("org.Dm.eg.db")
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene )
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
require(clusterProfiler)
# 得到 oldBedFiles 文件目录下所有 WT 的 peak 文件，构成一个 list
f = list.files(path = 'oldBedFiles/',pattern = 'WT',full.names = T)
f
[1] "oldBedFiles/Cg_WT.narrowPeak.bed" "oldBedFiles/Ez_WT.narrowPeak.bed" "oldBedFiles/Pc_WT.narrowPeak.bed" "oldBedFiles/Ph_WT.narrowPeak.bed.bed" "oldBedFiles/Pho_WT.narrowPeak.bed" "oldBedFiles/Psc_WT.narrowPeak.bed" "oldBedFiles/Spps_WT.narrowPeak.bed"
# 结合上面我们编写的函数使用 lapply 进行批量注释
# 这里需要好好学习 apply() 家族函数的学习
# lapply 可以用 list
# 可以看到 Rstudio 右上方 tmp 文件是含有 7 的
tmp = lapply(f, anno_bed)
# 提取 list 中的第一个文件进行查看
# 我们可以清楚看到有 Peak 的坐标信息，以后注释到的基因的位置信息以及基因名，距离基因 TSS 有多远等。
head(tmp[[1]])
seqnames start end width strand V4 V5 V6 V7 V8 V9 V10 annotation geneChr geneStart geneEnd
1 chr2L 5104 5375 272 * Cg_WT_peak_1 45 . 2.29299 5.95519 4.56517 19 Promoter (2-3kb) 1 7529 9484
2 chr2L 5544 5962 419 * Cg_WT_peak_2 268 . 3.91194 28.67657 26.81658 194 Promoter (1-2kb) 1 7529 9484
3 chr2L 16625 17115 491 * Cg_WT_peak_3 953 . 8.46518 97.93490 95.38403 171 Promoter (1-2kb) 1 9839 18570
4 chr2L 17743 19205 1463 * Cg_WT_peak_4 1011 . 6.21549 103.78052 101.17944 887 Promoter (<=1kb) 1 9839 18570
5 chr2L 19439 19789 351 * Cg_WT_peak_5 654 . 6.04155 67.72665 65.43917 208 Promoter (<=1kb) 1 9839 18570
6 chr2L 20538 21944 1407 * Cg_WT_peak_6 1413 . 9.20858 144.25723 141.30544 889 Promoter (<=1kb) 1 9839 21376
geneLength geneStrand geneId transcriptId distanceToTSS ENTREZID SYMBOL GENENAME
1 1956 1 FBgn0031208 FBtr0300689 -2154 33155 CG11023 uncharacterized protein
2 1956 1 FBgn0031208 FBtr0300689 -1567 33155 CG11023 uncharacterized protein
3 8732 2 FBgn0002121 FBtr0078170 1455 33156 l(2)gl lethal (2) giant larvae
4 8732 2 FBgn0002121 FBtr0078170 0 33156 l(2)gl lethal (2) giant larvae
5 8732 2 FBgn0002121 FBtr0078171 -869 33156 l(2)gl lethal (2) giant larvae
6 11538 2 FBgn0002121 FBtr0078166 0 33156 l(2)gl lethal (2) giant larvae
# 然而我们这里绘制此图可以看出只需要 annotation 那一列的信息。
# 首先我们来统计一个 Peak 的注释分布数目
test = tmp[[1]]
num1 = length(grep('Promoter', test$annotation))
num2 = length(grep("5' UTR", test$annotation))
num3 = length(grep('Exon', test$annotation))
num4 = length(grep('Intron', test$annotation))
num5 = length(grep("3' UTR", test$annotation))
num6 = length(grep('Intergenic', test$annotation))
c(num1,num2,num3,num4,num5,num6 )
# [1] 12250 32 99 3291 148 1999
# 我们可以看到对于不同的 Peak，只需要换一下名字即可，那么我们就可以写成函数
df = do.call(rbind,lapply(tmp, function(x){
#table(x$annotation)
num1 = length(grep('Promoter',x$annotation))
num2 = length(grep("5' UTR",x$annotation))
num3 = length(grep('Exon',x$annotation))
num4 = length(grep('Intron',x$annotation))
num5 = length(grep("3' UTR",x$annotation))
num6 = length(grep('Intergenic',x$annotation))
return(c(num1,num2,num3,num4,num5,num6 ))
}))
df
# [,1] [,2] [,3] [,4] [,5] [,6]
#[1,] 12250 32 99 3291 148 1999
#[2,] 7057 15 63 1027 79 1113
#[3,] 3871 12 69 763 37 1110
#[4,] 10558 34 92 2819 67 1809
#[5,] 6776 9 55 1436 40 1343
#[6,] 10085 19 79 2255 46 1669
# [7,] 3691 5 18 507 11 410
# 赋值列名
colnames(df) = c('Promoter',"5' UTR",'Exon','Intron',"3' UTR",'Intergenic')
# 加上行名，即样本名
# 这里好好理解以下 unlist 函数
rownames(df) = unlist(lapply(f, function(x){
basename(strsplit(x,'\\.')[[1]][1])
}))
df
# Promoter 5' UTR Exon Intron 3' UTR Intergenic
#Cg_WT 12250 32 99 3291 148 1999
#Ez_WT 7057 15 63 1027 79 1113
#Pc_WT 3871 12 69 763 37 1110
#Ph_WT 10558 34 92 2819 67 1809
#Pho_WT 6776 9 55 1436 40 1343
#Psc_WT 10085 19 79 2255 46 1669
#Spps_WT 3691 5 18 507 11 410
# 由于是百分比，我们可以再转换成百分比数据，其实也可以不转
df1 <- apply(df, 1, function(x) x/sum(x))
df1
# Cg_WT Ez_WT Pc_WT Ph_WT Pho_WT Psc_WT Spps_WT
# Promoter 0.687468433 0.754436605 0.660354828 0.686520580 0.7015218967 0.712569773 0.795131409
# 5' UTR 0.001795836 0.001603592 0.002047083 0.002210807 0.0009317735 0.001342472 0.001077122
# Exon 0.005555867 0.006735087 0.011770727 0.005982183 0.0056941712 0.005581855 0.003877639
# Intron 0.184690499 0.109792602 0.130160355 0.183301905 0.1486696345 0.159330177 0.109220164
# 3' UTR 0.008305741 0.008445585 0.006311839 0.004356590 0.0041412154 0.003250194 0.002369668
# Intergenic 0.112183624 0.118986530 0.189355169 0.117627934 0.1390413086 0.117925528 0.088323998
好了，到目前位置，我们得到了绘制此图所需要的数据，接下来进行可视化
# 将数据转换成 ggpurb 所需要的格式
library(reshape2)
df2 = melt(df1)
colnames(df2) = c('dis','sample','fraction')
head(df2)
# dis sample fraction
#1 Promoter Cg_WT 0.687468433
#2 5' UTR Cg_WT 0.001795836
#3 Exon Cg_WT 0.005555867
#4 Intron Cg_WT 0.184690499
#5 3' UTR Cg_WT 0.008305741
#6 Intergenic Cg_WT 0.112183624
library(ggpubr)
ggbarplot(df2, "sample", "fraction",
fill = "dis", color = "dis", palette = "jco" )
想和文章一样的配色？还记得在重复第一个图的时候提到的 Takecolor ? 直接吸就完事
# 吸到的颜色。。。
my_color = c("#DD2A18", # 红色
"#41678B", # 蓝色
"#42A441", # 绿色
"#8C4498", # 紫色
"#FE7400", # 橙色
"#FEFE2D") # 黄色
ggbarplot(df2, "sample", "fraction",
fill = "dis", color = "black", palette = my_color ) +
guides(fill = guide_legend(nrow = 1))
最后的最后，虽然有些地方和文章不同，但是此文仅提供思想，比如这里的注释到底怎么划分是要按照大家各自的需求来的。所以参考就好。ggpurb 虽然便捷，但是如果你用 R 用的频繁，那么请您好好学习 ggplot2 或者 baseplot 等。
最后按照惯例，我们应该是有生信技能树的友情推广，但是我们长沙站招生已经满了，而且未来的两个月并没有外出巡讲计划，所以就不宣传了，大家可以继续添加小助手，说出自己的城市需求，我们会优先安排明年巡讲场次：广州专场（全年无休）GEO数据挖掘课，带你飞（1.11-1.12）和 生信入门课全国巡讲2019收官--长沙站