RNA-seq數(shù)據(jù)差異表達分析流程比較

趙延輝 陳少康 翟麗維 侍玉梅 原佳妮 盛熙暉 齊曉龍 郭 勇 王楚端 邢 凱*

(1.北京農(nóng)學(xué)院 動物科學(xué)技術(shù)學(xué)院,北京 102206;2.北京市畜牧總站,北京 100107;3.中國農(nóng)業(yè)大學(xué) 動物科學(xué)技術(shù)學(xué)院,北京 100193)

RNA-seq即轉(zhuǎn)錄組測序技術(shù),能夠全面快速地獲得某一物種特定組織或器官在某一狀態(tài)下的幾乎所有轉(zhuǎn)錄本序列信息。隨著高通量測序技術(shù)改進、成本下降,RNA-seq已經(jīng)廣泛應(yīng)用于多物種的研究中。與DNA微陣列技術(shù)相比,RNA-seq數(shù)據(jù)可重復(fù)性強,背景噪聲低,檢測范圍廣泛[1],并且可以在沒有參考基因組的情況下使用,包括從頭測序[2]、豐度估計[3]、檢測選擇性剪切[4]等。RNA-seq所選工作流程不同,分析的準(zhǔn)確性、速度、成本可能會大不相同[5]。

RNA-seq技術(shù)最常用的用途就是尋找特定條件下的差異表達基因,其工作流程首先將樣本RNA片段化并且反轉(zhuǎn)錄成cDNA進行測序,將測序獲得的短序列比對到參考基因組上,通過比對到基因組的reads數(shù)目來估算基因的相對表達水平,之后通過統(tǒng)計學(xué)方法檢驗組間基因差異表達[6]。目前,針對此流程已有許多工具被開發(fā)出來,如比對工具中TopHat2[7]、HISAT2[8]和STAR[9]被廣泛使用,TopHat2可以將測序產(chǎn)生的各種長度的reads,即使是在高度重復(fù)的基因組或存在假基因的情況下,也能產(chǎn)生敏感和準(zhǔn)確的比對。HISAT2是目前可用的最快系統(tǒng),與其他工具相比其精度更好或相同,另外HISAT2需要計算機運行內(nèi)存少,并且支持任何大小的基因組。STAR同樣擁有良好的對齊精確度和靈敏度,同時STAR的運行速度也是其優(yōu)勢之一。DESeq2、edgeR和limma是3種常用的基因差異表達分析工具,DESeq2是Love等[10]在DESeq的基礎(chǔ)上開發(fā)出的Bioconductor包,其對基因排序和可視化進行了改進,擁有較高的靈敏性和準(zhǔn)確性,并且控制了假陽性率。EdgeR是由Mark等[11]開發(fā)的Bioconductor包,可以在沒有生物學(xué)重復(fù)的情況下尋找差異表達基因。Limma最開始時只能處理來自微陣列的基因表達數(shù)據(jù)[12],自3.9.19版本開始可以用于分析RNA-seq數(shù)據(jù)。

少數(shù)研究比較了不同差異表達分析工具的性能[13-15],但是并未考慮基因組比對工具可能帶來的影響。此外對于哪種分析工具是最佳的,以及如何確保分析結(jié)果的準(zhǔn)確性和復(fù)現(xiàn)性,目前尚未達成共識。不同的工具進行組合產(chǎn)生不同的工作流程,不同的工作流程產(chǎn)生的結(jié)果可能大不相同。本研究比較了TopHat2、HISAT2和STAR 3種廣泛使用的比對工具以及DESeq2、edgeR和limma 3種常用的基因差異表達分析工具,并綜合KEGG通路分析結(jié)果對不同分析流程進行評價(圖1),選擇出合適的分析流程為后續(xù)研究工作提供參考。

圖1 RNA-Seq數(shù)據(jù)差異表達分析工作流程Fig.1 The workflow of differential expression analysis for RNA-Seq data

1 材料與方法

1.1 轉(zhuǎn)錄組數(shù)據(jù)

本研究使用轉(zhuǎn)錄組數(shù)據(jù)來自參考文獻[16],項目號為PRJNA234335和PRJNA287471。具體為選取松遼黑豬和長白豬各6頭,食物充足、飲水自由、飼養(yǎng)條件相同。在屠宰場宰殺后取背部脂肪組織,并使用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)對原始數(shù)據(jù)進行質(zhì)控。之后使用TopHat進行基因組比對,edgeR進行基因差異表達分析。

1.2 基因組比對工具分析

首先在Ensembl數(shù)據(jù)庫下載豬參考基因組文件及其注釋文件(https：∥asia.ensembl.org/index.html)。通過TopHat2、HISAT2和STAR 3種工具將質(zhì)控后數(shù)據(jù)比對到豬參考基因組Sus crofa11.1[17]上,對3種比對工具識別reads數(shù)、運行速度、比對率進行比較。并利用HT-seq軟件計算每個樣本的基因表達量[18],構(gòu)建基因表達矩陣。以上步驟均在Linux系統(tǒng)中進行操作。

