RNA-seq数据差异表达分析流程比较

赵延辉 陈少康 翟丽维 侍玉梅 原佳妮 盛熙晖 齐晓龙 郭 勇 王楚端 邢 凯*

(1.北京农学院 动物科学技术学院,北京 102206;2.北京市畜牧总站,北京 100107;3.中国农业大学 动物科学技术学院,北京 100193)

RNA-seq即转录组测序技术,能够全面快速地获得某一物种特定组织或器官在某一状态下的几乎所有转录本序列信息。随着高通量测序技术改进、成本下降,RNA-seq已经广泛应用于多物种的研究中。与DNA微阵列技术相比,RNA-seq数据可重复性强,背景噪声低,检测范围广泛[1],并且可以在没有参考基因组的情况下使用,包括从头测序[2]、丰度估计[3]、检测选择性剪切[4]等。RNA-seq所选工作流程不同,分析的准确性、速度、成本可能会大不相同[5]。

RNA-seq技术最常用的用途就是寻找特定条件下的差异表达基因,其工作流程首先将样本RNA片段化并且反转录成cDNA进行测序,将测序获得的短序列比对到参考基因组上,通过比对到基因组的reads数目来估算基因的相对表达水平,之后通过统计学方法检验组间基因差异表达[6]。目前,针对此流程已有许多工具被开发出来,如比对工具中TopHat2[7]、HISAT2[8]和STAR[9]被广泛使用,TopHat2可以将测序产生的各种长度的reads,即使是在高度重复的基因组或存在假基因的情况下,也能产生敏感和准确的比对。HISAT2是目前可用的最快系统,与其他工具相比其精度更好或相同,另外HISAT2需要计算机运行内存少,并且支持任何大小的基因组。STAR同样拥有良好的对齐精确度和灵敏度,同时STAR的运行速度也是其优势之一。DESeq2、edgeR和limma是3种常用的基因差异表达分析工具,DESeq2是Love等[10]在DESeq的基础上开发出的Bioconductor包,其对基因排序和可视化进行了改进,拥有较高的灵敏性和准确性,并且控制了假阳性率。EdgeR是由Mark等[11]开发的Bioconductor包,可以在没有生物学重复的情况下寻找差异表达基因。Limma最开始时只能处理来自微阵列的基因表达数据[12],自3.9.19版本开始可以用于分析RNA-seq数据。

少数研究比较了不同差异表达分析工具的性能[13-15],但是并未考虑基因组比对工具可能带来的影响。此外对于哪种分析工具是最佳的,以及如何确保分析结果的准确性和复现性,目前尚未达成共识。不同的工具进行组合产生不同的工作流程,不同的工作流程产生的结果可能大不相同。本研究比较了TopHat2、HISAT2和STAR 3种广泛使用的比对工具以及DESeq2、edgeR和limma 3种常用的基因差异表达分析工具,并综合KEGG通路分析结果对不同分析流程进行评价(图1),选择出合适的分析流程为后续研究工作提供参考。

图1 RNA-Seq数据差异表达分析工作流程Fig.1 The workflow of differential expression analysis for RNA-Seq data

1 材料与方法

1.1 转录组数据

本研究使用转录组数据来自参考文献[16],项目号为PRJNA234335和PRJNA287471。具体为选取松辽黑猪和长白猪各6头,食物充足、饮水自由、饲养条件相同。在屠宰场宰杀后取背部脂肪组织,并使用Trizol试剂盒提取12份背部脂肪总RNA,用1%的琼脂糖凝胶电泳初步检测总RNA的完整性。采用Nanodrop、Qubit 2.0和Aglient 2100方法检测各样品的浓度、纯度和完整性。选取RIN大于8的RNA样本进行cDNA文库制备,随后利用Illumina HiSeq 2000平台进行双末端测序,并利用FastQC(http：∥www.bioinformatics.babraham.ac.uk/projects/fastqc)对原始数据进行质控。之后使用TopHat进行基因组比对,edgeR进行基因差异表达分析。

1.2 基因组比对工具分析

首先在Ensembl数据库下载猪参考基因组文件及其注释文件(https：∥asia.ensembl.org/index.html)。通过TopHat2、HISAT2和STAR 3种工具将质控后数据比对到猪参考基因组Sus crofa11.1[17]上,对3种比对工具识别reads数、运行速度、比对率进行比较。并利用HT-seq软件计算每个样本的基因表达量[18],构建基因表达矩阵。以上步骤均在Linux系统中进行操作。

1.3 基因差异表达工具分析

