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

image

• 第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_cutoff a 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区间的相关基因注释，但是没有功能区域注释，即所有文章中最常见的这幅图：

image

对于左图，我们这里给大家介绍绘制m6A修饰位点在基因组上的分布特征图可视化包：Guitar包

Guitar包与exomePeak2开发者来自同一个课题组，熟悉这个课题组的应该知道他们在关于m6A的分析上开发了许多的包。

这个包目前维护比较少，为了避免踩坑，建议安装"2.6.0"的版本。

Guitar包出来的结果特征如下：

image

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上。

image

对于右图，一般用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()

注释结果如下：

image

如果不喜欢这种注释，当然还可以直接用bedtools进行注释~

下一期分享另类可视化方法，如何与文献中的图片一样好看！