m6ASeq RNA甲基化测序

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

2021-01-27  本文已影响0人  信你个鬼

文章信息

文献标题: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覆盖度高,在区域覆盖度 > 104 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主要有两点:

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分析流程后续会陆陆续续更新,请期待吧。

上一篇下一篇

猜你喜欢

热点阅读