RobustDEA：一種快速魯棒的RNA-Seq數(shù)據(jù)尋找差異表達(dá)基因方法

張 禮， 王嘉瑞， 吳東洋

(南京林業(yè)大學(xué) 信息科學(xué)與技術(shù)學(xué)院，南京 210016)

下一代高通量DNA測(cè)序技術(shù)已經(jīng)得到快速發(fā)展和普及，轉(zhuǎn)錄組測(cè)序技術(shù)(RNA-Sequencing，RNA-Seq)具有高通量，高信噪比，可重復(fù)性好，所需樣本少等特點(diǎn)，被廣泛應(yīng)用的轉(zhuǎn)錄組學(xué)的各個(gè)研究領(lǐng)域中[1].尋找差異基因表達(dá)分析作為RNA-Seq數(shù)據(jù)分析中最基本應(yīng)用之一，其目的是在不同對(duì)照組和控制組樣本中，比如正常狀態(tài)和藥物治療，健康和疾病個(gè)體，不同生物組織或者不同發(fā)育階段等，識(shí)別出具有差異表達(dá)的基因.從RNA-Seq測(cè)序?qū)嶒?yàn)中獲得海量讀段數(shù)據(jù)后，首先需要將讀段數(shù)據(jù)定位到參考基因組或者轉(zhuǎn)錄組上，從而統(tǒng)計(jì)出每個(gè)基因或者外顯子上的讀段數(shù)目.一旦獲得基因或者外顯子的讀段數(shù)目，就可以進(jìn)行差異表達(dá)分析.盡管RNA-Seq技術(shù)相比傳統(tǒng)基因芯片具有很多優(yōu)勢(shì)，但是由于RNA-Seq測(cè)序本身技術(shù)上限制產(chǎn)生了很多困難，比如多源讀段映射，測(cè)序偏好和小樣本等問題，導(dǎo)致差異表達(dá)分析仍然是一個(gè)極具挑戰(zhàn)的任務(wù).

針對(duì)尋找差異表達(dá)基因所面臨的種種挑戰(zhàn)，提出了大量統(tǒng)計(jì)方法，通?？梢苑譃閮深悾夯谧x段數(shù)目的方法和基于表達(dá)水平的兩步法.基于讀段數(shù)目的方法是比較在不同樣本條件中基因上映射的讀段數(shù)目.由于在不同樣本之間，基因上讀段數(shù)目變異性較大，若采用泊松分布會(huì)導(dǎo)致讀段分布存在過散步(overdispersion)問題.因此，基于讀段數(shù)目方法通常采用負(fù)二項(xiàng)(Negative binomial)分布，比如DESeq[2]，DESeq2[3]，baySeq[4]，edgeR[5]等.相比泊松分布的單參數(shù)，負(fù)二項(xiàng)分布增加一個(gè)參數(shù)，能單獨(dú)描述數(shù)據(jù)方差特性，從而能解決讀段數(shù)據(jù)的過散步問題.此外，Voom使用正太線性模型來擬合讀段數(shù)據(jù)，并且嵌入了數(shù)據(jù)分布的均值和方差的相關(guān)性信息[6].SAMSeq[7]，NOISeq[8]以及NPEBSeq[9]等非參數(shù)模型，無需假設(shè)基因讀段分布服從任何概率分布.但這類方法都忽略了讀段的多源映射和數(shù)據(jù)偏差等問題，因此提出了基于表達(dá)水平的兩步法：第一步利用現(xiàn)有表達(dá)水平估計(jì)方法計(jì)算出基因的表達(dá)水平，比如Cufflinks[10]，PGSeq[11]等，這些方法在計(jì)算過程中考慮了讀段多源映射和數(shù)據(jù)偏差等問題，從而獲得更為準(zhǔn)確的基因表達(dá)水平；第二步是利用獲得的基因表達(dá)水平進(jìn)行差異分析，比如CuffDiff2與Cufflinks，BDSeq[12]與PGSeq.由于考慮了更多的數(shù)據(jù)特征，基于表達(dá)水平的兩步法通常能獲得較好的結(jié)果.但由于兩步法需要預(yù)先計(jì)算基因表達(dá)水平，其時(shí)間效率要比基于讀段數(shù)目的方法低很多.因此，在面向?qū)ふ也町惐磉_(dá)基因的單一任務(wù)時(shí)，目前主流做法是選擇基于讀段數(shù)目的方法，在較短的時(shí)間內(nèi)獲得較高的準(zhǔn)確度[13].

