exomePeak使用过程疑问

使用这个软件的原因主要在于前面使用MACS2进行m6A peak calling，下游没有可用的现成的方法衔接此部分结果进行差异peak分析，可能得自己根据文献想办法使用什么统计学方法定义差异的peak了，目前也考虑去搜集可以衔接这部分结果做差异peak的文献方法。

这款软件输入样本的bam文件，一个函数就可以完成peak calling以及差异peak分析，可以说是非常简单方便了；缺点就是软件自发表后，更新维护非常少，发表完了就完事了。引用率还行，目前没有很多可用的专门针对RNA 甲基化分析的软件和包，大多数都会借鉴一部分ATAC-Seq和ChIP-Seq分析pipeline。

# 安装软件
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 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")
​
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
​
## Download Gene Annotation Directly from Internet
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：

1610520090244.png

一、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

文章中的说明部分：

1610519542927.png

然后看了联川的报告，感觉给的一句解释说了跟没说一样，并且表头还有错误。

英文报告来源：https://www.lcsciences.com/documents/sample_data/m6A_sequencing/m6A_seq_polyA_html_report.html#4.%20Results

1610517880572.png

中文版本：

1610517980263.png

后面继续更新~