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