針對(duì)RNA-Seq數(shù)據(jù)分析發(fā)現(xiàn)，同一個(gè)數(shù)據(jù)集中的基因，采用不同的尋找差異表達(dá)基因方法中會(huì)得到差異較大的結(jié)果[14].采用DESeq，baySeq，edgeR和Voom 4個(gè)方法來尋找MAQC數(shù)據(jù)集中的差異表達(dá)基因，表1展示了兩個(gè)基因在不同方法中的P值和差異重要性排序，說明不同方法在分析同一個(gè)數(shù)據(jù)集是存在不一致性.導(dǎo)致這種現(xiàn)象除了數(shù)據(jù)集本身的生物特異性差異，以及不同研究者發(fā)表的基因表達(dá)數(shù)據(jù)的質(zhì)量差異外，不同模型的假設(shè)不同或特點(diǎn)數(shù)據(jù)集偏好性，如NOISeq比較適合小樣本數(shù)據(jù)集[15].此外，不同方法的結(jié)果報(bào)告格式不一致通常會(huì)導(dǎo)致結(jié)果理解和推斷的困難.

表1 不同方法下基因載MAQC數(shù)據(jù)集上的結(jié)果

為了克服上述的問題，一種潛在可行方法是將當(dāng)前算法結(jié)合起來，以獲得更魯棒的差異基因列表，其特點(diǎn)在于具有更高的統(tǒng)計(jì)能力，更少的假陽性和假陰性差異基因.研究人員通常是結(jié)合不同方法在同一個(gè)數(shù)據(jù)集上分析結(jié)果，從而推斷出更穩(wěn)定的結(jié)論.IDEAMEX是一個(gè)在線差異表達(dá)基因檢測(cè)軟件[16]，該軟件首先預(yù)設(shè)閾值得到每個(gè)方法所識(shí)別出的差異表達(dá)基因，在采用簡(jiǎn)單的投票方式選出4個(gè)方法都公認(rèn)差異表達(dá)基因列表.consensusDE[17]和MultiRankSeq[13]方法都采用了基因排名求和的思想.首先不同方法根據(jù)p-value值大小對(duì)基因分別進(jìn)行排序，然后將基因在不同方法中的排名求和，從而得到基因最終排名序列，排名越靠前的基因其越差異表達(dá)的可能性就越大.不同之處在于，MultiRankSeq方法結(jié)合了DESeq，edgeR和baySeq 3個(gè)方法，數(shù)據(jù)驗(yàn)證顯示，DESeq和edgeR的差異基因集很接近，而baySeq找到的差異基因卻和上述兩個(gè)方法重疊的部分較少.圖1是4個(gè)方法前1 000個(gè)最有可能的差異基因的Venn圖，可以發(fā)現(xiàn)baySeq方法與其他3個(gè)方法的重疊基因很少，這將導(dǎo)致MultiRankSeq的結(jié)果易受到baySeq的影響.因此consensusDE方法采用了DESeq2，edgeR和Voom 3個(gè)方法，并在數(shù)據(jù)預(yù)處理步驟考慮消除不變基因的影響.上述方法通過投票或者排名求和的結(jié)合方式，找到被所結(jié)合方法共同認(rèn)為最有可能的差異基因，但這些方法容易受到結(jié)合方法中單個(gè)方法的影響，反而降低了方法的準(zhǔn)確性和魯棒性.為了減輕單個(gè)方法對(duì)融合結(jié)果的影響，PANDORA方法對(duì)DESeq，edgeR，Vomm，NBPSeq，NOISeq和baySeq 6個(gè)方法進(jìn)行加權(quán)融合[18].但是PANDORA方法需要預(yù)先產(chǎn)生與測(cè)試數(shù)據(jù)集特征相近的模擬數(shù)據(jù)集，通過模擬數(shù)據(jù)集計(jì)算出每個(gè)方法的權(quán)重值，此過程極其耗時(shí).此外，PANDORA和MltiRankSeq都采用了DESeq，而不是性能更優(yōu)越的DESeq2.