通过DESeq2、edgeR和limma 3种工具筛选差异表达基因,以|log2Fold Change|>1且P<0.05、P<0.01、FDR<0.05为筛选条件判断基因表达量在两组之间是否存在差异。统计3种工具在不同条件下筛选到的差异基因数目及差异基因的交集。

1.4 KEGG通路富集分析

通过在线工具KOBAS[19]对筛选的差异表达基因进行KEGG功能富集分析。将3种筛选结果的上调表达和下调表达基因分别代入到KOBAS,选取物种为猪,之后对KEGG富集结果进行比较,统计3种上调表达和下调表达基因富集通路及其交集。

2 结果与分析

2.1 基因组比对分析

脂肪转录组数据比对结果表明,STAR和HISAT2识别reads数目相同,且与TopHat2识别数目无显著差异,3种工具识别reads数目平均为2.59×107以上(图2(a))。3种工具唯一比对率均在87%以上,STAR拥有最高的唯一比对率(图2(b))。此外,在相同线程数下,HISAT2处理一个样本用时平均为3 min,STAR用时平均为16 min,TopHat2用时平均为155 min。HISAT2平均运行速度比STAR快3～5倍,比TopHat2快约60～80倍。

(a)3种工具识别reads数;(b)3种工具比对率。(a) Number of reads identified by three tools;(b) Mapping rate of three tools.图2 3种工具识别reads数及比对率Fig.2 Reads number and mapping rate identified by three tools

2.2 差异表达基因分析

选取HISAT2-HTseq流程获得的基因表达数据进行后续基因差异表达分析。以|log2Fold Change|>1且P<0.05为标准对基因进行筛选。DESeq2筛选到616个差异基因,其中361个在松辽猪背部脂肪组织中上调表达,255个下调表达;edgeR筛选到890个差异基因,其中232个在松辽猪背部脂肪组织中上调表达,658个下调表达;limma筛选到829个差异基因,其中558个在松辽猪背部脂肪组织中上调表达,271个下调表达(图3(a))。3者有246个差异基因重合(图3(b))。同时,本研究还比较在P<0.01和FDR<0.05条件下,3种工具的差异基因。当筛选条件变严格时,差异基因数目相应会减少。而edgeR的结果所受影响较小,limma的结果受影响较大,尤其在FDR<0.05时,只筛选到15个差异基因(图3(c)和(d))。最终本研究选择P<0.05时的筛选结果进行后续分析。

(a)3种工具筛选到差异基因数目;(b)3种工具差异基因Venn图。筛选条件为|log2 Fold Change|>1且P<0.05;(c)3种工具差异基因Venn图,筛选条件为|log2 Fold Change|>1且P<0.01;(d)3种工具差异基因Venn图,筛选条件为|log2 Fold Change|>1且FDR<0.05。(a) The number of differential genes was screened by three tools;(b) Venn map of differential genes of three tools.The filter condition is |log2 Fold Change|>1 and P<0.05;(c) Venn map of differential genes of three tools.The filter condition is |log2 Fold Change|>1 and P<0.01;(d) Venn map of differential genes of three tools.The filter condition is |log2 Fold Change|>1 and FDR<0.05.图3 3种工具筛选的差异表达基因比较Fig.3 Comparison of differentially expressed genes screened by three tools

2.3 KEGG通路富集结果

通过在线工具KOBAS对3种差异表达工具筛选到的差异基因进行KEGG通路富集分析,并对排名前10条通路进行制图(图4(a)和(b)),DESeq2、edgeR和limma的上调差异表达基因分别富集到110、108和142条通路,其中有72条通路重合,而下调差异表达基因分别富集到190、247和177条通路,其中有158条通路重合(图4(c)和(d))。鉴于所用转录组数据目的为筛选影响脂肪沉积的基因,本研究中3种工具筛选出的差异表达基因富集结果中都存在与脂肪沉积相关的通路,例如脂肪酸代谢(Fatty acid metabolism)、不饱和脂肪酸生物合成(Biosynthesis of unsaturated fatty acids)、脂肪细胞因子信号通路(Adipocytokine signaling pathway)、脂肪细胞分解调节(Regulation of lipolysis in adipocytes)以及PPAR信号通路(PPAR signaling pathway)等。

(a)上调表达基因的KEGG通路图;(b)下调表达基因的KEGG通路图;(c)上调表达基因的KEGG通路Venn图;(d)下调表达基因的KEGG通路Venn图。(a) KEGG pathway map of up regulated genes;(b) KEGG pathway map of down regulated expression genes;(c) Venn map of KEGG pathway up regulated expression genes;(d) Venn map of KEGG pathway down regulated expression genes.图4 3种工具KEGG通路比较Fig.4 Comparison diagram of KEGG pathway of three tools

