m6A-Seq数据质量评估：trumpet包
前面我提到了有学员提问：新的ngs流程该如何学习之m6A学习大纲，所以我给了他优秀学员想要MeRIP-Seq/m6A-seq教程，但是他似乎是并没有主动给我投稿，不过我们的转录组授课讲师也感兴趣了这个技术，想丰富一下课程内容，所以开启了这个系列学习哈！
文章信息
文献标题：trumpet: transcriptome-guided quality assessment of m6A-seq data
发表杂志：BMC Bioinformatics
发表时间：13 July 2018
原文doi：https://doi.org/10.1186/s12859-018-2266-3
关于作者-Jia Meng
最近做m6A分析研发，发现搜到的软件都是他们家的，就去搜了一下这个大佬
主页介绍如下：https://www.xjtlu.edu.cn/en/departments/academic-departments/biological-sciences/staff/jia-meng
实验室主页：https://www.xjtlu.edu.cn/zh/news/2019/02/m6a
其中关于Peak识别发表的软件一部分如下：
可以说是非常高产了，其中exomePeak的升级版本预计今年5月份见刊。
文献摘要
RNA甲基化免疫沉淀测序(MeRIP-seq或m6A-seq)已被广泛用于分析RNA N6-腺苷甲基化在转录组中的分布。然而，由于RNA分子的固有特性以及该技术复杂的操作过程，m6A-seq数据往往存在各种缺陷。对m6A-seq数据的质量进行评估需要一种方便、全面的工具，以确保它们适合后续的分析。
从技术方面，m6A-seq可以认为是ChIP-Seq和RNA-Seq的结合。因此，通过有效地结合两种技术的数据质量评估指标，我们开发了用于m6A-seq数据质量评估的trumpet R包。trumpet包从m6A-seq数据中获取比对产生的BAM文件以及转录组信息作为输入，生成HTML格式的质量评估报告。
除此之外，还可以用于其他RNA免疫沉淀测序技术数据评估如m1A-seq, CeU-Seq, Ψ-seq等。
主要评价指标
1.测序数据统计
这个地方主要通过计算read count来获得对样本的一个全面了解，这可能是检查样本质量的最基本方法。低reads count或比对到特定基因组区域的reads比例差异过大可能与低数据质量有关，这是由于多样本混库测序不平衡、DNA污染或实验过程中的其他偏差造成的。
下表来源于小鼠中脑3个Fto基因敲除和3个对照样本的统计结果（GSE47217），其中IP2样本3'UTR的reads相比其他样本过少，可能是由于样本制备过程中 3′ bias 造成的。
2.read coverage评价
由于基因差异表达， PCR artifacts和randomness，可能会造成全转录组read覆盖度的异质性。下表中，相对于其他样本，IP2样本的外显子区域read覆盖度高，在区域覆盖度 > 10^4^ reads 时更显著。这表明read覆盖度的异质性很高，这可能表明在样品制备或测序过程中存在潜在的PCR rtifacts。
这个现象在 FASTQC评估结果中也得到了验证：IP2 has highest Kmer content among the samples (fold enrichment of the most over-represented Kmer: 33.31 in IP2 vs 23.61 and 27.63 in IP1 and IP3)。
PCR artifacts 可能进一步加剧了m6 A-seq实验reads覆盖的异质性。
3.reads分布可视化
RNA m6A一般分布在起始密码子和终止密码子附近。这个结果展示了(5’UTR, CDS and 3’UTR)这三个区域的m6A分布。
4.使用ESES评估免疫沉淀反应效率
m6A-Seq数据的一个主要评价指标就是免疫沉淀反应效率，只要体现在免疫沉淀信号的富集程度。
为了评估IP样本中的m6A信号，trumpet包使用 ESEC：exome signal extraction scaling这个指标。该指来源于评估ChIP-seq数据的信号SES方法，ESES不同于SES主要有两点：
首先，将SES方法中ChIP-seq数据的基因组背景替换为标准化的MeRIP-seq数据的基因特异性外显子背景，排除不携带有意义信号的区域(内含子和非基因区域)的影响。 其次，根据基因的表达水平对MeRIP-seq数据的read覆盖进行归一化处理，以消除基因不同表达水平的影响。
我们可以看到第二个IP样本(IP2)与其他样本有很大的不同，这与之前的结果是一致的。
5.使用C-test评估m6A信号富集程度
此指标也显示IP2样本异常，与之前的评估结果一致。
6.对样本进行层次聚类和PCA分析
我感觉这个结果有点充数了。。。主要用于样本的一致性和可重复性评估。
7.特异基因的read覆盖异质性
在实际项目中，bam中比对的reads并不是均匀地分布在同一个基因上。IP 样本中read覆盖度的异质性可能来源于m6A位点信号的富集，
坐标轴分别表示：每个样本中每一个gene的read count的均值和标准差。与IP样本相比，input样本的这个曲线相对更加平缓。
实现代码
github：https://github.com/skyhorsetomoon/Trumpet
使用特别简单，输入样本的bam文件和一个gtf文件就行了
## All gone with one touch
rm(list=ls())
options(stringsAsFactors = F)
## 安装
#devtools::install_github("skyhorsetomoon/Trumpet")
library(Trumpet)
## 使用R包中的示例数据
f1 <- system.file("extdata", "IP1.bam", package="Trumpet")
f2 <- system.file("extdata", "IP2.bam", package="Trumpet")
f3 <- system.file("extdata", "IP3.bam", package="Trumpet")
f4 <- system.file("extdata", "IP4.bam", package="Trumpet")
f5 <- system.file("extdata", "Input1.bam", package="Trumpet")
f6 <- system.file("extdata", "Input2.bam", package="Trumpet")
f7 <- system.file("extdata", "Input3.bam", package="Trumpet")
f8 <- system.file("extdata", "treated_IP1.bam", package="Trumpet")
f9 <- system.file("extdata", "treated_Input1.bam", package="Trumpet")
ip_bam <- c(f1,f2,f3,f4)
input_bam <- c(f5,f6,f7)
contrast_ip_bam <- c(f8)
contrast_input_bam <- c(f9)
## gtf文件
gtf <- system.file("extdata", "hg19toy.gtf", package="Trumpet")
## 生成HTML报告
trumpet_report <- Trumpet_report(IP_BAM = ip_bam,
Input_BAM = input_bam,
contrast_IP_BAM = contrast_ip_bam,
contrast_Input_BAM = contrast_input_bam,
condition1 = "untreated",
condition2 = "treat",
GENE_ANNO_GTF = gtf)
分析走到这里，出现了一个报错：
上网查找原因，是因为这个trumpet包需要使用到Guitar包，额，然后我发现这个Guitar包也是他们课题组发的。这个包后来更新了，现在更新到了2.6.0，这个函数名字就变成了另外一个名字。这个问题也有人在他们新开发的exomePeak2里面提到了，见：https://github.com/ZW-xjtlu/exomePeak2/issues/1。
然后我想，那我用老版本的吧，找到了1.7.0版本：https://mirror.nju.edu.cn/bioconductor/3.2/bioc/html/Guitar.html，然后本地安装：
install.packages("E:/01-m6A/m6A分析流程搭建/code/Guitar_1.7.0.tar.gz", repos = NULL, type = "source")
发现还是不行，就联系了作者，作者给我了一个1.5.0的版本的：https://github.com/lzcyzm/Guitar，重新安装：
#devtools::install_github("lzcyzm/Guitar")
library(Guitar)
package.version("Guitar")
[1] "1.5.0"
然后顺利跑上了前面的代码，生成了质控报告。可以说非常曲折了。这个1.5.0版本的，如果不是作者给我链接，我感觉我还真找不到。这个包由于自发布之后，更新比较少，希望我们将这个包推荐给大家使用后，大家能多多给反馈，希望作者也能后续更新以便更多的人能使用上。
生成的网页报告结果还是很不错的。
下面来看一眼生成的网页版本的m6A-Seq质控报告吧，默认在输出路径下生成一个Trumpet_report.html。
报告上方有目录链接，其实从这个地方，我们也可以学习作者的源代码，看看怎么使用R语言生成html报告。
有比对结果统计：
read分布可视化：m6A-Seq的peak一般在5'UTR和终止密码子附近有富集【doi:10.1038/nature11112 ，doi:10.1016/j.cell.2012.05.003】。
还有一些其他比较有意思的图，大家可以自行去看看啦。
m6A分析流程后续会陆陆续续更新，请期待吧。