圖1 DESeq, edgeR, Voom和baySeq

針對(duì)上述方法存在的問題，文中提出了一種魯棒的尋找差異表達(dá)基因表達(dá)方法(Robust differential expression analysis based on RNA-Seq data，RobustDEA).RobustDEA方法選擇了不同數(shù)據(jù)集中表現(xiàn)較好且穩(wěn)定的4個(gè)方法：DESeq2， Voom， edgeR和NOISeq，通過加權(quán)方式結(jié)合多個(gè)尋找差異表達(dá)基因方法，其權(quán)值可快速?gòu)臄?shù)據(jù)集中學(xué)習(xí)獲到，并能有效的反映數(shù)據(jù)集的特點(diǎn)，從而使得RobustDEA方法可以在不同數(shù)據(jù)集上都能獲得更穩(wěn)定的結(jié)果.通過多個(gè)真實(shí)RNA-Seq測(cè)序數(shù)據(jù)集的驗(yàn)證表明，RobustDEA方法能在不同數(shù)據(jù)集上獲得更加魯棒的差異基因集.

1 方法

1.1 常見結(jié)合方法

基因芯片數(shù)據(jù)和RNA-Seq數(shù)據(jù)的差異表達(dá)基因分析中，多種結(jié)合方法被廣泛應(yīng)用到對(duì)不同方法檢測(cè)結(jié)果的融合.目前，應(yīng)用最為廣泛的兩種常規(guī)統(tǒng)計(jì)方法和專門針對(duì)RNA-Seq數(shù)據(jù)的最流行 PANDORA方法.首先，假設(shè)基因i使用尋找差異表達(dá)基因方法j獲得其對(duì)應(yīng)的P值pij，然后可以使用下面幾種結(jié)合方法來進(jìn)行不同方法結(jié)果的融合.

1.1.1 Fisher方法

根據(jù)Fisher方法，基因i通過m個(gè)方法估計(jì)的P值合并為：

(1)

合并后得到的統(tǒng)計(jì)量符合自由度為2m的卡方分布.χ2值越大，P值就越小，表示對(duì)應(yīng)的基因其差異性就越明顯.但是Fisher方法設(shè)計(jì)目的是合并具有相同零假設(shè)的統(tǒng)計(jì)檢驗(yàn)方法在不同數(shù)據(jù)集中產(chǎn)生的P值，并不是合并同一個(gè)數(shù)據(jù)集上由不同統(tǒng)計(jì)檢驗(yàn)方法產(chǎn)生的P值.

1.2.2 Stouffer方法

針對(duì)每個(gè)基因，標(biāo)準(zhǔn)的Stouffer方法首先將pj值轉(zhuǎn)換為Zj=Φ-1(1-pj)，其中Φ是標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)，每個(gè)基因的Z值為：

(2)

根據(jù)Z值大小對(duì)基因差異表達(dá)排序.標(biāo)準(zhǔn)Stouffer方法的優(yōu)點(diǎn)易于為不同方法賦予不同權(quán)重值，因此提出改進(jìn)的加權(quán)Stouffer方法：

(3)

式中：ωj為方法j的權(quán)重，其值可根據(jù)具體問題來選擇，比如數(shù)據(jù)集的樣本規(guī)模等.Z值越大，對(duì)應(yīng)P值就越小，表示該基因的差異性就越顯著.Stouffer方法和Fisher方法密切相關(guān)，其設(shè)計(jì)目的的假設(shè)都是合并具有相同零假設(shè)的統(tǒng)計(jì)檢驗(yàn)方法在不同數(shù)據(jù)集中產(chǎn)生的P值.但加權(quán)Stouffer方法能便于對(duì)不同方法進(jìn)行加權(quán)，而Fisher方法的卡方統(tǒng)計(jì)量，存在加權(quán)困難問題.

