使用這個(gè)軟件的原因主要在于前面使用MACS2進(jìn)行m6A peak calling框咙，下游沒有可用的現(xiàn)成的方法銜接此部分結(jié)果進(jìn)行差異peak分析询吴，可能得自己根據(jù)文獻(xiàn)想辦法使用什么統(tǒng)計(jì)學(xué)方法定義差異的peak了鬼吵，目前也考慮去搜集可以銜接這部分結(jié)果做差異peak的文獻(xiàn)方法项炼。
這款軟件輸入樣本的bam文件揭厚，一個(gè)函數(shù)就可以完成peak calling以及差異peak分析但绕，可以說是非常簡(jiǎn)單方便了救崔；缺點(diǎn)就是軟件自發(fā)表后，更新維護(hù)非常少捏顺，發(fā)表完了就完事了六孵。引用率還行，目前沒有很多可用的專門針對(duì)RNA 甲基化分析的軟件和包幅骄，大多數(shù)都會(huì)借鑒一部分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))
結(jié)果目錄，每個(gè)結(jié)果詳細(xì)解釋請(qǐng)參考github：
一拆座、github上的文檔：https://github.com/ZW-xjtlu/exomePeak
有詳細(xì)的結(jié)果描述主巍，有p值說明，有FDR值挪凑，但是沒有說FDR使用何種方法孕索。需要去看核心代碼。
二岖赋、bioconductor上的示例：https://bioconductor.riken.jp/packages/3.1/bioc/html/exomePeak.html
與github上一樣檬果，但是示例中IP的樣本數(shù)（4個(gè)）與INPUT的樣本數(shù)（三個(gè)）不一樣，難道不需要確定并保證每一個(gè)IP樣本都有與之對(duì)應(yīng)的Input樣本做對(duì)照參考識(shí)別m6A么唐断？說明書還有另外一個(gè)函數(shù)可以使用选脊。需要去看核心代碼。
三脸甘、發(fā)表的相應(yīng)的文章：doi:10.1016/j.ymeth.2014.06.008
1.文章中的示例：
需要去重復(fù)一下該示例并與使用數(shù)據(jù)的結(jié)果進(jìn)行比較
2.文章中討論了比對(duì)部分
認(rèn)為恳啥，多比對(duì)的reads應(yīng)該去除，認(rèn)為是來自參考基因組中的repeat element（占50%的比例）丹诀。但是我發(fā)現(xiàn)我的m6A數(shù)據(jù)多比對(duì)率非常高钝的，如果去除翁垂，那么剩下的有效數(shù)據(jù)太少了，我需要確認(rèn)這些多比對(duì)的reads都比對(duì)到參考基因的什么地方硝桩。多比對(duì)率如下：
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
文章中的說明部分：
然后看了聯(lián)川的報(bào)告沿猜，感覺給的一句解釋說了跟沒說一樣，并且表頭還有錯(cuò)誤碗脊。
英文報(bào)告來源：https://www.lcsciences.com/documents/sample_data/m6A_sequencing/m6A_seq_polyA_html_report.html#4.%20Results
中文版本：
后面繼續(xù)更新~
