新一代高通量RNA測序數(shù)據(jù)的處理與分析-

1、 新一代高通量RNA測序數(shù)據(jù)的處理與分析*王曦1汪小我1王立坤1,2馮智星1張學工1*(1生物信息學教育部重點實驗室、清華信息科學與技術國家實驗室(籌生物信息學研究部、清華大學自動化系,北京100084;2吉林大學計算機科學與技術學院,長春130012摘要隨著新一代高通量DNA測序技術的快速發(fā)展,RNA測序(RNA-seq已成為基因表達和轉錄組分析的新的重要手段. RNA-seq技術產(chǎn)生的海量數(shù)據(jù)為生物信息學帶來了新的機遇和挑戰(zhàn). 有效地對測序數(shù)據(jù)進行針對性的生物信息學處理和分析,成為RNA-seq技術能否在科學探索中發(fā)揮重大作用的關鍵. 本文以新一代Illumina/Solexa測序平臺所產(chǎn)

2、生的數(shù)據(jù)為例,在扼要介紹高通量RNA-seq測序流程的基礎上,對RNA-seq數(shù)據(jù)處理和分析的方法和現(xiàn)有軟件做一個較為全面的綜述,并對其中有待進一步研究的問題進行展望.關鍵詞 高通量RNA測序,轉錄組,基因表達,數(shù)據(jù)處理與分析,生物信息學學科分類號 Q5(生物化學,Q6(生物物理學,Q7(分子生物學*國家自然科學基金資助項目(60702002、60721003、30873464、60905013,東南大學生物電子學國家重點實驗室開放研究基金資助課題.*通訊聯(lián)系人.Tel: 010-*, E-mail: zhangxg收稿日期:(以本文投稿日期為準,接受日期:(此項由編輯部填寫- 2 -生物化學

3、與生物物理進展 Prog. Biochem. Biophys.近年來,新一代高通量測序技術得到了突飛猛進的發(fā)展,在此基礎上,高通量RNA測序即RNA-seq1-5也迅速發(fā)展. 與基因芯片技術相比,RNA-seq無需設計探針,能在全基因組范圍內(nèi)以單堿基分辨率檢測和量化轉錄片段,并能應用于基因組圖譜尚未完成的物種6,具有信噪比高、分辨率高、應用范圍廣等優(yōu)勢,正成為研究基因表達和轉錄組的重要實驗手段.RNA-seq為組學的研究帶來了高分辨率的海量數(shù)據(jù),如何有效處理和分析這些海量數(shù)據(jù)成為這一新技術能否帶來新的科學發(fā)現(xiàn)的關鍵,一些生信息學方法與軟件也應運而生. 本文針對當前RNA-seq應用的現(xiàn)實情況,

4、嘗試以Illumina/Solexa測序平臺產(chǎn)生的mRNA-seq數(shù)據(jù)為例,對RNA測序數(shù)據(jù)的產(chǎn)生過程,數(shù)據(jù)處理和分析的基本流程、關鍵方法和現(xiàn)有軟件進行較全面的介紹,并討論RNA-seq數(shù)據(jù)分析中存在的挑戰(zhàn). 1高通量測序技術簡介誕生于20世紀70年代的Sanger法是最早廣泛應用的DNA測序技術7,也是完成人類基因組計劃的基礎. 但是,由于它測序通量低,費時費力,科學家們一直在尋求通量更高、速度更快、價格更便宜、自動化程度更高的測序技術. 自2005年以來,以Roche公司的454技術、Illumina公司的Solexa技術和ABI公司的SOLiD技術為標志的新一代測序技術相繼誕生8. 新一

5、代測序技術又稱作深度測序技術,主要特點是測序通量高、測序時間和成本顯著下降9.把這種高通量測序技術應用到由mRNA逆轉錄生成的cDNA上,從而獲得來自不同基因的mRNA片段在特定樣本中的含量,這就是mRNA測序或mRNA-seq. 同樣原理,各種類型的轉錄本都可以用深度測序技術進行高通量定量檢測,統(tǒng)稱作RNA-seq或RNA測序. 目前,在已經(jīng)推出的幾種新一代測序平臺中, Illumina/Solexa測序平臺上的RNA-seq應用最廣,我們以此為例來綜述RNA-seq數(shù)據(jù)處理和分析的生物信息學問題和方法.Illumina/Solexa測序技術的基本原理是邊合成邊測序(Sequencing b

6、y Synthesis, SBS10-12,即測序過程是以DNA單鏈為模板,在生成互補鏈時,利用帶熒光標記的dNTP發(fā)出不同顏色的熒光來確定不同的堿基. 新加入的dNTP的末端被可逆的保護基團封閉,既保證單次反應只能加入一個堿基,又能在該堿基讀取完畢后,將保護基團除去,使得下一個反應可繼續(xù)進行. 為了增加熒光強度,使之更易被成像系統(tǒng)所采集,該技術在測序之前還需要對待測片段做橋式擴增(bridge amplification13( 初期的Illumina/Solexa測序技術只能在較短的測序讀長上(2030堿基保證較高的正確率. 隨著技術的改進,目前的讀長已經(jīng)增加到100堿基以上. 同時,隨著雙

7、端測序(paired-end, PE技術的成熟,測序長度更可達到單端測序的兩倍,測序通量也隨之增加. 這種測序技術是Solexa公司發(fā)展起來的,2007年被Illumina公司收購,因此現(xiàn)在通常被稱為Illumina/Solexa測序技術. 近來兩年來,Illumina/Solexa測序平臺不斷升級,相繼推出了GA(Genome Analyzer、GA IIx、HiSeq 2000等測序儀. 更多關于高通量測序平臺的介紹,可以查閱相關文獻9, 14-16.2RNA-seq測序文庫制備和測序平臺數(shù)據(jù)輸出本小節(jié)針對Illumina/Solexa測序平臺,對RNA測序文庫制備標準和平臺底層數(shù)據(jù)產(chǎn)生做

8、一個簡單的介紹.2.1 RNA-seq測序文庫制備對于mRNA-seq實驗,從總RNA到最終的cDNA文庫制備完成主要包括以下步驟. 首先,用Poly(T寡聚核苷酸從總RNA中抽取全部帶Poly(A尾的RNA,其中的主要部分就是編碼基因所轉錄的mRNA. 將所得RNA隨機打斷成片段,再用隨機引物和逆轉錄酶從RNA片段合成cDNA片段. 然后,對cDNA片段生物化學與生物物理進展 Prog. Biochem. Biophys.- 3 -進行末端修復并連接測序接頭(adapter,得到將用于測序的cDNA. 在以上過程,將RNA隨機片段化和采用隨機引物進行反轉錄,都是為了使所得cDNA片段較均勻地

9、取自各個轉錄本. 為提高測序效率,一般還需要用電泳切膠法獲取長度范圍在200 bp(±25 bp的cDNA片段,再通過RCR擴增,得到最終的cDNA 文庫.在上述文庫制備過程中,如果不是只抽取帶Poly(A尾的RNA,而是使用全部的RNA,則RNA-seq測得的就是細胞中的全部轉錄本;如果把帶Poly(A尾的RNA過濾掉,也可以得到非編碼的RNA轉錄本;如果從總RNA中只提取長度為 2123個堿基左右的RNA,則得到全部的miRNA(microRNA轉錄本,相應的方法也稱作miRNA-seq.樣品制備最終得到的是雙鏈cDNA文庫. 在后續(xù)測序中,測得的每個讀段(read隨機地來自雙鏈

10、cDNA 的某一條鏈,從讀段序列本身無法得知它是與RNA方向相同還是倒轉互補,在后續(xù)的讀段定位時需要兩個方向都考慮. 在新基因識別等應用中,轉錄本的方向對基因注釋尤為重要,需要在文庫制備和測序中保留RNA的方向信息. 最近有文獻報道了保留方向信息的RNA-seq樣品制備方法17-20.2.2 測序平臺數(shù)據(jù)輸出將RNA-seq測序文庫加入流動槽(flow cell中的各通道(lane,在橋式PCR擴增后,就可以進行測序了. 測序過程中,計算機軟件同步地對熒光圖像數(shù)據(jù)進行處理,通過分析熒光信號來確定被測堿基,并給出質量評分. 按照圖像上的位置坐標,計算機程序將同一位置測得的堿基根據(jù)測序順序連成讀段

11、(read. 由于熒光圖像文件所占有的磁盤空間很大,通常GA IIx平臺一次實驗能就產(chǎn)生上太字節(jié)(TB的圖像文件,所以一般情況下不予保留原始的熒光圖像數(shù)據(jù),而是只保留程序讀出的讀段數(shù)據(jù)及對應的質量分值,這就是多數(shù)實驗室委托測序中心進行RNA-seq測序后得到的最原始的數(shù)據(jù).為了便于測序數(shù)據(jù)的發(fā)布和共享,高通量測序數(shù)據(jù)以FASTQ格式來記錄所測的堿基讀段和質量分數(shù). 如圖1所示,FASTQ格式以測序讀段為單位存儲,每條讀段占4行,其中第一行和的第三行由文件識別標志和讀段名(ID組成(第一行以“”開頭而第三行以“+”開頭;第三行中ID可以省略,但“+”不能省略,第二行為堿基序列,第四行為對應的測序

12、質量分數(shù). 關于FASTQ格式更多地介紹可參考文獻21. 為方便保存和共享各實驗室產(chǎn)生的高通量測序數(shù)據(jù),NCBI、EBI、DDBJ等數(shù)據(jù)中心建立了大容量的數(shù)據(jù)庫SRA (Sequence Read Archive, /Traces/sra來存放共享的測序數(shù)據(jù)22, 23. 圖1 讀段FASTQ數(shù)據(jù)格式示例Fig.1 FASTQ format examples3RNA-seq數(shù)據(jù)的基本處理RNA-seq的基本應用是測量一個樣本中的基因表達或轉錄組. 有實驗表明,新一代高通量測序的技術重復數(shù)據(jù)之間的相關度較高(R20.961, 2,因此,如果對同

13、一樣本在多個通道上進行了RNA測序的技術重- 4 -生物化學與生物物理進展 Prog. Biochem. Biophys.復,我們建議可以把幾個通道的數(shù)據(jù)進行合并,這樣等效地增加了測序深度. 本節(jié)討論單個樣本RNA測序數(shù)據(jù)的基本處理流程,如圖2(a所示. 圖2 RNA-seq數(shù)據(jù)處理和分析流程圖Fig.2 The flowchart for RNA-seq data processing and analysis (a為RNA-seq數(shù)據(jù)的基本處理,其方法介紹見正文第3小節(jié). (b為兩類樣本RNA-seq數(shù)據(jù)比較分析的框架,對應于正文的第4小節(jié). (b中虛線框內(nèi)為(a所示的流程,虛線箭頭表示可

14、選輸入. 3.1 讀段定位獲得RNA-seq的原始數(shù)據(jù)后,首先需要將所有測序讀段通過序列映射(mapping定位到參考基因組上, 這是所有后續(xù)處理和分析的基礎. 在讀段定位之前,有時還需要根據(jù)測序數(shù)據(jù)情況對其做某些基本的預處理. 例如,過濾掉測序質量較差的讀段,對miRNA測序讀段數(shù)據(jù)去除接頭序列等. 高通量測序的海量數(shù)據(jù)對計算機算法的運行時間提出了很高的要求. 針對諸如Illumina/Solexa等測序平臺得到的讀段一般較短、且插入刪除錯誤較少等特點,人們開發(fā)了一些短序列定位算法. 這些算法主要采用空位種子索引法(spaced-seed indexing或Burrows-Wheeler轉換

15、(Burrows-Wheeler Transform,BWT技術來實現(xiàn)24. 空位種子索引法首先將讀段切分,并選取其中一段或幾段作為種子建立搜索索引,再通過查找索引、延展匹配來實現(xiàn)讀段定位,通過輪換種子考慮允許出現(xiàn)錯配(mismatch的各種可能的位置組合.BWT 方法通過B-W轉換25將基因組序列按一定規(guī)則壓縮并建立索引,再通過查找和回溯來定位讀段,在查找時 可通過堿基替代來實現(xiàn)允許的錯配.表1列出了目前可免費下載使用的部分短序列定位軟件. 其中采用空位種子片段索引法的代表是Maq26,而采用Burrows-Wheeler轉換的代表是Bowtie27. 總的來說,采用BWT 的定位算法在時間

16、效率上要優(yōu)于空位種子片段索引法24, 28. 隨著讀長的增加,允許讀段序列中存在插入刪除(indel的定位變得可行而重要. 由于以上兩類方法對序列中插入刪除的處理較為困難,近來人們開發(fā)了一些基于改進的Smith-Waterman動態(tài)規(guī)劃算法29的序列比對工具,如BFAST30、SHRiMP31、Mosaik 生物化學與生物物理進展 Prog. Biochem. Biophys.- 5 -(/marthlab/Mosaik等,但算法速度較慢,大多需采用計算機并行編程技術來解決運行時間的問題.表1 適用于Illumina/Solexa測序平臺的

17、讀段定位軟件Table 1 Mappers/aligners for Illumina/Solexa sequencing data名稱SAM* 質量* 主要采用技術網(wǎng)址MAQ26否是空位種子Bowtie27是是BWT http:/bowtie-BWA32是是BWT http:/bio-ZOOM33否否空位種子ELAND 否否空位種子 SOAP234否否BWT RazerS35否否q-grams過濾http:/www.seqan.de/projects/razers.htmlNovoalign 是是Needleman-Wunsch算法SHRiMP31否是空位種子q-grams過濾Smith-W

18、aterman算法/shrimpBFAST30是是Smith-Waterman算法并行編程/index.php/BFASTMosaik 是是Smith-Waterman算法并行編程/marthlab/Mosaik*SAM:是否能以SAM格式輸出;*質量:是否提供讀段定位質量信息;BWT:Burrows-Wheeler轉換.在RNA測序數(shù)據(jù)的基因組定位中,一個特殊的問題是跨越兩個外顯子接合區(qū)的讀段(junction reads的定

19、位. 在真核生物中,成熟的mRNA是經(jīng)過由mRNA前體中的外顯子經(jīng)過剪接形成的. 如果一個讀段跨越了兩個外顯子,那么就無法將這個讀段完整地定位到基因組序列上. 而同時,這種跨兩個外顯子的讀段在分析轉錄本的剪接形式和研究選擇性剪接中有重要的作用. 為了解決這一問題,人們采取兩種典型的策略來進行接合區(qū)讀段的定位:一是根據(jù)已知的基因外顯子注釋,構建所有可能的外顯子接合區(qū)序列,與基因組序列一并作為定位的參考基因組;二是不依賴基因注釋,而是先利用能完整定位到基因組的讀段得到粗略的外顯子區(qū)域,并結合剪接位點序列構建出可能的剪接位點,然后將不能完整定位的讀段分段定位到兩個外顯子可能的結合區(qū)域. Illumi

20、na/Solexa平臺提供的RNA-seq軟件分析包GApipeline采用了第一種策略. 采用第二種策略的軟件有Tophat36和G-Mo.R-Se37等,最新的Tophat軟件增加了利用已知外顯子邊界注釋信息的選項.不論是哪種測序平臺,測序中都不可避免地存在一定的錯誤,基因組中又存在單核苷酸多態(tài)性等引起的序列變化,所以在讀段定位時通常允許一定數(shù)量的錯配,可以根據(jù)不同應用調(diào)節(jié)允許錯配的程度. 另一方面,由于基因組中重復序列和高相似度序列的影響,某些讀段會出現(xiàn)定位到基因組多個位置的情況. 這些因素影響了各個讀段到基因組的定位質量,在一些新的讀段定位算法中,同時給出每個讀段與基因組匹配質量. 通

21、常在后續(xù)處理前,人們將多定位的讀段都過濾掉,也有人嘗試用適當?shù)牟呗园讯喽ㄎ蛔x段“分配”到其中某些位置上2, 38.讀段定位到基因組后推薦采用SAM(Sequence Alignment/Map格式或其二進制版本BAM格式39來存儲. 二進制版本可大大節(jié)省存儲空間,但不能直接用普通文本編輯工具顯示. 關于SAM格式的詳細介紹,可查閱(3.2 基因表達水平估計在深度測序技術出現(xiàn)之前,高通量測量不同基因表達水平的主要手段是基因芯片,在此基礎上可以對不同組織或者不同發(fā)育階段的基因表達差異和模式進行分析. mRNA-seq 數(shù)據(jù)最基本的應用也是檢測基因的表達水平,與基因芯片數(shù)據(jù)相比,RNA 測序得到的是

22、數(shù)字化的表達信號,具有靈敏度高、分辨率高、無飽和區(qū)等優(yōu)勢40-42.RNA 測序數(shù)據(jù)是對提取出的RNA 轉錄本中隨機進行的短片段測序,如果一個轉錄本的豐度高,則測序后定位到其對應的基因組區(qū)域的讀段也就多,可以通過對定位到基因外顯子區(qū)的讀段計數(shù)來估計基因表達水平. 很顯然,讀段計數(shù)除了與基因真實表達水平成正比,還與基因長度成正比,同時也與測序深度即測序實驗中得到的總讀段數(shù)正相關. 為了保持對不同基因和不同實驗間估計的基因表達值的可比性,人們提出了RPM 和RPKM 的概念2. RPM (Reads Per Million reads 即每百萬讀段中來自于某基因的讀段數(shù),考慮了測序深度對讀段計數(shù)的

23、影響. RPKM (Reads Per Kilo bases per Million reads 是每百萬讀段中來自于某基因每千堿基長度的讀段數(shù),公式表示為:910RPKM =××基因區(qū)讀段計數(shù)基因長度測序深度. RPKM 不僅對測序深度作了歸一化,而且對基因長度也作了歸一化,使得不同長度的基因在不同測序深度下得到的基因表達水平估計值具有了可比性,是目前最常用的基因表達估計方法. 軟件rSeq 43、DEGseq 軟件包44和Cufflinks 45等都提供了用上述方法進行基因表達水平計算的功能.根據(jù)RNA-seq 文庫制備標準,在不考慮基因結構的理想情況下,讀段會均勻的分

24、布在基因上. 而實際上,通過對實際數(shù)據(jù)的可視化分析很容易發(fā)現(xiàn),讀段在基因上的分布有著自身的一些模式,呈現(xiàn)出不均勻性(見圖3所示的例子. 這一問題已經(jīng)引起很多學者的關注46-48. 造成讀段分布出現(xiàn)偏好的原因可能有多個方面:在制備cDNA 文庫時,反轉錄所采用的隨機引物對RNA 序列具有一定的偏好性,使得cDNA 片段不能夠完全均勻地取自各轉錄本;在PCR 擴增中,擴增效率與序列的GC 含量等特征相關,可導致GC 含量高的cDNA 片段在文庫中拷貝數(shù)增加超過其他片段;舍棄多定位的讀段也可能導致讀段的非均勻分布;等等. 如果能對讀段分布的不均勻性進行建模并加以校正,可以提高RNA-seq 推斷基因

25、表達量的準確度. 但根據(jù)對實際數(shù)據(jù)的觀察,對于較長轉錄本,讀段非均勻分布帶來的誤差很大程度上可相互抵消,用RPKM 來估計基因的表達水平可以得到比較滿意的結果.3.3 選擇性剪接事件識別和剪接異構體表達水平推斷在真核生物中,選擇性剪接現(xiàn)象普遍存在. 基因轉錄形成的mRNA 前體(pre-mRNA 在剪接過程中因去掉不同的內(nèi)含子區(qū)域或保留不同的外顯子區(qū)域,可形成不同的剪接異構體. 根據(jù)RNA-seq 原理,只要測序深度足夠深,就能檢測到所有轉錄本的全部序列,包括來自剪接接合區(qū)的序列. 利用考慮到接合區(qū)的讀段定位方法,就有可能系統(tǒng)地研究某一組織或某一條件下的基因選擇性剪接事件.前面已經(jīng)提到,Top

26、hat 等軟件定位剪接接合區(qū)讀段的策略能標定出剪接事件中的兩個剪接位點:供體位點和受體位點. 通過比較供體位點和受體位點的組合,就能識別選擇性剪接事件4, 49. 圖3中包含了選擇性剪接識別的一個例子. 進一步,通過對供體和受體位點的讀段計數(shù),結合外顯子其他區(qū)域的讀段數(shù)據(jù),還能定量地計算選擇性剪接事件之間的比例50, 51.對于每一個剪接異構體,RNA-seq 數(shù)據(jù)能在一定程度上推斷其表達水平. 比如,可以根據(jù)已知外顯子組成和各外顯子長度對剪接異構體建立數(shù)學模型,在測序讀段在轉錄本上均勻分布的假設下,利用各外顯子上的讀段數(shù)和接合區(qū)讀段數(shù)求解異構體的表達值. Jiang 等人的方法43及軟件Is

27、oInfer 52和cufflinks 45都采用了這種思路來實現(xiàn)剪接異構體的表達推斷. 需要指出的是,某些形式的剪接異構體表達水平在這種方法框架下不可辨識53.3.4 新基因的檢測在對RNA-seq 數(shù)據(jù)的分析中,人們發(fā)現(xiàn)往往不是所有讀段都能定位到已有注釋的基因區(qū),說明除了轉錄噪聲或測序錯誤等的影響外,可能還存在尚未被注釋的基因. 這里,我們把這種尚未注釋的基因稱為新基因,包括新的蛋白質編碼基因和非編碼RNA基因. 能檢測新基因,尤其是低表達基因是RNA-seq技術優(yōu)于基因芯片的特點之一,因為它不需要利用已知基因注釋來設計檢測探針.RNA-seq技術靈敏度高,但樣品污染、測序錯誤等仍可能帶來

28、背景噪聲. 從基因組未注釋區(qū)域的RNA 測序讀段信號中檢測新基因是典型的信號檢測問題. 如何控制新基因識別的誤發(fā)現(xiàn)率(FDR是檢測方法的關鍵. Useq軟件包54將ChIP-seq數(shù)據(jù)分析的方法移植到RNA-seq數(shù)據(jù)上,用滑窗的方法來識別測序讀段定位富集的區(qū)域,給出反映滑窗所在區(qū)域讀段富集顯著程度的p值(p-value及新基因誤發(fā)現(xiàn)率,通過設定p 值或誤發(fā)現(xiàn)率的閾值,可篩選出讀段富集的區(qū)域,再將相鄰區(qū)域合并或根據(jù)剪接接合區(qū)讀段將相應區(qū)域連接,完成新基因的檢測.3.5 讀段的可視化及注釋對于復雜的組學數(shù)據(jù),能盡可能方便地直接觀察數(shù)據(jù)對于數(shù)據(jù)的分析和解釋都非常重要,對新一代測序數(shù)據(jù)的可視化和交互

29、展示是一個非常重要但容易被人忽視的問題. 不深入考查數(shù)據(jù)的細節(jié),而是滿足于對數(shù)據(jù)的統(tǒng)計分析,是高通量數(shù)據(jù)應用中經(jīng)常容易陷入的誤區(qū),方便有效地可視化工具能夠幫助避免這樣的誤區(qū). 表2列出了部分適用于RNA-seq數(shù)據(jù)的全基因組瀏覽器,其中比較具有代表性的有UCSC Genome Browser、CisGenome Browser和IGV(Integrative Genomics Viewer等. 這些瀏覽器具有如下特點:(1能在不同尺度下顯現(xiàn)單個或多個讀段在基因組上的位置,包括來源于剪接接合區(qū)的讀段;(2能在不同尺度下顯示不同區(qū)域的讀段豐度,以反映不同區(qū)域的轉錄水平或測序效率;(3能顯示基因及其

30、剪接異構體的注釋信息;(4能顯示其他注釋信息,例如物種間基因組序列保守性、序列GC含量等;(5能直接或間接支持SAM/BAM讀段定位數(shù)據(jù)存儲格式. UCSC Genome Browser55屬于基于網(wǎng)絡模式的全基因組瀏覽器,所有數(shù)據(jù)都需要上載到遠程服務器,經(jīng)過處理后將圖形返回客戶端顯示. 圖3中的例子就是從UCSC Genome Browser的顯示截取的. CisGenome Browser56是典型的本地版基因組瀏覽器,所有讀段數(shù)據(jù)、注釋信息都存于本地文件,因此不需要網(wǎng)絡連接,方便內(nèi)部考查數(shù)據(jù)用. IGV(/igv可以說是以上兩種模式

31、的融合,既可以從遠程服務器端下載各種注釋信息,又可以從本地加載注釋信息.表2 適用于mRNA-seq數(shù)據(jù)的全基因組瀏覽器Table 2 The genome browsers/viewers applicable to mRNA-seq data viewing名稱支持的數(shù)據(jù)格式網(wǎng)址IGV GFF/GFF3, BED, SAM/BAM,WIG, ./igvUCSC Genome Browser55BED, bigBed, BEDGRAPH,GFF, GTF, WIG, bigWig,BAM, Ci

32、sGenomeBrowser56BAR, BED, refFlat, FA, /jiangh/browser MochiView Wig, BED, GFF, FASTA, /sj/mochiview-startSeqMonk ELAND, GFF, BED, MAQ,SAMhttp:/www.bioinformatics.bbsrc.ac.uk/projects/seqmonkGambit BAM, BED, GFF2, GFF3,FASTA, VCFGenomeGraphs57Expre

33、ssion data, Annotationtracks/packages/release/bioc/html/GenomeGraphs.html以上列出的全基因組瀏覽器均可在Windows、Linux和蘋果公司的Mac OS等計算機平臺下運行.除對讀段的可視化外,用描述統(tǒng)計學方法對實驗數(shù)據(jù)進行分類別統(tǒng)計也十分重要. 例如,統(tǒng)計讀段在各個染色體上的分布情況和在注釋的外顯子、內(nèi)含子、剪接接合區(qū)、基因間區(qū)的分布情況等. 目前,已經(jīng)有一些用于測序數(shù)據(jù)注釋的生物信息學軟件,比如SAMtools39、BEDtools58等,但由于測序技術發(fā)展迅速,用戶需求因人

34、而異,用戶經(jīng)常還需要根據(jù)需求編寫一定的程序或腳本完成或完善注釋分析的任務. 對于熟悉圖形用戶界面的研究人員,還可以利用UCSC Table Browse56和Galaxy59, 60來配合完成注釋分析. 由于UCSC Table Browser集成了大量基因組尺度上的注釋信息,而Galaxy又為用戶提供了書寫簡單、接口明晰和直觀的數(shù)據(jù)處理流程,這種方法十分方便有效,也為很多學者在展示研究成果時所采用. 圖3 mRNA-seq數(shù)據(jù)可視化示例Fig.3 An example for mRNA-seq data visualization.圖示區(qū)域為人類基因CBX7. 圖中紅色表示樣本A的數(shù)據(jù),藍色

35、表示樣本B. 各軌道(track依次為:基因組坐標,樣本A的剪接接合區(qū),樣本A的讀段分布,樣本B的剪接接合區(qū),樣本B的讀段分布,UCSC基因注釋. 圖中還標識了,1因受體位點不同而形成的選擇性剪接;2基因的5端出現(xiàn)讀段的非均勻分布;3在兩個樣本中,差異表達基因的讀段信號強度不同;4在基因標注的內(nèi)含子(intron區(qū)域存在少量不連續(xù)的讀段.以上描述了對于基因組已知的物種進行RNA-seq數(shù)據(jù)處理基本流程,圖2(a給出了其主要步驟. 若研究對象尚未完成基因組測序,則需要采用讀段的從頭拼裝(de novo assembly6, 61, 62來代替讀段定位,后續(xù)流程也須做相應的調(diào)整. 若RNA-seq

36、實驗在文庫制備時保留了RNA的方向信息,則應分別研究來自正鏈和反鏈的轉錄產(chǎn)物,并通過與基因注釋比較來檢測反義轉錄本63. 最近,RNA-MATE64軟件在其分析流程中加入如此的處理策略. 此外,通過分析定位到外顯子接合區(qū)的讀段,還可以獲取轉錄本結構,這為研究基因的剪接調(diào)控機理提供了重要信息5;而利用RNA-seq數(shù)據(jù)提供的序列信息,通過與DNA序列的細致比較可分析轉錄組的序列差異(如SNP等65,從而研究等位基因的表達模式66, 67及RNA編輯68等。最后需要指出,由于miRNA在序列和結構上具有一定的特點,miRNA-seq數(shù)據(jù)的基本處理流程也與本節(jié)所述有所不同,感興趣的讀者可參考軟件工具

37、miRDeep69提供的處理策略.4多類樣本mRNA-seq數(shù)據(jù)間的比較分析很多RNA-seq實驗的目的是為了比較兩種或多種樣本中基因表達或整個轉錄組的差異,如比較癌癥組織和正常組織的轉錄組差異等. 這些差異既包括通常意義下的差異表達基因,也主要包括選擇性剪接模式的差異、剪接異構體表達的差異、非編碼轉錄本的差異等. 這些差異一般可以用一些統(tǒng)計假設檢驗方法檢測,但這種檢驗有時會受到測序深度、基因長度等因素的影響70, 71,需要對結果進行仔細分析,消除可能的混雜因素,必要時可以用讀段的絕對表達值倍數(shù)變化(fold-change來作為補充. 圖2(b給出了兩類樣本數(shù)據(jù)分析的框架.4.1 差異表達基

38、因的識別雖然新一代測序相對第一代測序的單位成本大大降低,但是,利用RNA測序進行基因表達研究的成本仍很高,因此,很多實驗室沒有條件進行樣本重復. 如果兩類樣本均沒有生物重復,例如只對兩個細胞系各進行一次mRNA樣本測序,則可以用隨機采樣模型通過假設檢驗來分析差異表達. 對于某個基因,如果一個讀段來自于這個基因,我們稱事件A發(fā)生. 對于一次RNA-seq實驗,事件A發(fā)生的概率可以用這個基因上的讀段數(shù)n除以所有基因上的讀段總數(shù)N來估計,即RPM. 事件A發(fā)生的概率反應了這個基因的表達水平. 如果要判斷某個基因在兩個樣本中的表達水平是否一致,就可以通過檢驗事件A在兩種條件下發(fā)生的概率是否一致來實現(xiàn),

39、采用似然比檢驗1、Fisher精確檢驗72以及基于MA圖的統(tǒng)計檢驗方法44等. 同樣,也可用RPKM作為統(tǒng)計量來進行假設檢驗分析,由于是比較同一個基因在兩個樣本間的差異,基因長度的影響被抵消,用RPKM和用RPM得到的結果相似. 對無生物重復的RNA-seq數(shù)據(jù)進行差異表達基因分析,已經(jīng)有幾個公開發(fā)表的軟件,包括DEGseq44、Useq54、Cufflinks45中的Cuffdiff模塊等. 圖4展示了我們開發(fā)的DEGseq軟件提供的多種差異表達基因識別方法的應用例子. 圖4 用DEGseq軟件包識別差異表達基因的結果Fig.4 The results given by DEGseq for

40、 differentially expressed gene identification.(a各基因在樣本A和樣本B中表達水平的散點圖. (b、(c、(d圖中紅點表示分別用FET、LRT和MARS方法得到的差異表達基因. FET:Fishers Exact Test,Fisher精確檢驗. LRT:Likelihood Ratio Test,似然比檢驗. MARS:MA-plot-based method with Random Sampling model,基于MA圖的隨機采樣模型.如果每一類樣本都包含了若干生物重復,如病人和正常人對照研究,則可以沿用基因芯片數(shù)據(jù)分析中的很多方法. 比如,

41、可以用T檢驗結合倍數(shù)變化的方法來分析差異表達. 如果兩類樣本具有配對的信息,也可以通過整合每對樣本分析結果來實現(xiàn). 其步驟為,先在每對樣本中識別出差異表達的基因,再尋找這若干組差異表達基因之間的相同者,或用投票的方法來為基因的差異程度打分. 針對某些RNA-seq數(shù)據(jù)生物樣本量小,R軟件包DEGseq44和edgeR73等還專門提供了基于改進模型的統(tǒng)計方法. 此外,一類將分類器與特征選擇包裹在一起的方法也同樣適用于此類問題(見 差異表達剪接異構體的識別差異表達剪接異構體的識別方法與差異表達基因的識別相似. 如果把剪接異構體看成是獨立的基因,那么前面討論的用于識別差異表達基因的方法

42、對剪接異構體完全適用. 但是,注意到來自于同一個基因的剪接異構體并不獨立,某些假設檢驗的基本條件并不滿足,得到的結果就不一定正確. 此外,由于現(xiàn)在剪接異構體表達推斷的方法還不夠成熟,加之在基因結構不可辨識的剪接異構體上作表達推斷會出現(xiàn)病態(tài)結果,差異表達剪接異構體的識別問題還處于探索的階段. 目前,在剪接異構體表達水平可辨識且讀段覆蓋度較高的基因上,BASIS74方法通過貝葉斯模型來推斷差異表達的剪接異構體.換一個角度,剪接異構體由選擇性剪接造成,如果剪接異構體的表達有差異,那么導致這些異構體的選擇性剪接事件及異構體特異的外顯子的表達也會有差異. 因此,對差異表達剪接異構體的識別可以轉變?yōu)榉治鲞x

43、擇性剪接事件和外顯子表達的差異75. 外顯子表達差異的分析可以完全利用基因表達差異的分析方法. 而剪接接合區(qū)也可以看成是一個較短的“外顯子”(長度一般與測序長度相當. 不過,由于外顯子長度較基因的長度短,對應的讀段數(shù)量較少,差異識別的敏感度會有所下降. Solas方法75就是根據(jù)類似的原理,采用統(tǒng)計學假設檢驗的方法來識別差異表達的剪接異構體.4.3 對樣本的分類分析通過統(tǒng)計方法識別出來的差異表達基因及剪接異構體能否有效地區(qū)別兩類樣本,可以通過分類分析進一步證實. 如果把每個基因(或剪接異構體的表達值作為特征,則差異表達基因(或剪接異構體的選取也就是特征篩選的過程. 把前面用統(tǒng)計方法等檢測出來的

44、差異表達基因(或剪接異構體用于分類分析,常被稱為過濾法. 另一類基于分類器的包裹法,例如R-SVM76、SVM-RFE77等,可以根據(jù)每個特征在分類器中所占的權重來篩選特征,因此也可以用于差異表達基因(或剪接異構體的識別. 分類的性能可以用交叉驗證(Cross-Validation, CV方法來評估. 需要特別注意的是,交叉驗證應該包括對特征選擇步驟的交叉驗證,防止發(fā)生信息泄露而導致評估結果過于樂觀. 具體做法是:將樣本按一定的策略分成兩份,一份(通常是樣本數(shù)多的一份用于特征選取和分類器訓練,而用余下的樣本進行分類器性能的估計;重復以上步驟多次,就得到交叉驗證錯誤率. 必要時還可以用隨機置換檢

45、驗(permutation test來推斷所得錯誤率的統(tǒng)計顯著性76. 當樣本數(shù)較小時,可以采用留一法交叉驗證(Leave-One-Out Cross-Validation, LOOCV.4.4 其他高層分析方法檢測差異表達的基因或差異表達異構體是人們認識所研究的生物問題機理的第一步,接下來需要從功能上研究這些差異轉錄現(xiàn)象的分子機理. 這與在基因芯片應用中所面臨的是同樣的生物問題,對芯片數(shù)據(jù)分析結果的后續(xù)處理方法,都可以借鑒到測序數(shù)據(jù)上來. 如何進一步地從機理來解釋結果,還需結合已知生物學知識進行后續(xù)分析.人們對基因芯片得到的基因表達數(shù)據(jù)進行分析的很多方法都可以用到RNA-seq數(shù)據(jù)上來,比如

46、利用機器學習方法進行分類和特征選擇,對差異表達的基因進行GO(gene ontology78類別富集分析、信號通路富集分析等,一些常用的分析工具包括GoMiner79、DAVID80和VisANT81等.需要說明的是,在各種以差異表達基因為基礎的分析中,由于基因表達水平都是通過讀段計數(shù)來估計的,表達水平較高或轉錄本較長的基因擁有更多的讀段,更容易被多數(shù)統(tǒng)計方法識別為差異表達基因70. 這種偏好可能對后續(xù)分析帶來影響. 以GO類別富集分析為例,這種偏好將導致長基因占主導的功能類別更有可能被識別為富集的功能. 這將對生物機理的研究帶來誤導. 最近,Young等人發(fā)展了一種GOseq方法71,針對這

47、一偏好對GO類別富集分析做了改進.5RNA-seq數(shù)據(jù)處理中的生物信息學挑戰(zhàn)高通量測序技術的發(fā)展十分迅速,這要求相應的數(shù)據(jù)處理與分析方法快速跟進. 正是這些方法,架起了高通量實驗數(shù)據(jù)與科學問題之間的橋梁. 這種橋梁作用正日趨重要,也為生物信息學帶來了挑戰(zhàn)7, 71. 這里,我們重點討論兩方面的挑戰(zhàn):一、如何實現(xiàn)剪接接合區(qū)讀段的準確定位?二、在數(shù)據(jù)處理各階段中,如何對RNA-seq數(shù)據(jù)的系統(tǒng)誤差和固有偏好建模或補償,以消除它們可能帶來的錯誤推斷及結論?5.1 剪接接合區(qū)讀段的定位測序技術的一個發(fā)展趨勢是測序長度不斷增加. 隨著讀長的增加,RNA-seq中來自剪接接合區(qū)的讀段會越來越多. 我們粗略

48、估算,按照人類基因組refSeq基因注釋,一般情況下,如果測序讀長為50個堿基,則約有10%的讀段來自剪接接合區(qū);而當測序長度達到100個堿基時,這個比例將達到四分之一左右. 對這些剪接接合區(qū)讀段的分析,將使我們能夠更準確地檢測剪接事件和推斷剪接異構體的表達水平,大大推進人們對選擇性剪接的研究.在RNA-seq出現(xiàn)的早期,人們沒有意識到剪接接合區(qū)讀段的重要性. 因為當時的讀長只有2030個堿基,來自剪接接合區(qū)的讀段所占比例甚小. 當時讀段定位的通常做法是,先將讀段與全基因組序列做映射定位,再考慮不能定位的讀段否來自于剪接接合區(qū)2. 這種做法雖然在一定程度上保證了讀段定位的比率,但由于基因組中重

49、復序列和相似序列的存在,部分接合區(qū)讀段有可能在容許錯配的情況下被定位到基因組上其他位置,從而失去了定位到正確的剪接接合區(qū)的機會.在讀段定位時,如果要同時考慮基因組序列和剪接接合區(qū)序列,就要利用已知的剪接事件注釋. 這是目前軟件通用的方法. 然而,包括人類在內(nèi)的各物種的基因注釋信息都還有待完善,也沒有較完整的剪接組(splicome數(shù)據(jù)庫,能夠不依賴注釋信息和對剪接機理的現(xiàn)有認識,高效、準確地定位所有已知和未知的接合區(qū)讀段,仍然是對讀段映射定位算法的一個挑戰(zhàn).5.2 系統(tǒng)噪聲和偏好的分析雖然深度測序技術的準確性較以前的技術有了很大提高,但仍然存在錯誤和噪聲. 比如從圖3中可以看到,內(nèi)含子區(qū)內(nèi)有一

50、些不連續(xù)的讀段,很可能由系統(tǒng)噪聲造成,如樣品污染、測序錯誤和不恰當?shù)淖x段定位策略等. 從圖3還能看出,外顯子區(qū)域內(nèi)的讀段信號分布也很不均勻. 有文獻報道,序列組成尤其是GC 含量46、RNA二級結構2等也有可能是導致讀段不均勻分布的原因. 這些噪聲和分布偏好將影響新基因的識別和對剪接異構體形式和表達水平推斷.合理地建模RNA-seq數(shù)據(jù)中的系統(tǒng)噪聲和偏好是解決上述問題最有效的辦法. 基本的思路可以是:首先根據(jù)實驗原理尋找可能產(chǎn)生系統(tǒng)噪聲或偏差的因素,并盡可能將這些因素轉化成可量化的特征,如序列特征、二級結構等;然后,將用實驗數(shù)據(jù)對這些特征做統(tǒng)計分析,構造和訓練模型,用模型來對數(shù)據(jù)進行校正. 需

51、要注意的是,某些偏好是由當前的測序技術和分析方法共同造成的,難以完全消除71. 在這種情況下,后續(xù)處理和解釋時需要充分意識到這種偏好可能對生物學結論帶來的影響,必要時通過補充其他實驗來驗證和修正通過高通量測序得到的生物結論.6總結與展望本文以Illumina/Solexa測序平臺為例,嘗試對新一代測序技術的RNA-seq數(shù)據(jù)處理和分析方法做了較為全面的梳理,并對各個環(huán)節(jié)上可用的軟件進行了匯總. 高通量測序是正在飛速發(fā)展的技術,相應的生物信息學方法也在快速發(fā)展,這里討論的是RNA-seq中的一些代表性的方法和問題,希望能對正在或即將采用RNA-seq實驗進行科學研究的學者和進行RNA測序數(shù)據(jù)處理

52、的同行提供參考.RNA測序和基因芯片有很多共同的應用領域,盡管相對還不是很成熟,RNA-seq技術已經(jīng)在很多方面已經(jīng)表現(xiàn)出了優(yōu)勢,有人甚至預言基因芯片時代即將結束36. 但也有報道認為,RNA-seq數(shù)據(jù)在基因表達水平的估計上和基因芯片相比沒有明顯的優(yōu)勢60,加上測序的成本目前還遠高于芯片實驗的成本,所以更多人認為測序和基因芯片將長期共存,以各自不同的特點在現(xiàn)代組學研究中發(fā)揮作用.新一代高通量測序技術的應用面非常廣82,RNA-seq只是其中一個方面,除此之外,基因組的從頭測序和重測序83, 84、染色質免疫沉淀測序(ChIP-seq85, 86、甲基化測序(Methyl-seq87, 88等

53、技術都同樣有著廣泛的應用. 尤其是,用ChIP-seq研究蛋白質與DNA的相互作用,能夠得到高分辨率的轉錄因子結合數(shù)據(jù)和組蛋白修飾等表觀遺傳學數(shù)據(jù). 發(fā)展有效的生物信息學方法,將ChIP-seq數(shù)據(jù)與RNA-seq得到的轉錄組數(shù)據(jù)進行綜合分析,將大大推進人們對復雜的基因轉錄調(diào)控系統(tǒng)的認識.致謝感謝本實驗室劉霖曦、謝芃、孟璐等同學對本工作有意義的討論,感謝斯坦福大學Wing H. Wong教授、Hui Jiang博士和Jun Li同學等的討論和幫助.參 考 文 獻1 Marioni J C, Mason C E, Mane S M, et al., RNA-seq: an assessment

54、of technical reproducibility andcomparison with gene expression arrays. Genome Res, 2008, 18(9: 1509-1517.2 Mortazavi A, Williams B A, McCue K, et al., Mapping and quantifying mammalian transcriptomes by RNA-Seq.Nat Methods, 2008, 5(7: 621-628.3 Nagalakshmi U, Wang Z, Waern K, et al., The transcript

55、ional landscape of the yeast genome defined by RNAsequencing. Science, 2008, 320(5881: 1344-1349.4 Sultan M, Schulz M H, Richard H, et al., A global view of gene activity and alternative splicing by deepsequencing of the human transcriptome. Science, 2008, 321(5891: 956-960.5 Wang E T, Sandberg R, L

56、uo S, et al., Alternative isoform regulation in human tissue transcriptomes. Nature,2008, 456(7221: 470-476.6 Birzele F, Schaub J, Rust W, et al., Into the unknown: expression profiling without genome sequenceinformation in CHO by next generation sequencing. Nucleic Acids Res, 2010, doi:10.1093/nar/gkq116.7 Sanger F, Nicklen S, and Coulson A R, DNA sequencing with chain-terminating inhibitors. Proc Natl Acad SciU S A, 1977, 74(12: 5463-5467.8 Margulies M, Egholm M, Altman W E, et al., Genome sequencing in microfabricated high-density picolitrereactors. Nature, 2005, 437(7