1.2.3 PANDORA方法

PANDORA方法的模型為：

(4)

式中：ωj為方法j的權(quán)重，其值自動(dòng)從模擬數(shù)據(jù)集中獲得[18].此模擬數(shù)據(jù)集是基于正研究的數(shù)據(jù)集而產(chǎn)生，因此其權(quán)值能更為貼近研究的數(shù)據(jù)集.但模擬數(shù)據(jù)集的產(chǎn)生極其耗時(shí)，與結(jié)合方法的個(gè)數(shù)，測(cè)序樣本數(shù)，基因個(gè)數(shù)等諸多因素有直接關(guān)系.

1.2 RobustDEA方法

針對(duì)不同尋找差異表達(dá)基因方法檢測(cè)得到的差異基因集存在不一致情況，基于統(tǒng)計(jì)學(xué)中元分析的思想，文中提出一個(gè)快速魯棒的RNA-Seq數(shù)據(jù)尋找差異表達(dá)基因方法RobustDEA，結(jié)合多個(gè)不同的差異表達(dá)基因方法，其模型為：

(5)

式中：ωj為方法j的權(quán)重值.RobustDEA方法采用自動(dòng)加權(quán)方式，其權(quán)重值從數(shù)據(jù)集中自動(dòng)計(jì)算得到.結(jié)合的每個(gè)方法先分別對(duì)數(shù)據(jù)集進(jìn)行分析，獲得每個(gè)基因的P值，可計(jì)算出每個(gè)方法在此數(shù)據(jù)集上的接收者操作特征曲線(receiver operating characteristic curve,ROC)下的面積(area under curve,AUC)值，從而自動(dòng)估計(jì)出每個(gè)方法所對(duì)應(yīng)的權(quán)重值為：

(6)

式中：AUCj為方法j的AUC值.通過自動(dòng)加權(quán)方式，可以獲得與當(dāng)前分析數(shù)據(jù)集最合適的權(quán)重值，有利于提高RobustDEA方法的精確度和魯棒性.

在RNA-Seq數(shù)據(jù)集，計(jì)算各個(gè)方法的AUC值需要大量基準(zhǔn)基因(即已知差異表達(dá)的基因)，基準(zhǔn)基因的差異表達(dá)值需要通過實(shí)時(shí)熒光是量聚合酶鏈反應(yīng)技術(shù)(quantitative real-time PCR,qRT-PCR)生物實(shí)驗(yàn)來確定.而現(xiàn)實(shí)中包含PCR數(shù)據(jù)的RNA-Seq數(shù)據(jù)集是極其少的.因此，RobustDEA方法采用了一種相互交叉差異基因集的驗(yàn)證方式來計(jì)算各個(gè)方法的AUC值，從而避開對(duì)PCR數(shù)據(jù)的依賴.針對(duì)每個(gè)方法，先計(jì)算并排序每個(gè)基因的P值，選取前N個(gè)基因(N通常選擇2 000)，這部分基因被認(rèn)為最有可能的差異表達(dá)基因.再計(jì)算多個(gè)方法的差異基因的交集.交集中基因是被共同認(rèn)可的差異基因，以此作為該數(shù)據(jù)集的基準(zhǔn)基因，用來計(jì)算自身的AUC值，其算法的具體過程如算法1.

算法1:相互交叉差異基因集的驗(yàn)證AUC值計(jì)算方法輸入:m個(gè)方法對(duì)差異基因列表(根據(jù)P值排序)Repeat:對(duì)于方法j 1. 分別選取其它方法的前N個(gè)差異表達(dá)基因; 2. 前N個(gè)差異表達(dá)基因的交集; 3. 以交集中基因作為基準(zhǔn)基因; 4. 計(jì)算AUCj值Until:所有方法的AUC值計(jì)算完畢輸出:根據(jù)各個(gè)方法的AUC值計(jì)算權(quán)重大小