1.3 基因差異表達工具分析

通過DESeq2、edgeR和limma 3種工具篩選差異表達基因,以|log2Fold Change|>1且P<0.05、P<0.01、FDR<0.05為篩選條件判斷基因表達量在兩組之間是否存在差異。統(tǒng)計3種工具在不同條件下篩選到的差異基因數(shù)目及差異基因的交集。

1.4 KEGG通路富集分析

通過在線工具KOBAS[19]對篩選的差異表達基因進行KEGG功能富集分析。將3種篩選結(jié)果的上調(diào)表達和下調(diào)表達基因分別代入到KOBAS,選取物種為豬,之后對KEGG富集結(jié)果進行比較,統(tǒng)計3種上調(diào)表達和下調(diào)表達基因富集通路及其交集。

2 結(jié)果與分析

2.1 基因組比對分析

脂肪轉(zhuǎn)錄組數(shù)據(jù)比對結(jié)果表明,STAR和HISAT2識別reads數(shù)目相同,且與TopHat2識別數(shù)目無顯著差異,3種工具識別reads數(shù)目平均為2.59×107以上(圖2(a))。3種工具唯一比對率均在87%以上,STAR擁有最高的唯一比對率(圖2(b))。此外,在相同線程數(shù)下,HISAT2處理一個樣本用時平均為3 min,STAR用時平均為16 min,TopHat2用時平均為155 min。HISAT2平均運行速度比STAR快3～5倍,比TopHat2快約60～80倍。

(a)3種工具識別reads數(shù);(b)3種工具比對率。(a) Number of reads identified by three tools;(b) Mapping rate of three tools.圖2 3種工具識別reads數(shù)及比對率Fig.2 Reads number and mapping rate identified by three tools

2.2 差異表達基因分析

選取HISAT2-HTseq流程獲得的基因表達數(shù)據(jù)進行后續(xù)基因差異表達分析。以|log2Fold Change|>1且P<0.05為標(biāo)準(zhǔn)對基因進行篩選。DESeq2篩選到616個差異基因,其中361個在松遼豬背部脂肪組織中上調(diào)表達,255個下調(diào)表達;edgeR篩選到890個差異基因,其中232個在松遼豬背部脂肪組織中上調(diào)表達,658個下調(diào)表達;limma篩選到829個差異基因,其中558個在松遼豬背部脂肪組織中上調(diào)表達,271個下調(diào)表達(圖3(a))。3者有246個差異基因重合(圖3(b))。同時,本研究還比較在P<0.01和FDR<0.05條件下,3種工具的差異基因。當(dāng)篩選條件變嚴(yán)格時,差異基因數(shù)目相應(yīng)會減少。而edgeR的結(jié)果所受影響較小,limma的結(jié)果受影響較大,尤其在FDR<0.05時,只篩選到15個差異基因(圖3(c)和(d))。最終本研究選擇P<0.05時的篩選結(jié)果進行后續(xù)分析。