3 讨 论

本研究通过对RNA-seq技术不同流程分析表明,工具的选择对分析的准确性及运行时间有很大影响,HISAT2在3种比对工具中运行速度最快,这是由于HISAT2使用两种类型的索引进行比对：全局索引、全基因组索引和数万个小的本地索引[20]。两种索引都使用一种称为Burrows-Wheeler转换(BWT)的数据结构,该结构能够以高度压缩的形式存储参考基因组。并且使用一种称为Ferragina-Manzini (FM)索引[21]的特殊索引方案,这些使得HISAT2能够极其迅速地搜索一个基因组,从而获得以每小时数百万次读取为单位的比对速度。HISAT2的一些开发者同时是TopHat2的开发者,因而HISAT2作为TopHat2的继承者拥有更快的运行速度,并且所需要的内存更少,能够在传统台式机上运行。STAR拥有最高的唯一比对率,这与Sahraeian 等[5]研究结果一致,同时该研究中还表示STAR得出的数据整体质量较低。Costa等[22]研究表明不同比对工具对差异表达分析影响较小。通过综合考虑,本研究选取HISAT2数据进行后续差异表达基因筛选分析。

差异表达基因筛选结果表明DESeq2、edgeR和limma 3种工具筛选到差异表达基因有较大差异。当筛选条件为|log2Fold Change|>1且P<0.05时,DESeq2筛选到差异表达基因数目最少,edgeR筛选到差异基因数目最多,limma与edgeR数量相近,三者共同筛选到的基因有246个。这可能是由于3种工具对数据归一化的方法以及对低表达丰度基因进行过滤的方法不同造成的,这与Stupnikov等[23]的观点一致。EdgeR默认归一化处理方法为表达量对数值的加权平均值(Trimmed mean of M-values,TMM)[24],其假设对照组和处理组间绝大多数基因表达不发生差异,并比较每个样本的CPM(Counts per million)的上四分位数与全部样本的CPM的平均上四分位数之间的差值,选择差值最小的样本作为参考样本,之后以参考样本为基准进行校正。根据limma用户指南,本次研究中limma采用edgeR包的TMM标准化方法。而DESeq2归一化方法基于DESeq的相对对数表达式[3],具体为以所有样本的几何平均值为参考样本,计算每个样本与参考样本的比值,之后以给定样本的所有比例的中值作为该样本的标准化因子。DEseq2和edgeR都是基于负二项分布对原始数据进行处理,而limma通过voom函数将均值-方差关系转换为精度权重,从而提供了与基于负二项式的软件包同样的性能,并且对于大型数据具有更可靠的速度和可靠性[12,25-26]。EdgeR发现的差异基因数目最多,并且当筛选条件更加严格时,所受影响较小,但这可能会引入更多的假阳性。

通过KEGG通路富集分析发现3种差异表达基因筛选结果富集到通路条目相近,并且重合通路在各自富集结果中所占比例较大,经过对通路观察发现3组差异表达基因都能富集到与脂肪合成、分解、代谢相关的重要通路,富集到这些通路的基因很有可能是研究所需的目的基因。这表明虽然不同工具筛选到差异表达基因不同,但都能筛选到研究所需的目的基因。

每种差异分析方法都有其优点,可能适用于特定的RNA序列数据集。在RNA-Seq数据差异表达分析中,需要考虑生物学重复,测序深度等问题,当研究不存在生物学重复时,推荐使用edgeR,可以根据具体数据设置生物变异系数(BCV)值。为减少分析过程中的假阳性,可以选择2个或者多个工具的差异表达基因的交集。另外存在其他差异基因筛选工具,如NOIseq[27]、SAMseq[28]、Cuffdiff[29]等,进行多数据、多工具的组合比较可能得出更合适的工作流程,得到更好地研究结果。

4 结 论

不同RNA-seq步骤分析表明,分析工具的选择对结果的准确性及运行时间有很大影响,比对工具中HISAT2运行速度最快,因此本研究推荐使用HISAT2进行基因组比对。当研究不存在生物学重复时,推荐使用edgeR进行差异表达基因筛选。而为了减少分析过程中的假阳性,可以选择DESeq2或者2个及以上工具的差异表达基因的交集。本研究方法将有助于研究人员从转录组数据中获得更好、更全面的生物学见解。