|
American Journal of Botany 99(2): 248–256. 2012. A COMPARISON OF STATISTICAL METHODS FOR DETECTING DIFFERENTIALLY EXPRESSED GENES FROM RNA-SEQ DATA
作為一個階段性總結(jié),，該文出現(xiàn)得很及時。請看摘要： “RNA-seq技術(shù)正在快速革新基因組學(xué)研究，而RNA-seq數(shù)據(jù)的統(tǒng)計方法正在不斷發(fā)展,。適時地回顧和比較最近提出的統(tǒng)計方法可以提供一個有用的指南,，以便選擇合適的方法進行數(shù)據(jù)分析。人們對檢測基因的差異表達的能力情有獨鐘,。這里我們通過基于不同的分布模型或真實數(shù)據(jù)進行的一系列模擬,，比較了四種近期提出的統(tǒng)計方法，edgeR,，DESeq,，baySeq和一個兩步Poisson模型(TSPM)的方法。我們按照基因的顯著性排序和假陽性率控制比較了這些方法檢測差異表達基因的能力,。所有進行比較的方法都用可免費獲取的軟件來實現(xiàn),。我們還討論了這些軟件當(dāng)前可獲取版本的可用性和功能?！?/p>
首先提綱挈領(lǐng)地回顧歷史，綜述研究背景,，突出該文章的研究目的和重要性,。 1、RNA-seq與微陣列比較 “相比于基于雜交的微陣列技術(shù),，RNA-seq有若干優(yōu)勢,，包括更大的表達水平范圍，更多的信息來檢測等位基因特異的表達,，新啟動子,，新亞型，更低噪音,，更高通量,。因此，RNA-seq正準(zhǔn)備在未來幾年取代微陣列技術(shù)成為研究基因表達的主要平臺,?！?/font> 2、RNA-seq實驗與分析技術(shù)概要 “典型的RNA-seq實驗中,，一個RNA樣品被轉(zhuǎn)換成一個cDNA片段文庫,，然后在高通量商用平臺上測序，這樣的平臺諸如Illumina的Genome Analyzer,，Helicos BioSciences的HeliScope,，Applied Biosystems的SOLiD，Pacific Biosciences的SMRT,，以及Roche的454Life Sciences測序系統(tǒng),。原始數(shù)據(jù)由大量的DNA片段序列（稱為reads）組成，這些數(shù)據(jù)要經(jīng)歷一系列的分析步驟。Oshlack等(2010)提供了一個極好的分析流程的評論,，流程包括映射reads,，匯總每個基因的reads計數(shù)，歸一化和檢測差異表達基因,。其中的表1提供了每一步分析的軟件列表,。通常，RNA-seq研究產(chǎn)生的reads要基于映射到目標(biāo)基因組或de novo組裝的轉(zhuǎn)錄組的情況分配給基因或其他分類單元,。有一些RNA-seq數(shù)據(jù)的基因表達水平定量方法仍在研究中,。可變剪切轉(zhuǎn)錄本和亞型表達的復(fù)雜性使得其成為一個活躍的研究領(lǐng)域,。因為亞型檢測不是本文的重點,，有興趣的讀者可參考Hiller等(2009)和Salzman等(2011)及其估計RNA-seq數(shù)據(jù)中的亞型豐度的參考文獻?！颉俏覀冐灤┍疚氖Ｓ嗖糠侄疾捎玫囊粋€廣義術(shù)語,，可以指一個基因模型的單外顯子或外顯子的子集合或所有外顯子?；虮磉_用映射到一個基因的reads數(shù)來度量,。因此，RNA-seq產(chǎn)生了一種基因表達的離散度量,，這與微陣列技術(shù)中可被視為連續(xù)變量的熒光強度度量不同,。因此，用于分析微陣列數(shù)據(jù)的統(tǒng)計方法不能直接應(yīng)用,，而迫切需要發(fā)展合適的統(tǒng)計方法來處理海量的RNA-seq數(shù)據(jù),。” 3,、聚焦于差異表達 “檢測跨處理/條件的差異表達基因是一個關(guān)鍵步驟,，而且有時是RNA-seq數(shù)據(jù)統(tǒng)計分析的主要目標(biāo)。差異表達基因的確定有助于我們闡明基因功能,，當(dāng)細胞響應(yīng)不同的處理和條件時,。此外，檢測差異表達基因是聚類基因表達譜或檢驗基因集富集性的事先步驟,。由于RNA-seq歷史尚短并且在不斷發(fā)展,，目前還沒有可用的標(biāo)準(zhǔn)方法基于這些數(shù)據(jù)來檢測差異表達基因。很多統(tǒng)計工作者在為此而努力,。一些文章已經(jīng)發(fā)表,，更多的可能還在研究中?！?/font> 4,、文章研究內(nèi)容與結(jié)果預(yù)告 “本文中我們首先回顧了目前可用的檢測差異表達基因的方法,，包括edgeR，DESeq,，baySeq和一個基于兩步Poisson模型的方法,。我們提供了如何下載對應(yīng)的包或代碼以在R軟件中應(yīng)用這些方法的信息。然后我們在模擬真實數(shù)據(jù)的各種設(shè)置下通過模擬研究,，比較了它們在基因的顯著性排序上的表現(xiàn),。我們還檢查了不同過程的假陽性率控制，對這樣一個基因組數(shù)據(jù)分析中的高維檢驗問題來說是一個必要的步驟,?！?/font> “我們的結(jié)果表明baySeq在低假陽性條件下具有最高的真陽性率。還發(fā)現(xiàn)相比之下,，TSPM在樣本大小為2時不能像其他方法一樣有效執(zhí)行,。在這些方法的FDR控制中，我們發(fā)現(xiàn)真實的FDR事實上比期望值大很多,。我們對這些方法的優(yōu)劣的研究是透明的,，這可能對科學(xué)家分析將來從RNA-seq研究中得到的數(shù)據(jù)是有用的?！?/font>
第二部分是基于RNA-seq數(shù)據(jù)檢測差異表達基因的若干方法的回顧 “要檢測基因的差異表達,，統(tǒng)計假設(shè)檢驗已經(jīng)準(zhǔn)備好了?；谡龖B(tài)分布已經(jīng)發(fā)展了很多統(tǒng)計方法用于歸一化微陣列的基因表達度量。例如,，廣泛使用的修正t檢驗就是R包limma基于正態(tài)假設(shè)實現(xiàn)的,。如前所述，RNA-seq技術(shù)產(chǎn)生了基因表達的離散度量,，因此基于正態(tài)分布開發(fā)的統(tǒng)計方法不能直接應(yīng)用,。對數(shù)變換可能會使高度扭曲的離散RNA-seq數(shù)據(jù)更接近正態(tài)分布。不過,，這就不得不添加一個任意的小數(shù)給一些樣本中的那些零計數(shù)基因以完成對數(shù)轉(zhuǎn)換,。即便如此，轉(zhuǎn)換后的數(shù)據(jù)仍可能不太吻合正態(tài)分布,。研究者不再聚焦于尋求轉(zhuǎn)換以使現(xiàn)存微陣列分析方法可以被應(yīng)用,，而是基于可直接對基因計數(shù)建模的離散分布開發(fā)出了一些方法?！?/p> 有三個離散概率分布被提出來對RNA-seq研究的計數(shù)數(shù)據(jù)進行建模：二項分布,、Poisson分布、負二項分布,。數(shù)學(xué)上可以證明：如果reads數(shù)充分大（對RNA-seq來說是真的）,，而且一條read映射到一個給定基因的概率足夠小,，那么二項分布可以用Poisson分布很好滴逼近。在早期使用單來源RNA的RNA-seq研究報道中,，對大部分基因來說技術(shù)重復(fù)之間計數(shù)的分布用Poisson分布擬合得很好,。但是Poisson分布的一個性質(zhì)是方差等于均值。當(dāng)有生物學(xué)重復(fù)時,，RNA-seq數(shù)據(jù)會表現(xiàn)出比Poisson分布期望的更高的變異性,，對相當(dāng)多的基因來說方差可能超過均值。這種現(xiàn)象叫做過離散,。對過離散數(shù)據(jù),，基于Poisson的分析容易因取樣誤差低估而產(chǎn)生高假陽性率。這里可用擬似然方法,，因為它引入了一個關(guān)于方差的尺度因子以允許它不同于均值,。假設(shè)一個負二項模型而不是Poisson模型是處理過離散數(shù)據(jù)的另一方法，因為負二項分布可以設(shè)定方差大于均值,。因為在研究生物學(xué)上有意義的結(jié)果時生物學(xué)重復(fù)是至關(guān)重要的,，我們希望所有試驗都能設(shè)計成包含生物學(xué)重復(fù)。因此,，我們只回顧目前能處理過離散的可用方法,。 應(yīng)該提到的是，為Serial analysis of gene expression（SAGE）數(shù)據(jù)而開發(fā)的方法可應(yīng)用于RNA-seq數(shù)據(jù)分析,。例如,，edgeR方法最初就是為SAGE數(shù)據(jù)分析開發(fā)的，現(xiàn)在用于RNA-seq數(shù)據(jù)分析,。對用于分析SAGE數(shù)據(jù)的方法的綜合分析超出了本文范圍,，有興趣的讀者可參考Lu等(2005)，Robinson & Smyth（2008）和Baggerly等(2004). 基于Poisson分布 最近,，Auer和Doerge(2011)提出了一種基于一個兩步Poisson模型的方法,。它們的原理是一些基因可能有過離散而另一些可能沒有。因此,，該方法的第一步是檢驗每個基因的過離散,。如果檢驗表明有過離散，則在第二階段應(yīng)用一個準(zhǔn)Poisson似然方法來檢驗差異表達,。否則,，在第二階段應(yīng)用基于Poisson模型的檢驗。他們對兩個基因列表分別控制FDR,，因為評估差異表達顯著性的兩種不同方法被應(yīng)用于不同的基因,。 Srivastava & Chen(2010)有一篇很有意思的文章，用了一個廣義Poisson分布對位置水平的reads計數(shù)進行建模,。GPseq實現(xiàn)的這個方法考慮了進行差異表達分析時潛在的位置偏倚,，這與我們回顧的所有其他方法不同,。我們可以獲得的真實數(shù)據(jù)集均未提供位置水平的計數(shù)，而GPseq不能處理基因水平的計數(shù),。因此,，我們的分析不包括GPseq。
基于負二項分布 有三種R包（edgeR,，DESeq,，baySeq）實現(xiàn)的方法是基于負二項模型的。edgeR所用方法是最先提出的,，原本是為SAGE數(shù)據(jù)而開發(fā)的,，SAGE數(shù)據(jù)可視為小規(guī)模的RNA-seq數(shù)據(jù)。負二項分布被用作Poisson分布的自然推廣,，只要加上一個散度參數(shù)而比Poisson分布允許額外的變異,。因為重復(fù)很昂貴，這就導(dǎo)致了RNA-seq研究的很小樣本量,，所以散度參數(shù)的估計是一個很有挑戰(zhàn)性的問題,。Robinson和Smyth(2008)提出對所有基因使用一個公共散度來達到對散度參數(shù)的更好估計。如果散度參數(shù)(度量了相比于均值額外的變異)對于所有基因相同的假設(shè)成立,，則公共散度參數(shù)可以非常精確地估計出來,，因為很多數(shù)據(jù)可用于該估計。但是,，所有基因有公共散度這一點在實踐中很少是一個合適的假設(shè),。可能更好的策略是允許不同的基因有不同的散度參數(shù),，而這些散度參數(shù)的估計可以用一些合適的統(tǒng)計方法借助基因間的信息來改進,。一些策略用在微陣列數(shù)據(jù)分析中，發(fā)展出很多檢驗借助基因間信息來更好地估計方差或平均表達與方差,。這些檢驗在與不借助這些信息的檢驗進行比較時,，被證明有很好的性能,。遵循類似的策略,，一種修正了的檢驗被提出來用于RNA-seq數(shù)據(jù)，可以用edgeR包來實現(xiàn),。 Anders和Huber(2010)也嘗試了借助基因間信息來更好地估計散度參數(shù),。他們假設(shè)了一個方差和平均表達水平之間的局部線性關(guān)系。該假設(shè)允許使用具有相似表達水平的混合數(shù)據(jù)來估計方差（或等價地過離散參數(shù)）,。這種方法在R/Bioconductor包DESeq中實現(xiàn),。edgeR和DESeq都提供基于精確檢驗或精確檢驗的逼近的檢驗p-value。 Hardcastle和Kelly (2010)提出的方法也假設(shè)數(shù)據(jù)服從負二項分布,，但在顯著性估計上與另兩種方法不同,。他們遵循一個經(jīng)驗Bayes方法,，將基因按照后驗概率估計進行排序，該模型對每個基因定義了差異表達,。關(guān)于負二項分布參數(shù)的先驗概率是利用數(shù)據(jù)找到的,。相互之間表現(xiàn)相似的樣本應(yīng)該對潛在的基因參數(shù)具有相同的先驗分布，而表現(xiàn)不同的樣本應(yīng)該具有不同的先驗分布,。該方法在R/Bioconductor包中已實現(xiàn),。
歸一化 對于使用RNA-seq數(shù)據(jù)來比較樣本間的表達，進行歸一化來調(diào)整不同的測序深度和潛在的其他重復(fù)間的技術(shù)性影響,，在上面的四種方法中都需要歸一化,。一個例子是用每個樣本的總reads數(shù)和基因長度來歸一化reads計數(shù)。該方法用RPKM來定量轉(zhuǎn)錄本水平,。但是在進行樣本間相同基因的差異表達分析而不在基因之間比較時,，相對于基因長度的歸一化就不重要了（也就是說，這種偏倚對不同樣本中會以同樣的方式影響相同的基因）,。如果不考慮基因長度,，目前可用的歸一化方法可以用一些相對于平均表達水平的尺度因子來進行。最簡單最常用的歸一化因子是文庫的總reads數(shù),，這是考慮到樣本測序越深則每個基因分配到的reads越多,。但是，通常reads總數(shù)主要是一小部分大量表達的基因貢獻的,。如果這小部分基因是差異表達的,，則使用reads總數(shù)就會極大地影響檢測差異表達的結(jié)果。Bullard等(2010)比較了幾種歸一化方法并發(fā)現(xiàn)在每個lane內(nèi)使用非零計數(shù)分布的75th分位數(shù)作為歸一化因子是一種更穩(wěn)健的選擇,，相比于標(biāo)準(zhǔn)的總計數(shù)歸一化來說,，總體表現(xiàn)是所研究的這些方法中最優(yōu)的。R包DESeq通過scaled counts的中位數(shù)來估計歸一化因子,，這是一種與75th分位數(shù)歸一化類似的思想,。R/Bioconductor包edgeR使用了對數(shù)表達比率的加權(quán)截斷均值（trimmed mean of M values，TMM）,，這是另一種穩(wěn)健的歸一化方法,。根據(jù)我們的經(jīng)驗，75th分位數(shù)和TMM方法的性能相似,。
第三部分是模擬結(jié)果 模擬數(shù)據(jù)的好處是我們知道產(chǎn)生數(shù)據(jù)的真實潛在機制,；因此，我們能評估結(jié)果,，例如,，一個宣稱的陽性事實上是正確的（真陽性）。 為baySeq和TSPM,，我們遵循Bullard等(2010),、估計歸一化因子為計數(shù)的第三個四分位數(shù),。而DESeq和edgeR自有方法估計歸一化因子。 我們基于兩個標(biāo)準(zhǔn)比較和評價不同統(tǒng)計方法的分析結(jié)果,。首先,，看基因的顯著性排序。選取不同的閾值得到TPR,、FPR對,，作出receiver operating characteristic (ROC) curve。其次,，比較FDR,。其中baySeq沒有FDR，不參與此項比較,。 文章中共有四次模擬,，每次模擬都超過100個數(shù)據(jù)集，每個數(shù)據(jù)及不少于1萬個gene,。四次模擬各有側(cè)重,，分別是一半Poisson混合一半overdispersed Poisson模擬數(shù)據(jù)、基于玉米研究真實Illumina數(shù)據(jù)生成NB參數(shù)再產(chǎn)生模擬數(shù)據(jù),、同樣基于玉米真實數(shù)據(jù)用經(jīng)驗分布得到NB參數(shù),、基于幼成淋巴細胞系真實數(shù)據(jù)模擬差異表達數(shù)據(jù)。 總體上都表明baySeq表現(xiàn)最好,。
第四部分是數(shù)據(jù)分析 edgeR和DESeq吻合度較高（86%,、91%），而與TSPM一致性較差（70?G未被DESeq和edgeR檢測到）,。
第五部分是討論 模擬數(shù)據(jù)的分析能給事件中選擇合適的方法分析RNA-seq數(shù)據(jù)提供有益指導(dǎo),。ROC曲線結(jié)果表明baySeq優(yōu)于edgeR、DESeq和TSPM,，而edgeR和DESeq性能相似,，且與baySeq接近。TSPM方法可變性大,，對小樣本產(chǎn)生的結(jié)果不好,，但是樣本量增加時性能會改善很多。 FDR控制結(jié)果常常比真實的FDR要高,。需要研究為什么沒控制好,，開發(fā)更好地FDR控制方法,。建議實踐者采用更嚴格的FDR控制方法來避免太多的假發(fā)現(xiàn),。 對于處理不同的實驗設(shè)計，所用R包的靈活性各異,。都能比較完全隨機化設(shè)計的兩組比較,。而有一些提供允許更復(fù)雜的實驗設(shè)計,，還有的允許不同的估計模式。baySeq可以分析涉及多個處理組的實驗設(shè)計,。edgeR可用于兩組或多組,，至少一個組有重復(fù)。但是差異分析只能成對比較,。edgeR有兩種方法估計散度參數(shù)——普通的和tag層次的,。除了DESeq可以沒有重復(fù)外，所有方法都要求至少一個重復(fù),。但是不推薦無重復(fù)的實驗設(shè)計,。 鑒于這些選項，方法的選取還取決于實驗設(shè)計,。 |
|
|
