孟佳课题组|Peak Calling软件exomePeak以及exomePeak2使用指北
前面我提到了有学员提问： 新的ngs流程该如何学习之m6A学习大纲，所以我给了他优秀学员想要MeRIP-Seq/m6A-seq教程，但是他并没有主动给我投稿，不过我们的转录组授课讲师也感兴趣了这个技术，想丰富一下课程内容，所以开启了这个系列学习哈！
使用这个软件的原因主要在于它是专门针对6a-Seq数据的甲基化处理软件。
这款软件输入样本的bam文件，一个函数就可以完成peak calling以及差异peak分析，可以说是非常简单方便了；缺点就是软件自发表后，更新维护非常少，发表完了就完事了。引用率还行，目前很多可用的专门针对RNA 甲基化分析的软件和包，大多数都会借鉴一部分ATAC-Seq和ChIP-Seq分析pipeline。
# 安装软件，服务器 conda环境
conda install -c bioconda bioconductor-exomepeak -y
# 使用示例
rm(list=ls())
options(stringsAsFactors = F)
#devtools::install_github("ZW-xjtlu/exomePeak")
library(exomePeak)
gtf <- system.file("extdata", "example.gtf", package="exomePeak")
f1 <- system.file("extdata", "IP1.bam", package="exomePeak")
f2 <- system.file("extdata", "IP2.bam", package="exomePeak")
f3 <- system.file("extdata", "IP3.bam", package="exomePeak")
f4 <- system.file("extdata", "IP4.bam", package="exomePeak")
f5 <- system.file("extdata", "Input1.bam", package="exomePeak")
f6 <- system.file("extdata", "Input2.bam", package="exomePeak")
f7 <- system.file("extdata", "Input3.bam", package="exomePeak")
result <- exomepeak(GENE_ANNO_GTF=gtf, IP_BAM=c(f1,f2,f3,f4), INPUT_BAM=c(f5,f6,f7))
recommended_peaks <- result$con_peaks # consistent peaks (Recommended!)
peaks_info <- mcols(recommended_peaks) # information of the consistent peaks
head(peaks_info)
## to get all the peak detected (some of them do not consistently appear on all replicates):
all_peaks <- result$all_peaks # get all peaks
peaks_info <- mcols(all_peaks) # information of all peaks
head(peaks_info)
## 差异peak分析
## Peak Calling and Differential Methylation Analysis
gtf <- system.file("extdata", "example.gtf", package="exomePeak")
f1 <- system.file("extdata", "IP1.bam", package="exomePeak")
f2 <- system.file("extdata", "IP2.bam", package="exomePeak")
f3 <- system.file("extdata", "IP3.bam", package="exomePeak")
f4 <- system.file("extdata", "IP4.bam", package="exomePeak")
f5 <- system.file("extdata", "Input1.bam", package="exomePeak")
f6 <- system.file("extdata", "Input2.bam", package="exomePeak")
f7 <- system.file("extdata", "Input3.bam", package="exomePeak")
f8 <- system.file("extdata", "treated_IP1.bam", package="exomePeak")
f9 <- system.file("extdata", "treated_Input1.bam", package="exomePeak")
# 本地指定gtf文件
result <- exomepeak(GENE_ANNO_GTF=gtf,
IP_BAM=c(f1,f2,f3,f4),
INPUT_BAM=c(f5,f6,f7),
TREATED_IP_BAM=c(f8),
TREATED_INPUT_BAM=c(f9))
diff_peaks <- result$diff_peaks # consistent differential peaks (Recommended!)
peaks_info <- mcols(diff_peaks) # information of the consistent peaks
head(peaks_info[,1:3]) # peak calling information
head(peaks_info[,4:6]) # differential analysis information
## 联网使用gtf文件，推荐上面的本地指定
result <- exomepeak(GENOME="hg19",
IP_BAM=c(f1,f2,f3,f4),
INPUT_BAM=c(f5,f6,f7),
TREATED_IP_BAM=c(f8),
TREATED_INPUT_BAM=c(f9))
结果目录，每个结果详细解释请参考github：
一、github上的文档：https://github.com/ZW-xjtlu/exomePeak
有详细的结果描述，有p值说明，有FDR值，但是没有说FDR使用何种方法。需要去看核心代码。
二、bioconductor上的示例：https://bioconductor.riken.jp/packages/3.1/bioc/html/exomePeak.html
与github上一样，但是示例中IP的样本数（4个）与INPUT的样本数（三个）不一样，难道不需要确定并保证每一个IP样本都有与之对应的Input样本做对照参考识别m6A么？说明书还有另外一个函数可以使用。需要去看核心代码。
三、发表的相应的文章：doi:10.1016/j.ymeth.2014.06.008
1.文章中的示例：
需要去重复一下该示例并与使用数据的结果进行比较
2.文章中讨论了比对部分
认为，多比对的reads应该去除，认为是来自参考基因组中的repeat element（占50%的比例）。但是我发现我的m6A数据多比对率非常高，如果去除，那么剩下的有效数据太少了，我需要确认这些多比对的reads都比对到参考基因的什么地方。多比对率如下：
11065604 (53.68%) aligned >1 times
5432119 (75.88%) aligned >1 times
22413822 (73.37%) aligned >1 times
5764062 (19.73%) aligned >1 times
24660269 (73.18%) aligned >1 times
22220016 (81.88%) aligned >1 times
14252127 (51.49%) aligned >1 times
25307217 (85.24%) aligned >1 times
文章中的说明部分：
使用感受：由于是R语言，因此真实项目数据跑起来耗时比较久；其次，peak calling与diff peak是两个独立的过程，结果会有一些不一致的地方。
后来，作者又出了一个exomePeak2，并不是在原有的exomePeak基础上进行的更新。
主页：https://rdrr.io/github/ZhenWei10/exomePeak2/man/exomePeak2.html
exomePeak2与exomePeak最大的区别在于exomePeak2可以输出中间结果，应该还有其他不一样的地方，来看看吧。
功能：exomePeak2使用bam文件进行peak calling以及peak统计，它整合了meRIP-seq数据分析的常规分析内容：
1.使用scanMeripBAM检查BAM的index 2.使用exomePeakCalling识别外显子区域的被修饰的peaks 3.normalizeGC计算GC偏倚 4.glmM或glmDM构建线性模型来计算差异位点 5.exportResults输出peak结果
看一下函数
exomePeak2(
bam_ip = NULL,
bam_input = NULL,
bam_treated_ip = NULL,
bam_treated_input = NULL,
txdb = NULL,
bsgenome = NULL,
genome = NA,
gff_dir = NULL,
mod_annot = NULL,
paired_end = FALSE,
library_type = c("unstranded", "1st_strand", "2nd_strand"),
fragment_length = 100,
binding_length = 25,
step_length = binding_length,
peak_width = fragment_length/2,
pc_count_cutoff = 5,
bg_count_cutoff = 50,
p_cutoff = 1e-05,
p_adj_cutoff = NULL,
log2FC_cutoff = 1,
parallel = FALSE,
background_method = c("all", "Gaussian_mixture", "m6Aseq_prior", "manual"),
manual_background = NULL,
correct_GC_bg = TRUE,
qtnorm = FALSE,
glm_type = c("DESeq2", "Poisson", "NB"),
LFC_shrinkage = c("apeglm", "ashr", "Gaussian", "none"),
export_results = TRUE,
export_format = c("CSV", "BED", "RDS"),
table_style = c("bed", "granges"),
save_plot_GC = TRUE,
save_plot_analysis = FALSE,
save_plot_name = "",
save_dir = "exomePeak2_output",
peak_calling_mode = c("exon", "full_tx", "whole_genome")
)
与exomePeak相比，增了加双端数据paired_end ，文库类型library_type ，count阈值pc_count_cutoff ，bg_count_cutoff 等参数。还有一些输出结果的更改。
没有了那个要求在group中call peak一致性的参数。这个参数应该是在分步计算里面。
在运行之前，示例代码使用USCS中的注释文件，最好先安装好那个注释包
BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")
library(BSgenome.Hsapiens.UCSC.hg19)
指定进行peak calling的样本，
rm(list = ls())
options(stringsAsFactors = F)
#BiocManager::install("exomePeak2")
library(exomePeak2)
# 设置随机种子
set.seed(1)
GENE_ANNO_GTF = system.file("extdata", "example.gtf", package="exomePeak2")
f1 = system.file("extdata", "IP1.bam", package="exomePeak2")
f2 = system.file("extdata", "IP2.bam", package="exomePeak2")
f3 = system.file("extdata", "IP3.bam", package="exomePeak2")
f4 = system.file("extdata", "IP4.bam", package="exomePeak2")
IP_BAM = c(f1,f2,f3,f4)
f1 = system.file("extdata", "Input1.bam", package="exomePeak2")
f2 = system.file("extdata", "Input2.bam", package="exomePeak2")
f3 = system.file("extdata", "Input3.bam", package="exomePeak2")
INPUT_BAM = c(f1,f2,f3)
分步骤计算过程，大致分为7个步骤：
## Peak Calling and Visualization in Multiple Steps
#The exomePeak2 package can achieve peak calling and peak statistics calulation with multiple functions.
## 1. 检查bam文件
MeRIP_Seq_Alignment <- scanMeripBAM(
bam_ip = IP_BAM,
bam_input = INPUT_BAM,
paired_end = FALSE)
# 同时含有处理组和非处理组
MeRIP_Seq_Alignment <- scanMeripBAM(
bam_ip = IP_BAM,
bam_input = INPUT_BAM,
bam_treated_input = TREATED_INPUT_BAM,
bam_treated_ip = TREATED_IP_BAM,
paired_end = FALSE)
# str函数可以方便快速的查看s4对象的结构和内容
str(MeRIP_Seq_Alignment,max.level = 3)
#2. 使用bam文件进行 peak calling
SummarizedExomePeaks <- exomePeakCalling(merip_bams = MeRIP_Seq_Alignment,
gff_dir = GENE_ANNO_GTF,
genome = "hg19")
#可选，用来评估数据
SummarizedExomePeaks <- exomePeakCalling(merip_bams = MeRIP_Seq_Alignment,
gff_dir = GENE_ANNO_GTF,
genome = "hg19",
mod_annot = MOD_ANNO_GRANGE)
# 查看peak结果
str(SummarizedExomePeaks, max.level = 4)
#3. 计算size factors用来对GC偏倚进行矫正
SummarizedExomePeaks <- normalizeGC(SummarizedExomePeaks)
#4. 使用glmM构造peak统计量
SummarizedExomePeaks <- glmM(SummarizedExomePeaks)
# 可选，如果有差异分析，就分析此步骤
SummarizedExomePeaks <- glmDM(SummarizedExomePeaks)
#5. Generate the scatter plot between GC content and log2 Fold Change (LFC).
p <- plotLfcGC(SummarizedExomePeaks)
# 点的大小有点小，取出来数据重新加工
library(ggplot2)
data <- p$data
head(data)
p1 <- ggplot(data,aes(x=GC_idx,y=Log2FC,color=Label)) + geom_point(size=2) + theme_classic()
p1
#6. Generate the bar plot for the sequencing depth size factors.
plotSizeFactors(SummarizedExomePeaks)
#7. Export the modification peaks and the peak statistics with user decided format.
exportResults(SummarizedExomePeaks, format = "BED")
GC含量散点图与SizeFactor图
使用示例结果中的一个差异peak在IGV中查看peak分布
相对于exomePeak，作者还输出了每个peak在每个样本中的read count数，就用第一个peak示例吧，差异FC比较大，p值也显著。peak大概100个bp这么长。
至于具体效果如何，大家自行判断吧。
请等待后续更新。(因为转录组讲师有自己的工作，所以这个m6a系列没办法更新那么快，而且也不清楚受众多不多，如果你比较喜欢这方面教程，可以文末留言催更哈！)