通過AUC值計(jì)算的相互交叉驗(yàn)證，可避免PANDORA方法中需要預(yù)先使用模擬數(shù)據(jù)集去計(jì)算各個(gè)方法權(quán)重的過程，不僅可大幅度的提高權(quán)重的計(jì)算效率，同時(shí)也能保留數(shù)據(jù)集的特點(diǎn).

1.3 評(píng)估準(zhǔn)則

選取了多種評(píng)估準(zhǔn)則來評(píng)估不同方法和驗(yàn)證RobustDEA方法的性能.當(dāng)RNA-Seq數(shù)據(jù)集包含PCR驗(yàn)證基因，或者將結(jié)合方法共同認(rèn)可的差異基因作為基準(zhǔn)基因，選擇F1值來評(píng)估方法性能.F1值是對(duì)精確率(Precision)和召回率(Recall)的加權(quán)調(diào)和平均，其定義為：

此外，進(jìn)一步采用ROC曲線和AUC值評(píng)估按方法的預(yù)測(cè)準(zhǔn)確性.ROC曲線圖是反映橫坐標(biāo)假陽性率(false positive rate, FPR)與縱坐標(biāo)真陽性率(true positive rate, TPR)之間關(guān)系的曲線，其中：

AUC值表示ROC曲線下面積，其值越大表示正樣本越有可能被排在負(fù)樣本之前，即方法的預(yù)測(cè)性越好.

1.4 RobustDEA實(shí)現(xiàn)

為了便于用戶使用，RobustDEA方法采用R語言實(shí)現(xiàn)，且選擇的結(jié)合方法都在Bioconductor中提供了R語言包，避免了用戶額外安裝其它基于Unix的軟件，比如Cuffdiff2等.

RobustDEA方法的R語言包和詳細(xì)文檔可以從網(wǎng)址https://github.com/BIO-LAB/RobustDEA下載，供全球生物學(xué)家免費(fèi)使用.

2 實(shí)驗(yàn)結(jié)果與討論

2.1 數(shù)據(jù)集

為了評(píng)估RobustDEA方法，文中選擇了ReCount數(shù)據(jù)庫中的3個(gè)真實(shí)RNA-Seq數(shù)據(jù)集[19]，分別是人類大腦數(shù)據(jù)集MAQC，老鼠數(shù)據(jù)集Katz和小鼠數(shù)據(jù)集Hammer.

人類大腦MAQC數(shù)據(jù)集來自美國(guó)藥品監(jiān)督局(FDA)聯(lián)合全球多所研究機(jī)構(gòu)進(jìn)行的“生物芯片質(zhì)量控制”項(xiàng)目[20].MAQC數(shù)據(jù)集包含2個(gè)條件，正常大腦組織(HBR)和病變大腦組織(UHR)，每個(gè)條件下包含7個(gè)重復(fù)試驗(yàn)樣本，總計(jì)包含52 580個(gè)基因.過濾掉在所有樣本中讀段數(shù)目都是零的基因，獲得11 907個(gè)具有表達(dá)的基因.此外，MAQC項(xiàng)目提供了1 000個(gè)qRT-PCR驗(yàn)證基因.通過篩選具有顯著變化的基因(FC<=0.5或FC>=2.0)，獲得305個(gè)驗(yàn)證基因作為標(biāo)準(zhǔn)基因，用來評(píng)估不同方法尋找差異表達(dá)基因的準(zhǔn)確性.

大鼠數(shù)據(jù)集Katz包含兩個(gè)實(shí)驗(yàn)條件，正常成肌細(xì)胞和敲掉CUGBP1的成肌細(xì)胞，每個(gè)條件包含2個(gè)生物性重復(fù)實(shí)驗(yàn)，總計(jì)包含36 536個(gè)基因.過濾掉在所有樣本中讀段數(shù)目都是零的基因，獲得9 501個(gè)具有表達(dá)的基因.

