m6A图文复现07-Peak结果以及分布特征图
上一期提到使用exomePeak2进行m6A的Peak Calling，结果如下，每个样本一个结果：
├── ADDInfo
│ ├── ADDInfo_GLM_allDesigns.csv：Peak识别过程中的一些模型统计参数值
│ ├── ADDInfo_ReadsCount.csv：每个Peak的count值
│ └── ADDInfo_RPKM.csv：每个peak的RPKM值
├── Mod.bed ：peak的bed格式
├── Mod.csv：peak的csv格式
├── Mod.rds：peak的Rdata格式
└── RunInfo.txt：运行过程中的一些参数文件
其中Mod.csv与Mod.bed前12列相同，bed为bed12，具体每一列如下，exomePeak2的结果不仅直接对每一个peak进行了基因注释，还输出了对应Peak的IP与Input样本的count值与RPKM值。筛选Peak主要根据pvalue或者padj。默认为pvalue<1e-05
第1列 chr：染色体编号
第2列 chromStart：Peak起始位点
第3列 chromEnd：Peak终止位点
第4列 name：Peak名字
第5列 score：the -log2 p value of the peak
第6列 strand：Peak在参考基因组上的链信息，+表示正链，-表示负链
第7列 thickStart：同第二列
第8列 thickEnd: 同第三列
第9列 itemRgb: the column for the RGB encoded color in BED file visualization.
第10列 blockCount: the block (exon) number within the peak.
第11列 blockSizes: the widths of the blocks.
第12列 blockStarts: the start positions of the blocks.
第13列 geneID: Peak区间内注释到的基因ID
第14列 ReadsCount.input：the total read count of input samples.
第15列 ReadsCount.IP: the total read count of IP samples.
这一列作为重点给予以下说明：
第16列 log2FoldChange: the estimate of IP over input log2 fold change (coefficient estimates of βi,1βi,1 in GLM), the Bayesian estimation implemented in the bioconductor package apeglm[3] will be returned if the regularized NB GLMs are fitted using DESeq2.
这一列很多人会有个误解，主要跟作者给的注释也有一定关系
包中的注释是这个样子的：
log2FC_cutoffa numeric value for the cutoff on log2 IP over input fold changes in peak calling; default = 0.
看着就是每一个Peak的IP样本中count/ Input样本中的count之后的log2值，即倍数关系
但是，当用户设置了这个参数，比如log2FC_cutoff=1的时候，发现结果中依然会有小于1的结果，后来给作者写信，给出答复如下，供大家参考：
回复1：
函数中的p_cutoff和log2FC_cutoff其实是peak calling时设定的阈值，exomePeak2的算法是这样的：
先在外显子组上生成划窗，对每个划窗收集到的IP / input 样本count做统计检验，并对其p-value (one-sided) 和 log2FC进行cut （即p_cutoff和log2FC_cutoff所指定的)，显著的划窗会被merge成为peak region，并再次count和计算统计值。
所以在exomePeak2最终输出的表格中，其实是基于merge的peak所计算的统计值，本身已经包含了一些positive bias了。而差异甲基化也是基于所有样本peak calling后得到的peak进行的统计。
回复2：
那个log2FC是peak calling时针对sliding window用的filter，而输出的结果是在peak 层面上重新计算的统计值，并未和sliding window共用filter。为了避免给后来使用者造成困惑，我在更新的版本中将默认的peak calling log2FC cutoff改成0了 （在peak calling中log2FC并不影响很大，主要的影响还是p value)。
第17列pvalue: the Wald test p value on the βi,1βi,1 coefficient in the GLM. 第18列padj: 对pvalue进行BH矫正后的fdr值
这个结果中已经有了Peak区间的相关基因注释，但是没有功能区域注释，即所有文章中最常见的这幅图：
对于左图，我们这里给大家介绍绘制m6A修饰位点在基因组上的分布特征图可视化包：Guitar包
Guitar包与exomePeak2开发者来自同一个课题组，熟悉这个课题组的应该知道他们在关于m6A的分析上开发了许多的包。
这个包目前维护比较少，为了避免踩坑，建议安装"2.6.0"的版本。
Guitar包出来的结果特征如下：
Figure 3: m6 A sites on mRNA and lncRNA. In mRNA, the strongest binding sites (“IP/input ratio” larger than 8) are highly enriched near stop codon side of 3 UTR and deficient on TSS (transcription starting site) side of 5 UTR and the phenomena are more prominent than lowly methylated sites. In contrast, the m6 A sites are almost uniformly distributed on lncRNA despite the “IP/input ratio” specified. Please note that, in this figure, the size of 5 UTR, CDS, and 3 UTR reflects their true width within the transcriptome, so the 5 UTR region is much shorter compared with the other two components. This result is based on peaks called on human HepG2 dataset [10].
参考：Biomed Research International, 28 Apr 2016, 2016:8367534
代码如下：
rm(list=ls())
options(stringsAsFactors = F)
# load package
library(Guitar)
package.version("Guitar")
# [1] "2.6.0"
stBedFiles <- list("data/KO1/Mod.bed")
txdb <- makeTxDbFromGFF(file = "data/KO1/Mus_musculus.GRCm38.101.gtf.gz",
format="gtf",
dataSource="Ensembl",
organism="Mus musculus")
p <- GuitarPlot(txTxdb = txdb,
stBedFiles = stBedFiles,
headOrtail = FALSE,
enableCI = FALSE,
mapFilterTranscript = TRUE,
pltTxType = c("mrna"),
stGroupName = "KO1")
p <- p + ggtitle(label = "Distribution on mRNA") +
theme(plot.title = element_text(hjust = 0.5)) + theme_bw()
png(file="data/KO1_guitarPlot_mRNA.pdf",width=8,height=6)
print(p)
dev.off()
png(filename="data/KO1_guitarPlot_mRNA.png",width = 1000,height = 800, res = 180)
print(p)
dev.off()
结果图如下：Peaks主要富集在终止密码子以及3'UTR上。
对于右图，一般用ChIPseeker软件来进行区间注释。
rm(list=ls())
options(stringsAsFactors = F)
#BiocManager::install("ChIPseeker",ask = F)
library(ChIPseeker)
library(org.Hs.eg.db)
library(clusterProfiler)
library(GenomicFeatures)
# annotation file gtf
txdb <- makeTxDbFromGFF(file = "data/KO1/Mus_musculus.GRCm38.101.gtf.gz", format="gtf",
dataSource="Ensembl", organism="Mus musculus")
# overlap
peak_file <- readPeakFile("data/KO1/Mod.bed")
## peak annotation
#region select："Promoter", "5UTR", "3UTR", "Exon", "Intron", "Downstream", "Intergenic"
peak_anno <- annotatePeak(peak_file,
tssRegion = c(-3000, 3000),
TxDb = txdb,
assignGenomicAnnotation = TRUE,
genomicAnnotationPriority = c("5UTR", "3UTR", "Exon","Intron","Intergenic"),
addFlankGeneInfo = TRUE,
flankDistance = 5000)
pdf(file = "KO1.Anno.Pie.pdf",width = 6,height = 5)
plotAnnoPie(peak_anno)
dev.off()
png(filename="KO1.Anno.Pie.png" ,width=1000, height=850, res=200)
plotAnnoPie(peak_anno)
dev.off()
注释结果如下：
如果不喜欢这种注释，当然还可以直接用bedtools进行注释~
下一期分享另类可视化方法，如何与文献中的图片一样好看！
文末友情推荐
做教学我们是认真的，如果你对我们的马拉松授课（直播一个月互动教学）有疑问，可以看完我们从2000多个提问互动交流里面精选的200个问答！ 2021第二期_生信入门班_微信群答疑整理，以及 2021第二期_数据挖掘班_微信群答疑笔记
与十万人一起学生信，你值得拥有下面的学习班：