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

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

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

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

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

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

1 材料與方法

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

本研究使用轉(zhuǎn)錄組數(shù)據(jù)來(lái)自參考文獻(xiàn)[16],項(xiàng)目號(hào)為PRJNA234335和PRJNA287471。具體為選取松遼黑豬和長(zhǎng)白豬各6頭,食物充足、飲水自由、飼養(yǎng)條件相同。在屠宰場(chǎng)宰殺后取背部脂肪組織,并使用Trizol試劑盒提取12份背部脂肪總RNA,用1%的瓊脂糖凝膠電泳初步檢測(cè)總RNA的完整性。采用Nanodrop、Qubit 2.0和Aglient 2100方法檢測(cè)各樣品的濃度、純度和完整性。選取RIN大于8的RNA樣本進(jìn)行cDNA文庫(kù)制備,隨后利用Illumina HiSeq 2000平臺(tái)進(jìn)行雙末端測(cè)序,并利用FastQC(http：∥www.bioinformatics.babraham.ac.uk/projects/fastqc)對(duì)原始數(shù)據(jù)進(jìn)行質(zhì)控。之后使用TopHat進(jìn)行基因組比對(duì),edgeR進(jìn)行基因差異表達(dá)分析。

1.2 基因組比對(duì)工具分析

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

1.3 基因差異表達(dá)工具分析

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

1.4 KEGG通路富集分析

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

2 結(jié)果與分析

2.1 基因組比對(duì)分析

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

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

2.2 差異表達(dá)基因分析

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

2.3 KEGG通路富集結(jié)果

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

(a)上調(diào)表達(dá)基因的KEGG通路圖;(b)下調(diào)表達(dá)基因的KEGG通路圖;(c)上調(diào)表達(dá)基因的KEGG通路Venn圖;(d)下調(diào)表達(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 討 論

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

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

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

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

4 結(jié) 論

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