小鼠數(shù)據(jù)集Hammer是關(guān)于小鼠脊柱神經(jīng)發(fā)育的研究，包含胚胎發(fā)育后2周和2月兩個(gè)時(shí)間點(diǎn).每個(gè)時(shí)間點(diǎn)包含2個(gè)條件，正常脊柱神經(jīng)細(xì)胞和進(jìn)行脊柱神經(jīng)結(jié)扎的神經(jīng)細(xì)胞.每個(gè)條件包含2個(gè)生物性重復(fù)樣本，總計(jì)包含29 516個(gè)基因，過濾掉在不同時(shí)間點(diǎn)時(shí)樣本中讀段數(shù)目都是零的基因，兩個(gè)時(shí)間點(diǎn)分別獲得17 125和18 463個(gè)具有表達(dá)的基因.

2.2 結(jié)合方法的選擇

目前，提出了大量尋找差異基因表達(dá)方法，結(jié)合方法的選擇會(huì)對(duì)RobustDEA方法的穩(wěn)定性和性能影響較大.因此，在不同數(shù)據(jù)集中性能穩(wěn)定且時(shí)間效率高的方法作為RobustDEA的結(jié)合方法的最有選擇.根據(jù)多個(gè)綜述性評(píng)估不同方法的精確度和時(shí)間效率(如表2)，文中選擇DESeq2(v1.26.0) ，edgeR(v3.28.1)，NOISeq(v2.30.0) 和Voom(v3.42.2)4種方法作為RobustDEA的原始方法.PANDORA方法來自R包metaseqR(v1.26.0).

表2 尋找差異表達(dá)基因方法在不同數(shù)據(jù)集上的性能

2.3 MAQC數(shù)據(jù)集上qRT-PCR標(biāo)準(zhǔn)基因驗(yàn)證

經(jīng)過qRT-PCR生物實(shí)驗(yàn)驗(yàn)證過的基因往往被當(dāng)做標(biāo)準(zhǔn)基因，用來評(píng)估生物信息學(xué)中各種方法的準(zhǔn)確度.通過注釋文件和倍數(shù)比率的篩選，MAQC數(shù)據(jù)集獲得了305個(gè)qRT-PCR基因，被用來評(píng)估不同尋找差異方法的準(zhǔn)確性.圖2顯示不同方法在MAQC數(shù)據(jù)集上的ROC曲線.其中圖2(a)為RobustDEA方法和其他4個(gè)原始方法的比較.結(jié)果表明，RobustDEA方法和Voom方法的ROC曲線接近重合，但是要明顯要好于DESeq2，edgeR和NOISeq 3個(gè)方法.這表明在RobustDEA方法在結(jié)合不同方法的結(jié)果之后，至少可獲得不低于原始方法中最優(yōu)性能.圖2(b)為RobustDEA方法與其他3種結(jié)合方法的比較結(jié)果.從圖中可以看出，RobustDEA方法的ROC曲線要明顯高于其他3個(gè)結(jié)合方法.Fisher，Stouffer和PANDORA方法得到比原始方法更差的結(jié)果，說明這些方法不夠穩(wěn)定，容易受到原始方法中較差結(jié)果的影響.

圖2 不同方法在MAQC數(shù)據(jù)集上的ROC曲線

圖3為原始方法和結(jié)合方法在MAQC數(shù)據(jù)集上的AUC值和F1值，從這兩個(gè)值均可說明RobustDEA方法獲得了最好準(zhǔn)確度，表明該方法相比于其他結(jié)合方法，不易受到較差性能的原始方法影響，具有更好的魯棒性.

圖3 不同方法在MAQC數(shù)據(jù)集上的AUC值和F1值

2.4 Katz和Hammer數(shù)據(jù)集的大規(guī)?；蝌?yàn)證

