Alicia Oshlack, Mark D Robinson, Matthew D Young.
From RNA-seq reads to differential expression results.Genome Biology 2010, 11:220 http://genomebiology.com/2010/11/12/220
這篇評論文章摘要只有一句話:有很多可用的方法和工具進行預處理高通量RNA-seq數(shù)據(jù)和檢測差異表達。
高通量測序技術現(xiàn)已普遍用于生物學中。這些技術產生了數(shù)以百萬計的短序列reads，常規(guī)地應用于基因組、表觀基因組和轉錄組。對樣品中穩(wěn)定狀態(tài)的RNA 進行測序，稱為RNA-seq，從先前技術的很多限制中解放出來，如對先驗的物種知識的依賴，這是微陣列和PCR所須的。此外，RNA-seq允許我們闡 明先前難以見到的轉錄組復雜性，如等位點特異的表達和新啟動子、亞型。但是，產生的數(shù)據(jù)集是大而復雜的，解釋不是直接的。正如任何高通量技術一樣，分析方 法對解釋數(shù)據(jù)是至關重要的，而RNA-seq分析過程一直在不斷演變。因此，是時候回顧當前可用的分析方法和評論問來研究方向了。
理解RNA-seq數(shù)據(jù)依賴于感興趣的科學問題。例如，決定等位點表達的差異需要精確確定轉錄的SNPs的廣泛存在。另一方面融合基因或癌癥樣品中的畸變可以通過尋找RNA-seq數(shù)據(jù)中的新轉錄本來檢測。過去一年（即2009年），一些方法涌現(xiàn)出來，用RNA-seq數(shù)據(jù)進行豐度估計，可變剪接、RNA編輯和新轉錄本的檢測。然而，很多生物學研究的基本對象是樣品間的基因表達譜。因此，在本評論中我們聚焦于可以檢測樣品間基因表達水平差異的可用方法。這種分析與控制實驗尤其相關，如比較同一組織的野生型和突變株的表達，比較處理與未處理的細胞，癌癥的和正常的細胞等等。我們在此列出用于檢測RNA-seq差異表達的處理流程，并檢查可以執(zhí)行該分析的可用方法和開源軟件。我們還突出了需要進一步研究的一些方面。
多數(shù)RNA-seq實驗取一個純化了的RNA樣品，切碎，轉換成cDNA并在高通量平臺如Illumina GA/HiSeq、SOLiD或Roche 454上測序。該過程產生了來自cDNA片段的一端的、數(shù)以百萬計的reads（25~300bp）。該過程一個常用的變式是生成雙末端reads，即paired-end reads。各平臺在化學和處理步驟上有本質不同，但忽略精確的細節(jié)后，原始數(shù)據(jù)都是由帶有質量值的短序列的一個長列表組成；這就形成了本評論的切入點。
圖1列出了差異表達分析的典型RNA-seq流程概覽。首先，reads映射到基因組或轉錄組。其次，每個樣品映射的reads依實驗目的而組裝成基因水 平、外顯子水平或轉錄本水平的表達概括。接下來，匯總的數(shù)據(jù)進行歸一化以與差異表達的統(tǒng)計檢驗相協(xié)調，產生了一個帶有P-value和倍數(shù)變化的、排好序 的基因列表。最后，執(zhí)行系統(tǒng)生物學方法從這些列表中獲得生物學的見解，就像在微陣列上進行的那樣。我們批判了下列目前可用的RNA-seq數(shù)據(jù)分析方法的 每一步。我們聚焦于普遍可用的開源軟件而不是提供一個所有工具的完整列表。
映射
為了用RNA-seq數(shù)據(jù)比較樣品間表達水平，必須把短reads轉換成表達定量。這個過程的第一步就是read映射或比對。最簡單的，映射工作就是找到 短read與參考序列已知的唯一位置。然而，真實情形中參考序列從來不是所測序RNA的實際生物源的完美表示。除了樣品特異的屬性如SNPs和 indels之外，還要考慮來自剪接過的轉錄組而非基因組的reads。而且，短read有時完美地比對到多個位置，也可能包含不得不考慮的測序錯誤。因 此，真正的任務是找到短read最佳匹配到參考序列的位置，其中允許錯誤和結構變異。
盡管對如何最佳比對reads到參考序列的研究還在進行，但是所有的解決辦法都涉及在算法的計算需求和允許匹配參考序列的模糊性之間一定的妥協(xié)。幾乎所有 的短read比對器都采用了首先通過啟發(fā)式匹配的策略，這迅速找打了可能位置的一個簡化列表，接著對候選位置進行全面評估，通過一個復雜的局部比對算法。 如果不做預先的啟發(fā)式搜索來約減潛在的比對位置數(shù)，那么在目前的硬件上執(zhí)行百萬級短reads的局部比對會是計算上不可能的。
目前的比對器能用hash表或Burrows Wheeler變換（BWT）進行快速啟發(fā)式匹配。hash表比對器對于檢測read和參考序列的復雜差異有易于擴展的優(yōu)點，在不斷增加計算需求的代價 下。而BWT比對器可以很有效地映射很接近匹配參考序列的reads，但是一旦考慮更復雜的錯配就會大幅度慢下來。這些技術的詳細說明可參考文獻 23,26-30.
比對器在怎樣處理多映射方面也很不一樣。多數(shù)比對器要么忽略多映射、隨機定位它們，要么基于局部覆蓋度的估計來定位，盡管結合比對分數(shù)的方法也已經提出。PE reads減少了多映射問題，使得多映射的模糊性在大多數(shù)情況下都可以解決。
當考慮reads來自基因組DNA時，所有要做的就是映射到一個相關的參考基因組上。但是，RNA-seq是測序轉錄組片段。這個差異可以用幾種方式處 理。既然轉錄組是建立在基因組之上，那么最常用的（至少是最初的）方法就是用基因組自身做參考序列。這有容易而不偏向任何已知注釋的好處。但是跨外顯子邊 界的reads不會映射到參考序列上。因此，用基因組做參考序列會給出較少外顯子的轉錄本以更高的覆蓋度。越長的reads越可能跨外顯子邊界，因此引起 接合reads比例增加。
為考慮接合reads，通常的實踐是建立外顯子接合位點庫，其中用注視了的外顯子的邊界構建了參考序列。為了不依賴現(xiàn)存注釋的跨外顯子邊界，可用數(shù)據(jù)集自 身來從頭檢測剪接位點。另一個選擇是轉錄組從頭組裝。所有的從頭方法都能鑒定新轉錄本，并且對沒有參考基因組或注釋的物種來說是唯一的選項。但是從頭方法 是計算密集的，需要長PE reads和高覆蓋度以可靠地進行計算。
常用的轉錄組映射方法是逐漸增加映射策略的復雜性以處理未比對上的reads。
匯總映射的reads
已經盡可能多的獲得了reads的基因組位置，下一步任務就是在一些生物學意義單位上匯總和合計這些reads，如外顯子、轉錄本、基因等水平。最簡單的 最常用的方法是計數(shù)與基因的外顯子重疊的reads。但是有相當部分reads映射到基因組上已注釋外顯子以外的區(qū)域，即使是良好注釋的物種，如小鼠、 人。
一個可選的匯總是包含沿基因全長的reads，從而結合了內含子reads。這就在匯總中包含了未注釋的外顯子并考慮了注釋不太好或可變的外顯子邊界。但 是，包含內含子也會捕獲到重疊轉錄本——它們共享一個基因組位置但是源于不同基因。還有其他很多可能的變體用于匯總，如只包含映射到編碼序列的reads 或者只匯總從頭預測的外顯子的reads。接合reads也可添加到基因匯總計數(shù)中或用于對剪接亞型的豐度進行建模。這些不同的可能性在圖2b中圖示說 明。在這些選項下，匯總的選擇可能大幅改變每個基因的reads計數(shù)，甚至比映射策略的選擇影響要大。盡管如此，很少有研究實現(xiàn)了哪種匯總方法是最適合差 異表達檢測的。
歸一化
歸一化使得可以比較樣品間和樣品內的表達水平。已經證明，歸一化是RNA-seq數(shù)據(jù)的差異表達分析的一個關鍵步驟。文庫內和文庫間比較的歸一化方法是不同的。
文庫內歸一化允許定量每個基因相對于樣品內其他基因的表達水平。因為越長的轉錄本有越高的reads計數(shù)，文庫內歸一化的常用方法是用基因長度去除匯總計 數(shù)。廣泛使用的RPKM在樣品內比較中同時解釋了文庫大小和基因長度的影響。為了驗證該方法，Mortazavi等引入了一些阿拉伯芥RNAs到小鼠的組 織樣品中，跨過一系列基因長度和表達水平。這些非天然的RNAs稱為spike-ins，說明了RPKM給出了基因間表達水平的精確比較。然而，已經證明 表達的轉錄本的read覆蓋深度是不一致的，因為序列內容和RNA制備方法，如隨機六聚體引發(fā)。把這些認識結合到文庫內歸一化方法中可能會改進比較表達水 平的能力。使用RNA-seq數(shù)據(jù)來估計樣品中轉錄本的絕對數(shù)目也是可能的，但是這需要RNA標準品和額外信息，如總細胞數(shù)和RNA制備產出率。
在樣本間檢驗單個基因的差異表達時，技術偏倚如基因長度與核酸組成，大部分會抵消，因為用于匯總的基礎序列在樣本之間是相同的。然而，樣本間歸一化對于相 對不同的文庫的比較計數(shù)來說仍是必要的。最簡單最常用的歸一化通過文庫的總reads進行調整，考慮了測序深度的影響。但是已經證明需要更聰明的歸一化來 考慮組成的影響，或這說事實上小部分高表達基因會占總序列數(shù)的相當部分。為了對這些特征進行說明，可以從數(shù)據(jù)中估計出尺度因子，用于差異表達檢驗的統(tǒng)計模 型。對于后續(xù)分析來說，尺度因子比原始計數(shù)有優(yōu)勢。另一方面，分位數(shù)歸一化和一種用匹配指數(shù)律分布的方法也被提出用于RNA-seq的樣本間歸一化。這些 變換的非線性去除了數(shù)據(jù)的計數(shù)本質，使得不清楚怎樣合適地進行差異表達檢驗。目前，分位數(shù)歸一化似乎并未改善差異表達檢測到合適的尺度因子那樣的程度，也 不清楚指數(shù)律分布應用于所有數(shù)據(jù)集的情況。
差異表達
差異表達分析的目標是突出在不同實驗條件下豐度顯著變化的基因。一般地，這意味著得到每個文庫的匯總計數(shù)數(shù)據(jù)表并進行感興趣樣本間的統(tǒng)計檢驗。
很多方法開發(fā)出來以進行微陣列數(shù)據(jù)的差異表達分析。然而，RNA-seq給出每個基因的離散度量而微陣列的強度給出了一個連續(xù)的強度分布。盡管微陣列的強 度通常是對數(shù)變換過的，而作為正態(tài)分布的隨機變量進行分析，但是計數(shù)數(shù)據(jù)的轉換并不能用連續(xù)分布很好的逼近，特別是對低計數(shù)范圍和小樣本。因此，合適于計 數(shù)數(shù)據(jù)的統(tǒng)計模型對于抽取RNA-seq數(shù)據(jù)的大部分信息是很重要的。
一般地，Poisson分布構成了對RNA-seq數(shù)據(jù)進行建模的基礎。在任何使用單RNA源的RNA-seq研究中，在一個Illumina GA測序儀的多個lane上進行測序，擬合優(yōu)度檢驗表明多數(shù)基因在lane之間的分布事實上是Poisson分布的。這一點被獨立的技術試驗驗證過，而且 已有可用的軟件工具執(zhí)行該分析。但是Poisson假設并未很好滴捕捉到生物學變異。因此對于有生物學重復的數(shù)據(jù)集進行基于Poisson的分析將傾向于 高的假陽性率，源于低估了取樣誤差。盡管RNA-seq平臺有低背景、高敏感性，但是帶有生物學重復的試驗設計對于將RNA豐度變化推廣到取樣群體中仍是 至關重要的。一般，RNA-seq試驗設計，包括分組、隨機化和重復都已經深入討論過。
為了解釋生物學變異性，為SAGE數(shù)據(jù)而開發(fā)的方法最近被用于RNA-seq數(shù)據(jù)。兩者的主要差異在于數(shù)據(jù)集的規(guī)模。為了解釋生物學變異性，使用負二項分 布作為Poisson分布的自然推廣，要估計一個額外的散度參數(shù)。一些基于負二項分布的計數(shù)數(shù)據(jù)差異表達分析變體涌現(xiàn)出來，包括普通散度模型，用加權似然 率共享所有基因的信息，均值方差關系的經驗估計和使用等價類的經驗Bayes實現(xiàn)。Poisson模型的包含散度的擴展也已提出，如廣義Poisson分 布或兩步Poisson模型（在取決于數(shù)據(jù)中散度證據(jù)的兩個模式下檢驗差異表達）。一些實時轉錄本發(fā)現(xiàn)和定量或可變亞型分析的工具也執(zhí)行差異表達分析。但 是值得注意的是這些方法要么用Poisson分布要么用Fisher 精確檢驗，顯然都不能處理上面討論的生物學變異。
當前的很多計數(shù)數(shù)據(jù)差異表達分析策略中有很多都限于簡單實驗設計，如成對或多組比較。據(jù)我們所知，還沒有針對更復雜設計的分析的通用方法提出來，如配對樣 本會時間過程實驗，在RNA-seq數(shù)據(jù)的語境下。缺少這樣的方法，研究者就把他們的計數(shù)數(shù)據(jù)轉換成有合適工具的連續(xù)數(shù)據(jù)。一般線性模型提供了對上述計數(shù) 數(shù)據(jù)的邏輯推廣，而也需要開發(fā)更聰明的策略來共享所有基因的信息；目前軟件工具提供了這些方法。進一步，上面討論的方法主要目的在于匯總已有注釋的表達水 平。以無目標方式檢測差異表達的方法近來也提出來了，如極大均值差異檢驗。
系統(tǒng)生物學：超越基因列表
在很多情形中，建立差異表達基因列表并非分析的最終步驟；可以通過尋找基因集的表達變化獲得對試驗系統(tǒng)的深入生物學見解。很多聚焦于基因集檢驗、網(wǎng)絡推斷 和知識庫的工具為分析微陣列數(shù)據(jù)集的差異表達基因而設計出來。然而，RNA-seq受到微陣列數(shù)據(jù)所沒有的一些偏倚所影響。例如，基因長度偏倚是RNA- seq數(shù)據(jù)的一個問題，其中越長的基因有越高的計數(shù)。這導致了對長的高表達基因來說，差異表達檢測有更高的檢定力。這些偏倚極大地影響下游分析結果，如 GO富集。為了能進行基因集分析，Bullard等建議通過除以基因長度的平方根來修改差異表達t-statictic以最小化差異表達中長度偏倚的影 響。另一方面，GOseq特別為RNA-seq數(shù)據(jù)開發(fā)的工具，可以把長度或總計數(shù)偏倚合并到基因集檢驗中。隨著對RNA-seq數(shù)據(jù)偏倚的認識深化，結 合了這種認識的系統(tǒng)生物學工具對于提取出生物學見解是至關重要的。
對于集成RNA-seq數(shù)據(jù)結果和其他生物學數(shù)據(jù)源已建立更完整的基因調控圖，有著廣泛的理解。例如，RNA-seq與基因分型數(shù)據(jù)結合以鑒定解釋個體間 基因表達變異的遺傳基因座（表達數(shù)量性狀位點，eQTLs）。而且，整合表達數(shù)據(jù)和轉錄因子結合、RNA干擾、組蛋白修飾以及DNA甲基化信息具有更好理 解各種調控機制的潛力。這種整合性分析的一些報道近來也出現(xiàn)了。例如，Lister等突出了基因體中CG和非CG甲基化水平與RNA-seq表達水平的顯 著差異。類似地，測序數(shù)據(jù)集的組合正開始提供單等位基因與表達、組蛋白修飾和DNA甲基化之間關聯(lián)性的見解。
Outlook
本評論中，我們列出了處理RNA-seq短reads以進行樣本間差異表達分析的主要步驟。簡言之，過程就是，映射并匯總短reads序列，然后樣本間歸 一化并執(zhí)行差異表達統(tǒng)計檢驗。進一步的生物學見解可以通過尋找基因集內表達變化模式和整合RNA-seq數(shù)據(jù)和其它來源的數(shù)據(jù)來獲得。
盡管這個流程的很多部分都是擴展研究的焦點，但仍有些領域存在進一步細化的可能。目前，很少有工作在研究那種匯總度量是最適合尋找樣本間差異表達基因的。 為了進行更復雜試驗設計的分析，還有擴展現(xiàn)有差異表達檢測統(tǒng)計方法的余地。而且，現(xiàn)有的很多方法的相對優(yōu)點應在進一步的研究中經受考驗，依照其分析各種研 究設計的靈活性，其在大大小小的研究中的性能，對測序深度的依賴和強加的假設（如均值方差關系）的準確性。進一步，盡管有很多用RNA-seq進行可變剪 接檢測的例子，但是仍有必要擴展當前方法以在生物學變異占主導時檢測基因亞型偏好的差異，也許使用上述基于計數(shù)的方法。
由于產生短reads的實驗協(xié)議之間有本質不同，正式比較RNA-seq平臺以及很多數(shù)據(jù)分析方法的相對優(yōu)點將是很重要的。這樣的研究會揭示平臺特異的差 異表達分析的好處并會促進更好的數(shù)據(jù)整合。該領域仍然相對年輕，我們希望在未來有更多的RNA-seq數(shù)據(jù)分析新方法和工具涌現(xiàn)出來。