(a)3種工具篩選到差異基因數(shù)目;(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通路富集結(jié)果

通過在線工具KOBAS對3種差異表達工具篩選到的差異基因進行KEGG通路富集分析,并對排名前10條通路進行制圖(圖4(a)和(b)),DESeq2、edgeR和limma的上調(diào)差異表達基因分別富集到110、108和142條通路,其中有72條通路重合,而下調(diào)差異表達基因分別富集到190、247和177條通路,其中有158條通路重合(圖4(c)和(d))。鑒于所用轉(zhuǎn)錄組數(shù)據(jù)目的為篩選影響脂肪沉積的基因,本研究中3種工具篩選出的差異表達基因富集結(jié)果中都存在與脂肪沉積相關(guān)的通路,例如脂肪酸代謝(Fatty acid metabolism)、不飽和脂肪酸生物合成(Biosynthesis of unsaturated fatty acids)、脂肪細胞因子信號通路(Adipocytokine signaling pathway)、脂肪細胞分解調(diào)節(jié)(Regulation of lipolysis in adipocytes)以及PPAR信號通路(PPAR signaling pathway)等。

(a)上調(diào)表達基因的KEGG通路圖;(b)下調(diào)表達基因的KEGG通路圖;(c)上調(diào)表達基因的KEGG通路Venn圖;(d)下調(diào)表達基因的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技術(shù)不同流程分析表明,工具的選擇對分析的準(zhǔn)確性及運行時間有很大影響,HISAT2在3種比對工具中運行速度最快,這是由于HISAT2使用兩種類型的索引進行比對：全局索引、全基因組索引和數(shù)萬個小的本地索引[20]。兩種索引都使用一種稱為Burrows-Wheeler轉(zhuǎn)換(BWT)的數(shù)據(jù)結(jié)構(gòu),該結(jié)構(gòu)能夠以高度壓縮的形式存儲參考基因組。并且使用一種稱為Ferragina-Manzini (FM)索引[21]的特殊索引方案,這些使得HISAT2能夠極其迅速地搜索一個基因組,從而獲得以每小時數(shù)百萬次讀取為單位的比對速度。HISAT2的一些開發(fā)者同時是TopHat2的開發(fā)者,因而HISAT2作為TopHat2的繼承者擁有更快的運行速度,并且所需要的內(nèi)存更少,能夠在傳統(tǒng)臺式機上運行。STAR擁有最高的唯一比對率,這與Sahraeian 等[5]研究結(jié)果一致,同時該研究中還表示STAR得出的數(shù)據(jù)整體質(zhì)量較低。Costa等[22]研究表明不同比對工具對差異表達分析影響較小。通過綜合考慮,本研究選取HISAT2數(shù)據(jù)進行后續(xù)差異表達基因篩選分析。

差異表達基因篩選結(jié)果表明DESeq2、edgeR和limma 3種工具篩選到差異表達基因有較大差異。當(dāng)篩選條件為|log2Fold Change|>1且P<0.05時,DESeq2篩選到差異表達基因數(shù)目最少,edgeR篩選到差異基因數(shù)目最多,limma與edgeR數(shù)量相近,三者共同篩選到的基因有246個。這可能是由于3種工具對數(shù)據(jù)歸一化的方法以及對低表達豐度基因進行過濾的方法不同造成的,這與Stupnikov等[23]的觀點一致。EdgeR默認(rèn)歸一化處理方法為表達量對數(shù)值的加權(quán)平均值(Trimmed mean of M-values,TMM)[24],其假設(shè)對照組和處理組間絕大多數(shù)基因表達不發(fā)生差異,并比較每個樣本的CPM(Counts per million)的上四分位數(shù)與全部樣本的CPM的平均上四分位數(shù)之間的差值,選擇差值最小的樣本作為參考樣本,之后以參考樣本為基準(zhǔn)進行校正。根據(jù)limma用戶指南,本次研究中l(wèi)imma采用edgeR包的TMM標(biāo)準(zhǔn)化方法。而DESeq2歸一化方法基于DESeq的相對對數(shù)表達式[3],具體為以所有樣本的幾何平均值為參考樣本,計算每個樣本與參考樣本的比值,之后以給定樣本的所有比例的中值作為該樣本的標(biāo)準(zhǔn)化因子。DEseq2和edgeR都是基于負(fù)二項分布對原始數(shù)據(jù)進行處理,而limma通過voom函數(shù)將均值-方差關(guān)系轉(zhuǎn)換為精度權(quán)重,從而提供了與基于負(fù)二項式的軟件包同樣的性能,并且對于大型數(shù)據(jù)具有更可靠的速度和可靠性[12,25-26]。EdgeR發(fā)現(xiàn)的差異基因數(shù)目最多,并且當(dāng)篩選條件更加嚴(yán)格時,所受影響較小,但這可能會引入更多的假陽性。

通過KEGG通路富集分析發(fā)現(xiàn)3種差異表達基因篩選結(jié)果富集到通路條目相近,并且重合通路在各自富集結(jié)果中所占比例較大,經(jīng)過對通路觀察發(fā)現(xiàn)3組差異表達基因都能富集到與脂肪合成、分解、代謝相關(guān)的重要通路,富集到這些通路的基因很有可能是研究所需的目的基因。這表明雖然不同工具篩選到差異表達基因不同,但都能篩選到研究所需的目的基因。

每種差異分析方法都有其優(yōu)點,可能適用于特定的RNA序列數(shù)據(jù)集。在RNA-Seq數(shù)據(jù)差異表達分析中,需要考慮生物學(xué)重復(fù),測序深度等問題,當(dāng)研究不存在生物學(xué)重復(fù)時,推薦使用edgeR,可以根據(jù)具體數(shù)據(jù)設(shè)置生物變異系數(shù)(BCV)值。為減少分析過程中的假陽性,可以選擇2個或者多個工具的差異表達基因的交集。另外存在其他差異基因篩選工具,如NOIseq[27]、SAMseq[28]、Cuffdiff[29]等,進行多數(shù)據(jù)、多工具的組合比較可能得出更合適的工作流程,得到更好地研究結(jié)果。

4 結(jié) 論

不同RNA-seq步驟分析表明,分析工具的選擇對結(jié)果的準(zhǔn)確性及運行時間有很大影響,比對工具中HISAT2運行速度最快,因此本研究推薦使用HISAT2進行基因組比對。當(dāng)研究不存在生物學(xué)重復(fù)時,推薦使用edgeR進行差異表達基因篩選。而為了減少分析過程中的假陽性,可以選擇DESeq2或者2個及以上工具的差異表達基因的交集。本研究方法將有助于研究人員從轉(zhuǎn)錄組數(shù)據(jù)中獲得更好、更全面的生物學(xué)見解。