MAQC數(shù)據(jù)集的僅提供了305個(gè)qRT-PCR驗(yàn)證基因，缺少大量基因的評(píng)估.但是由于qRT-PCR實(shí)驗(yàn)的費(fèi)用問題，現(xiàn)實(shí)中缺少包含大量qRT-PCR驗(yàn)證基因的數(shù)據(jù)集.為了進(jìn)一步在大規(guī)?；蛏显u(píng)估RobustDEA方法的性能，采用交互驗(yàn)證方式構(gòu)建大規(guī)模差異基因數(shù)據(jù)集[11].首先分別選擇4個(gè)原始方法各自認(rèn)為最可能是差異表達(dá)的前2 000個(gè)基因.將這些基因合并成一個(gè)驗(yàn)證基因集，其中4個(gè)方法共同認(rèn)定為差異表達(dá)基因作為“真實(shí)“差異表達(dá)基因，其余為非差異表達(dá)基因.

在Katz數(shù)據(jù)集中，通過交互驗(yàn)證方式構(gòu)建的驗(yàn)證基因集包含2 501個(gè)基因，其中1 517個(gè)基因被共同認(rèn)定為差異表達(dá)基因，其余基因即被劃分為非差異基因.通過此驗(yàn)證基因集的評(píng)估，不同方法的ROC曲線如圖4.

圖4 不同方法在Katz數(shù)據(jù)集中的ROC曲線

從圖4(a)中可以看出，RobustDEA方法比4個(gè)單獨(dú)原始方法的性能都要好.從圖4(b)可以獲得4個(gè)結(jié)合方法都獲得較高的準(zhǔn)確度(圖4(b)中，PANDORA和Fisher的ROC曲線高度重合).而相比3個(gè)結(jié)合方法，RobustDEA方法同樣表現(xiàn)出最好的性能.

Hammer數(shù)據(jù)集包含兩周后和兩月后兩個(gè)時(shí)間點(diǎn)，通過交互驗(yàn)證方式分別構(gòu)建驗(yàn)證基因集.兩周后時(shí)間點(diǎn)的驗(yàn)證基因集包含2 540個(gè)基因，其中1 484個(gè)基因被共同認(rèn)定為差異表達(dá)基因.兩月后時(shí)間點(diǎn)的驗(yàn)證基因集包含2 995個(gè)基因，其中1 088個(gè)基因被認(rèn)定為差異表達(dá)基因.在驗(yàn)證基因集中其余基因?yàn)榉遣町惐磉_(dá)基因.通過兩個(gè)驗(yàn)證基因集的評(píng)估，圖5(a)和(c)顯示RobustDEA方法的ROC曲線要高于其他4個(gè)單獨(dú)原始方法，表明RobustDEA方法通過結(jié)合不同方法能進(jìn)一步提升性能.圖5(b)和(d)中與其他3個(gè)結(jié)合方法進(jìn)行比較，可以看出RobustDEA和PANDORA方法表現(xiàn)穩(wěn)定，而Fisher和Stouffer方法性能有大幅度下降.

圖5 不同方法在Hammer數(shù)據(jù)集中2個(gè)時(shí)間點(diǎn)的ROC曲線

表3為不同方法在Katz和Hammer數(shù)據(jù)集上AUC值，結(jié)果表明在3個(gè)大規(guī)模基因數(shù)據(jù)集上，RobustDEA方法均可獲得最高的性能.雖然PANDORA方法在不同數(shù)據(jù)集中會(huì)比其中某個(gè)單獨(dú)原始方法性能要略低，但是在整體上能獲得較為穩(wěn)定的結(jié)果，不會(huì)受到性能較差的NOISeq方法的影響.Fisher和Stouffer方法表現(xiàn)極為不穩(wěn)定，在Katz數(shù)據(jù)集中能獲得極高的準(zhǔn)確度，但是在Hammer數(shù)據(jù)集的2個(gè)時(shí)間點(diǎn)上性能都大幅度下降.導(dǎo)致此問題的原因在于Fisher方法對(duì)不同原始方法賦予相同的權(quán)重，容易受到性能較差方法的干擾.而RobustDEA和PANDORA方法對(duì)根據(jù)原始方法的性能進(jìn)行加權(quán)，可以降低性能較差方法的影響，從而使性能更加表現(xiàn)更加穩(wěn)定，具有較高的魯棒性.

