m6A-Seq数据质量评估：trumpet包

文章信息

文献标题：trumpet: transcriptome-guided quality assessment of m6A-seq data
发表杂志：BMC Bioinformatics
发表时间：13 July 2018
原文doi： https://doi.org/10.1186/s12859-018-2266-3

image-20210125162326019.png

关于作者-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识别发表的软件一部分如下：

image-20210125162923656.png

可以说是非常高产了，其中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 造成的。

image-20210125165739190.png

2.read coverage评价

由于基因差异表达， PCR artifacts和randomness，可能会造成全转录组read覆盖度的异质性。下表中，相对于其他样本，IP2样本的外显子区域read覆盖度高，在区域覆盖度 > 10⁴ 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覆盖的异质性。

image-20210125170918589.png

3.reads分布可视化

RNA m6A一般分布在起始密码子和终止密码子附近。这个结果展示了(5’UTR, CDS and 3’UTR)这三个区域的m6A分布。

image-20210125172710463.png

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覆盖进行归一化处理，以消除基因不同表达水平的影响。

image-20210125175224898.png

我们可以看到第二个IP样本(IP2)与其他样本有很大的不同，这与之前的结果是一致的。

image-20210125175340959.png

5.使用C-test评估m6A信号富集程度

此指标也显示IP2样本异常，与之前的评估结果一致。

image-20210125175803640.png

6.对样本进行层次聚类和PCA分析

我感觉这个结果有点充数了。。。主要用于样本的一致性和可重复性评估。

image-20210125180006565.png

7.特异基因的read覆盖异质性

在实际项目中，bam中比对的reads并不是均匀地分布在同一个基因上。IP 样本中read覆盖度的异质性可能来源于m6A位点信号的富集，

坐标轴分别表示：每个样本中每一个gene的read count的均值和标准差。与IP样本相比，input样本的这个曲线相对更加平缓。

image-20210125180220455.png

实现代码

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)

分析走到这里，出现了一个报错：

image-20210127112720399.png

上网查找原因，是因为这个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("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报告。

image-20210127113720767.png

有比对结果统计：
这个地方直接用read数而不是M为单位的我感觉会更直观好一点。

image-20210127113855434.png

read分布可视化：
m6A-Seq的peak一般在5'UTR和终止密码子附近有富集【doi:10.1038/nature11112 ，doi:10.1016/j.cell.2012.05.003】。

image-20210127114043155.png

还有一些其他比较有意思的图，大家可以自行去看看啦。

m6A分析流程后续会陆陆续续更新，请期待吧。