表3 不同方法在Katz和Hammer數(shù)據(jù)集上AUC值

2.5 時(shí)間效率的分析

結(jié)合方法的評(píng)估過程分為原始方法計(jì)算，權(quán)重計(jì)算和結(jié)合評(píng)估3個(gè)部分.對(duì)于所有結(jié)合方法來說，原始方法計(jì)算所消耗時(shí)間都是一致，區(qū)別主要在于權(quán)重計(jì)算部分.PANDORA方法的權(quán)重估計(jì)是根據(jù)現(xiàn)有數(shù)據(jù)特點(diǎn)去生成模擬數(shù)據(jù)集，然后再通過計(jì)算不同原始方法在模擬數(shù)據(jù)集上的精度作為權(quán)重.這部分的計(jì)算需要花費(fèi)大量時(shí)間，導(dǎo)致PANDORA方法的時(shí)間效率大幅度下降.并且PANDORA方法的時(shí)間消耗與原始方法個(gè)數(shù)，模擬數(shù)據(jù)集中樣本量，基因個(gè)數(shù)等參因數(shù)有直接影響.RobustDEA方法的權(quán)重估計(jì)是基于原始方法計(jì)算的結(jié)果，因此能快速獲得不同方法的權(quán)重值，極大了提高了計(jì)算效率.表4為4個(gè)結(jié)合方法在MAQC數(shù)據(jù)集上的計(jì)算效率.PANDORA方法設(shè)置模擬數(shù)據(jù)集為4個(gè)測(cè)序樣本，1 000個(gè)基因.RobustDEA方法選擇N為1 000，其時(shí)間為同一工作站的CPU工作時(shí)間，性能為2.9 GHz雙核和8 G內(nèi)存，采用單核計(jì)算.從表中結(jié)果可以看出PANDORA方法花費(fèi)了整個(gè)時(shí)間894.6 s中的95%時(shí)間進(jìn)行權(quán)重計(jì)算，而RobustDEA方法只需要花費(fèi)5.5 s的時(shí)間估計(jì)權(quán)重值，大幅度提高了效率.雖然RobustDEA方法計(jì)算效率僅略高于Fisher和Stouffer方法，但其計(jì)算精度和穩(wěn)定性要明顯高于這兩個(gè)方法.

表4 不同方法在MAQC數(shù)據(jù)集上計(jì)算效率

3 結(jié)論

針對(duì)不同尋找差異表達(dá)基因方法檢測(cè)得到的差異基因集不一致問題，提出了一種快速魯棒的RNA-Seq數(shù)據(jù)尋找差異表達(dá)基因的結(jié)合方法.通過人類大腦MAQC數(shù)據(jù)集和多個(gè)老鼠數(shù)據(jù)集的評(píng)估，得出如下結(jié)論：

(1) RobustDEA方法結(jié)合現(xiàn)有尋找差異表達(dá)基因方法，通過全新的權(quán)重評(píng)估方式，可快速獲得與當(dāng)前分析數(shù)據(jù)集最合適的權(quán)重值.

(2) 相比于單個(gè)差異表達(dá)基因方法和現(xiàn)有結(jié)合方法，RobustDEA方法能獲得更加準(zhǔn)確和穩(wěn)定的結(jié)果.這表明該方法具有較好的魯棒性，受到性能較差的原始方式的影響較小.

(3) 相比PANDOR方法，RobustDEA方法采用全新的權(quán)重估計(jì)方式，可大幅度提高計(jì)算效率.

今后工作將RobustDEA方法實(shí)現(xiàn)為一個(gè)可視化在線平臺(tái)，用戶只需導(dǎo)入基因讀段數(shù)據(jù)，平臺(tái)可快速計(jì)算出差異基因集，并提供相應(yīng)的分析報(bào)告，以及更加直觀的圖表展示，從而讓用戶使用更加便利.