高通量rna-seq數(shù)據(jù)分析與分析的研究進(jìn)展

高通量rna-seq數(shù)據(jù)分析與分析的研究進(jìn)展

近年來，新一代高流量法測(cè)序技術(shù)取得了快速發(fā)展。在此基礎(chǔ)上，高流量?jī)?nèi)容量的研究序列（即na-seq）迅速發(fā)展。與基因芯明技術(shù)相比，ra-seq不需要設(shè)計(jì)和檢測(cè)，可以在整個(gè)重組范圍內(nèi)檢測(cè)單堿基分辨率，并量化修改段，并應(yīng)用于構(gòu)建的矩陣光譜中的后代。具有信噪比高、分辨率高、應(yīng)用范圍廣等優(yōu)點(diǎn)。已成為研究基因表達(dá)與移植組重要的實(shí)驗(yàn)手段。RNA-seq為基因組學(xué)的研究帶來了高分辨率的海量數(shù)據(jù),如何有效處理和分析這些海量數(shù)據(jù)成為這一新技術(shù)能否帶來新的科學(xué)發(fā)現(xiàn)的關(guān)鍵,一些生物信息學(xué)方法與軟件也應(yīng)運(yùn)而生.本文針對(duì)當(dāng)前RNA-seq應(yīng)用的現(xiàn)實(shí)情況,嘗試以Illumina/Solexa測(cè)序平臺(tái)產(chǎn)生的mRNA-seq數(shù)據(jù)為例,對(duì)RNA測(cè)序數(shù)據(jù)的產(chǎn)生過程及數(shù)據(jù)處理和分析的基本流程、關(guān)鍵方法和現(xiàn)有軟件進(jìn)行較全面的介紹,并討論RNA-seq數(shù)據(jù)分析中存在的挑戰(zhàn).1rna-seq/rna-seq/n-認(rèn)誕生于20世紀(jì)70年代的Sanger法是最早被廣泛應(yīng)用的DNA測(cè)序技術(shù),也是完成人類基因組計(jì)劃的基礎(chǔ).但是,由于它測(cè)序通量低,費(fèi)時(shí)費(fèi)力,科學(xué)家們一直在尋求通量更高、速度更快、價(jià)格更便宜、自動(dòng)化程度更高的測(cè)序技術(shù).自2005年以來,以Roche公司的454技術(shù)、Illumina公司的Solexa技術(shù)和ABI公司的SOLiD技術(shù)為標(biāo)志的新一代測(cè)序技術(shù)相繼誕生.新一代測(cè)序技術(shù)又稱作深度測(cè)序技術(shù),主要特點(diǎn)是測(cè)序通量高、測(cè)序時(shí)間和成本顯著下降.把這種高通量測(cè)序技術(shù)應(yīng)用到由mRNA逆轉(zhuǎn)錄生成的cDNA上,從而獲得來自不同基因的mRNA片段在特定樣本中的含量,這就是mRNA測(cè)序或mRNA-seq.同樣原理,各種類型的轉(zhuǎn)錄本都可以用深度測(cè)序技術(shù)進(jìn)行高通量定量檢測(cè),統(tǒng)稱作RNA-seq或RNA測(cè)序.目前,在已經(jīng)推出的幾種新一代測(cè)序平臺(tái)中,Illumina/Solexa測(cè)序平臺(tái)上的RNA-seq應(yīng)用最廣,我們以此為例來綜述RNA-seq數(shù)據(jù)處理和分析的生物信息學(xué)問題和方法Illumina/Solexa測(cè)序技術(shù)的基本原理是邊合成邊測(cè)序(sequencingbysynthesis,SBS),即測(cè)序過程是以DNA單鏈為模板,在生成互補(bǔ)鏈時(shí),利用帶熒光標(biāo)記的dNTP發(fā)出不同顏色的熒光來確定不同的堿基.新加入dNTP的末端被可逆的保護(hù)基團(tuán)封閉,既保證單次反應(yīng)只能加入一個(gè)堿基,又能在該堿基讀取完畢后,將保護(hù)基團(tuán)除去,使得下一個(gè)反應(yīng)可繼續(xù)進(jìn)行.為了增加熒光強(qiáng)度,使之更易被成像系統(tǒng)所采集,該技術(shù)在測(cè)序之前還需要對(duì)待測(cè)片段做橋式擴(kuò)增(bridgeamplification)2倫的seq搜索庫(kù)的準(zhǔn)備和搜索平臺(tái)數(shù)據(jù)的輸出本小節(jié)針對(duì)Illumina/Solexa測(cè)序平臺(tái),對(duì)RNA測(cè)序文庫(kù)制備標(biāo)準(zhǔn)和平臺(tái)底層數(shù)據(jù)產(chǎn)生做一個(gè)簡(jiǎn)單的介紹.2.1基于轉(zhuǎn)錄本的rna-seq基因測(cè)序?qū)τ趍RNA-seq實(shí)驗(yàn),從總RNA到最終的cDNA文庫(kù)制備完成主要包括以下步驟.首先,用Poly(T)寡聚核苷酸從總RNA中抽取全部帶Poly(A)尾的RNA,其中的主要部分就是編碼基因所轉(zhuǎn)錄的mRNA.將所得RNA隨機(jī)打斷成片段,再用隨機(jī)引物和逆轉(zhuǎn)錄酶從RNA片段合成cDNA片段.然后,對(duì)cDNA片段進(jìn)行末端修復(fù)并連接測(cè)序接頭(adapter),得到將用于測(cè)序的cDNA.在以上過程,將RNA隨機(jī)片段化和采用隨機(jī)引物進(jìn)行反轉(zhuǎn)錄,都是為了使所得cDNA片段較均勻地取自各個(gè)轉(zhuǎn)錄本.為提高測(cè)序效率,一般還需要用電泳切膠法獲取長(zhǎng)度范圍在200bp(±25bp)的cDNA片段,再通過RCR擴(kuò)增,得到最終的cDNA文庫(kù).在上述文庫(kù)制備過程中,如果不是只抽取帶Poly(A)尾的RNA,而是使用全部的RNA,則RNA-seq測(cè)得的就是細(xì)胞中的全部轉(zhuǎn)錄本,如果把帶Poly(A)尾的RNA過濾掉,也可以得到非編碼的RNA轉(zhuǎn)錄本,如果從總RNA中只提取長(zhǎng)度為21~23個(gè)堿基左右的RNA,則得到全部的miRNA(microRNA)轉(zhuǎn)錄本,相應(yīng)的方法也稱作miRNA-seq.樣品制備最終得到的是雙鏈cDNA文庫(kù).在后續(xù)測(cè)序中,測(cè)得的每個(gè)讀段(read)隨機(jī)地來自雙鏈cDNA的某一條鏈,從讀段序列本身無(wú)法得知它是與RNA方向相同還是倒轉(zhuǎn)互補(bǔ),在后續(xù)的讀段定位時(shí)需要兩個(gè)方向都考慮.在新基因識(shí)別等應(yīng)用中,轉(zhuǎn)錄本的方向?qū)蜃⑨層葹橹匾?需要在文庫(kù)制備和測(cè)序中保留RNA的方向信息.最近有文獻(xiàn)報(bào)道了保留方向信息的RNA-seq樣品制備方法.2.2熒光圖像數(shù)據(jù)處理將RNA-seq測(cè)序文庫(kù)加入流動(dòng)槽(flowcell)中的各通道(lane),在橋式PCR擴(kuò)增后,就可以進(jìn)行測(cè)序了.測(cè)序過程中,計(jì)算機(jī)軟件同步地對(duì)熒光圖像數(shù)據(jù)進(jìn)行處理,通過分析熒光信號(hào)來確定被測(cè)堿基,并給出質(zhì)量評(píng)分.按照?qǐng)D像上的位置坐標(biāo),計(jì)算機(jī)程序?qū)⑼晃恢脺y(cè)得的堿基根據(jù)測(cè)序順序連成讀段(read).由于熒光圖像文件所占有的磁盤空間很大,通常GAIIx平臺(tái)一次實(shí)驗(yàn)就能產(chǎn)生上太字節(jié)(TB)的圖像文件,所以一般情況下不予保留原始的熒光圖像數(shù)據(jù),而是只保留程序讀出的讀段數(shù)據(jù)及對(duì)應(yīng)的質(zhì)量分值,這就是多數(shù)實(shí)驗(yàn)室委托測(cè)序中心進(jìn)行RNA-seq測(cè)序后得到的最原始的數(shù)據(jù)..3單個(gè)樣本rna測(cè)序的基本處理流程RNA-seq的基本應(yīng)用是測(cè)量一個(gè)樣本中的基因表達(dá)或轉(zhuǎn)錄組.有實(shí)驗(yàn)表明,新一代高通量測(cè)序技術(shù)重復(fù)數(shù)據(jù)之間的相關(guān)度較高(R2≈0.96),因此,如果對(duì)同一樣本在多個(gè)通道上進(jìn)行了RNA測(cè)序的技術(shù)重復(fù),我們建議可以把幾個(gè)通道的數(shù)據(jù)進(jìn)行合并,這樣等效地增加了測(cè)序深度.本節(jié)討論單個(gè)樣本RNA測(cè)序數(shù)據(jù)的基本處理流程,如圖2a所示.3.1采用bwt的定位算法獲得RNA-seq的原始數(shù)據(jù)后,首先需要將所有測(cè)序讀段通過序列映射(mapping)定位到參考基因組上,這是所有后續(xù)處理和分析的基礎(chǔ).在讀段定位之前,有時(shí)還需要根據(jù)測(cè)序數(shù)據(jù)情況對(duì)其做某些基本的預(yù)處理.例如,過濾掉測(cè)序質(zhì)量較差的讀段、對(duì)miRNA測(cè)序讀段數(shù)據(jù)去除接頭序列等.高通量測(cè)序的海量數(shù)據(jù)對(duì)計(jì)算機(jī)算法的運(yùn)行時(shí)間提出了很高的要求.針對(duì)諸如Illumina/Solexa等測(cè)序平臺(tái)得到的讀段一般較短、且插入刪除錯(cuò)誤較少等特點(diǎn),人們開發(fā)了一些短序列定位算法.這些算法主要采用空位種子索引法(spaced-seedindexing)或Burrows-Wheeler轉(zhuǎn)換(Burrows-WheelerTransform,BWT)技術(shù)來實(shí)現(xiàn).空位種子索引法首先將讀段切分,并選取其中一段或幾段作為種子建立搜索索引,再通過查找索引、延展匹配來實(shí)現(xiàn)讀段定位,通過輪換種子考慮允許出現(xiàn)錯(cuò)配(mismatch)的各種可能的位置組合.BWT方法通過B-W轉(zhuǎn)換將基因組序列按一定規(guī)則壓縮并建立索引,再通過查找和回溯來定位讀段,在查找時(shí)可通過堿基替代來實(shí)現(xiàn)允許的錯(cuò)配.表1列出了目前可免費(fèi)下載使用的部分短序列定位軟件.其中采用空位種子片段索引法的代表是Maq,而采用Burrows-Wheeler轉(zhuǎn)換的代表是Bowtie.總的來說,采用BWT的定位算法在時(shí)間效率上要優(yōu)于空位種子片段索引法.隨著讀長(zhǎng)的增加,允許讀段序列中存在插入刪除(indel)的定位變得可行而重要.由于以上兩類方法對(duì)序列中插入刪除的處理較為困難,近來人們開發(fā)了一些基于改進(jìn)的SmithWaterman動(dòng)態(tài)規(guī)劃算法的序列比對(duì)工具,如BFAST、SHRiMP、Mosaik(/marthlab/Mosaik)等,但算法速度較慢,大多需采用計(jì)算機(jī)并行編程技術(shù)來解決運(yùn)行時(shí)間的問題.在RNA測(cè)序數(shù)據(jù)的基因組定位中,一個(gè)特殊的問題是跨越兩個(gè)外顯子接合區(qū)的讀段(junctionreads)定位.在真核生物中,成熟的mRNA是經(jīng)過由mRNA前體中的外顯子經(jīng)過剪接形成的.如果一個(gè)讀段跨越了兩個(gè)外顯子,那么就無(wú)法將這個(gè)讀段完整地定位到基因組序列上.而同時(shí),這種跨兩個(gè)外顯子的讀段在分析轉(zhuǎn)錄本的剪接形式和研究選擇性剪接中有重要的作用.為了解決這一問題,人們采取兩種典型的策略來進(jìn)行接合區(qū)讀段的定位:一是根據(jù)已知的基因外顯子注釋,構(gòu)建所有可能的外顯子接合區(qū)序列,與基因組序列一并作為定位的參考基因組;二是不依賴基因注釋,而是先利用能完整定位到基因組的讀段得到粗略的外顯子區(qū)域,并結(jié)合剪接位點(diǎn)序列構(gòu)建出可能的剪接位點(diǎn),然后將不能完整定位的讀段分段定位到兩個(gè)外顯子可能的結(jié)合區(qū)域.Illumina/Solexa平臺(tái)提供的RNA-seq軟件分析包GApipeline采用了第一種策略.采用第二種策略的軟件有Tophat和G-Mo.R-Se等,最新的Tophat軟件增加了利用已知外顯子邊界注釋信息的選項(xiàng).不論是哪種測(cè)序平臺(tái),測(cè)序中都不可避免地存在一定的錯(cuò)誤,基因組中又存在單核苷酸多態(tài)性等引起的序列變化,所以在讀段定位時(shí)通常允許一定數(shù)量的錯(cuò)配,可以根據(jù)不同應(yīng)用調(diào)節(jié)允許錯(cuò)配的程度.另一方面,由于基因組中重復(fù)序列和高相似度序列的影響,某些讀段會(huì)出現(xiàn)定位到基因組多個(gè)位置的情況.這些因素影響了各個(gè)讀段到基因組的定位質(zhì)量,在一些新的讀段定位算法中,同時(shí)給出每個(gè)讀段與基因組匹配質(zhì)量.通常在后續(xù)處理前,人們將多定位的讀段都過濾掉,也有人嘗試用適當(dāng)?shù)牟呗园讯喽ㄎ蛔x段“分配”到其中某些位置上.讀段定位到基因組后通常采用SAM(SequenceAlignment/Map)格式或其二進(jìn)制版本BAM格式來存儲(chǔ).二進(jìn)制版本可大大節(jié)省存儲(chǔ)空間,但不能直接用普通文本編輯工具顯示.關(guān)于SAM格式的詳細(xì)介紹,可查閱(/SAM1.pdf).3.2rpm和rpkm在深度測(cè)序技術(shù)出現(xiàn)之前,高通量測(cè)量不同基因表達(dá)水平的主要手段是基因芯片,在此基礎(chǔ)上可以對(duì)不同組織或者不同發(fā)育階段的基因表達(dá)差異和模式進(jìn)行分析.mRNA-seq數(shù)據(jù)最基本的應(yīng)用也是檢測(cè)基因的表達(dá)水平,與基因芯片數(shù)據(jù)相比,RNA測(cè)序得到的是數(shù)字化的表達(dá)信號(hào),具有靈敏度高、分辨率高、無(wú)飽和區(qū)等優(yōu)勢(shì).RNA測(cè)序數(shù)據(jù)是對(duì)提取出的RNA轉(zhuǎn)錄本中隨機(jī)進(jìn)行的短片段測(cè)序,如果一個(gè)轉(zhuǎn)錄本的豐度高,則測(cè)序后定位到其對(duì)應(yīng)的基因組區(qū)域的讀段也就多,可以通過對(duì)定位到基因外顯子區(qū)的讀段計(jì)數(shù)來估計(jì)基因表達(dá)水平.很顯然,讀段計(jì)數(shù)除了與基因真實(shí)表達(dá)水平成正比,還與基因長(zhǎng)度成正比,同時(shí)也與測(cè)序深度即測(cè)序?qū)嶒?yàn)中得到的總讀段數(shù)正相關(guān).為了保持對(duì)不同基因和不同實(shí)驗(yàn)間估計(jì)的基因表達(dá)值的可比性,人們提出了RPM和RPKM的概念.RPM(readspermillionreads)即每百萬(wàn)讀段中來自于某基因的讀段數(shù),考慮了測(cè)序深度對(duì)讀段計(jì)數(shù)的影響.RPKM(readsperkilobasespermillionreads)是每百萬(wàn)讀段中來自于某基因每千堿基長(zhǎng)度的讀段數(shù),公式表示為:RPKM不僅對(duì)測(cè)序深度作了歸一化,而且對(duì)基因長(zhǎng)度也作了歸一化,使得不同長(zhǎng)度的基因在不同測(cè)序深度下得到的基因表達(dá)水平估計(jì)值具有了可比性,是目前最常用的基因表達(dá)估計(jì)方法.軟件rSeq、DEGseq軟件包和Cufflinks等都提供了用上述方法進(jìn)行基因表達(dá)水平計(jì)算的功能.根據(jù)RNA-seq文庫(kù)制備標(biāo)準(zhǔn),在不考慮基因結(jié)構(gòu)的理想情況下,讀段會(huì)均勻地分布在基因上.而實(shí)際上,通過對(duì)實(shí)際數(shù)據(jù)的可視化分析很容易發(fā)現(xiàn),讀段在基因上的分布有著自身的一些模式,呈現(xiàn)出不均勻性(圖3).這一問題已經(jīng)引起很多學(xué)者的關(guān)注.造成讀段分布出現(xiàn)偏好的原因可能有多個(gè)方面:在制備cDNA文庫(kù)時(shí),反轉(zhuǎn)錄所采用的隨機(jī)引物對(duì)RNA序列具有一定的偏好性,使得cDNA片段不能夠完全均勻地取自各轉(zhuǎn)錄本;在PCR擴(kuò)增中,擴(kuò)增效率與序列的GC含量等特征相關(guān),可導(dǎo)致GC含量高的cDNA片段在文庫(kù)中拷貝數(shù)增加超過其他片段;舍棄多定位的讀段也可能導(dǎo)致讀段的非均勻分布;等等.如果能對(duì)讀段分布的不均勻性進(jìn)行建模并加以校正,可以提高RNA-seq推斷基因表達(dá)量的準(zhǔn)確度.但根據(jù)對(duì)實(shí)際數(shù)據(jù)的觀察,對(duì)于較長(zhǎng)轉(zhuǎn)錄本,讀段非均勻分布帶來的誤差很大程度上可相互抵消,用RPKM來估計(jì)基因的表達(dá)水平可以得到比較滿意的結(jié)果.3.3選擇性剪接事件的計(jì)算在真核生物中,選擇性剪接現(xiàn)象普遍存在.基因轉(zhuǎn)錄形成的mRNA前體(pre-mRNA)在剪接過程中因去掉不同的內(nèi)含子區(qū)域或保留不同的外顯子區(qū)域,可形成不同的剪接異構(gòu)體.根據(jù)RNA-seq原理,只要測(cè)序深度足夠深,就能檢測(cè)到所有轉(zhuǎn)錄本的全部序列,包括來自剪接接合區(qū)的序列.利用考慮到接合區(qū)的讀段定位方法,就有可能系統(tǒng)地研究某一組織或某一條件下的基因選擇性剪接事件.前面已經(jīng)提到,Tophat等軟件定位剪接接合區(qū)讀段的策略能標(biāo)定出剪接事件中的兩個(gè)剪接位點(diǎn):供體位點(diǎn)和受體位點(diǎn).通過比較供體位點(diǎn)和受體位點(diǎn)的組合,就能識(shí)別選擇性剪接事件.圖3中包含了選擇性剪接識(shí)別的一個(gè)例子.進(jìn)一步,通過對(duì)供體和受體位點(diǎn)的讀段計(jì)數(shù),結(jié)合外顯子其他區(qū)域的讀段數(shù)據(jù),還能定量地計(jì)算選擇性剪接事件之間的比例.對(duì)于每一個(gè)剪接異構(gòu)體,RNA-seq數(shù)據(jù)能在一定程度上推斷其表達(dá)水平.比如,可以根據(jù)已知外顯子組成和各外顯子長(zhǎng)度對(duì)剪接異構(gòu)體建立數(shù)學(xué)模型,在測(cè)序讀段轉(zhuǎn)錄本上均勻分布的假設(shè)下,利用各外顯子上的讀段數(shù)和接合區(qū)讀段數(shù)求解異構(gòu)體的表達(dá)值.Jiang等的方法及軟件IsoInfer和cufflinks都采用了這種思路來實(shí)現(xiàn)剪接異構(gòu)體的表達(dá)推斷.需要指出的是,某些形式的剪接異構(gòu)體表達(dá)水平在這種方法框架下不可辨識(shí).3.4新基因檢測(cè)在對(duì)RNA-seq數(shù)據(jù)的分析中,人們發(fā)現(xiàn),往往不是所有讀段都能定位到已有注釋的基因區(qū),說明除了轉(zhuǎn)錄噪聲或測(cè)序錯(cuò)誤等的影響外,可能還存在尚未被注釋的基因.這里,我們把這種尚未注釋的基因稱為新基因,包括新的蛋白質(zhì)編碼基因和非編碼RNA基因.能檢測(cè)新基因,尤其是低表達(dá)基因是RNA-seq技術(shù)優(yōu)于基因芯片的特點(diǎn)之一,因?yàn)樗恍枰靡阎蜃⑨寔碓O(shè)計(jì)檢測(cè)探針.RNA-seq技術(shù)靈敏度高,但樣品污染、測(cè)序錯(cuò)誤等仍可能帶來背景噪聲.從基因組未注釋區(qū)域的RNA測(cè)序讀段信號(hào)中檢測(cè)新基因是典型的信號(hào)檢測(cè)問題.如何控制新基因識(shí)別的誤發(fā)現(xiàn)率(FDR)是檢測(cè)方法的關(guān)鍵.Useq軟件包將ChIP-seq數(shù)據(jù)分析的方法移植到RNA-seq數(shù)據(jù)上,用滑窗的方法來識(shí)別測(cè)序讀段定位富集的區(qū)域,給出反映滑窗所在區(qū)域讀段富集顯著程度的P值(P-value)及新基因誤發(fā)現(xiàn)率,通過設(shè)定P值或誤發(fā)現(xiàn)率的閾值,可篩選出讀段富集的區(qū)域,再將相鄰區(qū)域合并或根據(jù)剪接接合區(qū)讀段將相應(yīng)區(qū)域連接,完成新基因的檢測(cè).3.5轉(zhuǎn)錄水平及測(cè)序效率對(duì)于復(fù)雜的組學(xué)數(shù)據(jù),能盡可能方便地直接觀察數(shù)據(jù)對(duì)于數(shù)據(jù)的分析和解釋都非常重要,對(duì)新一代測(cè)序數(shù)據(jù)的可視化和交互展示是一個(gè)非常重要但容易被人忽視的問題.不深入考查數(shù)據(jù)的細(xì)節(jié),而是滿足于對(duì)數(shù)據(jù)的統(tǒng)計(jì)分析,是高通量數(shù)據(jù)應(yīng)用中經(jīng)常容易陷入的誤區(qū),方便有效的可視化工具能夠幫助避免這樣的誤區(qū).表2列出了部分適用于RNA-seq數(shù)據(jù)的全基因組瀏覽器,其中比較具有代表性的有UCSCGenomeBrowser、CisGenomeBrowser和IGV(IntegrativeGenomicsViewer)等.這些瀏覽器具有如下特點(diǎn):a.能在不同尺度下顯現(xiàn)單個(gè)或多個(gè)讀段在基因組上的位置,包括來源于剪接接合區(qū)的讀段;b.能在不同尺度下顯示不同區(qū)域的讀段豐度,以反映不同區(qū)域的轉(zhuǎn)錄水平或測(cè)序效率;c.能顯示基因及其剪接異構(gòu)體的注釋信息;d.能顯示其他注釋信息,例如物種間基因組序列保守性、序列GC含量等;e.能直接或間接支持SAM/BAM讀段定位數(shù)據(jù)存儲(chǔ)格式.UCSCGenomeBrowser屬于基于網(wǎng)絡(luò)模式的全基因組瀏覽器,所有數(shù)據(jù)都需要上載到遠(yuǎn)程服務(wù)器,經(jīng)過處理后將圖形返回客戶端顯示.圖3中的例子就是從UCSCGenomeBrowser的顯示截取的.CisGenomeBrowser除對(duì)讀段的可視化外,用描述統(tǒng)計(jì)學(xué)方法對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分類別統(tǒng)計(jì)也十分重要.例如,統(tǒng)計(jì)讀段在各個(gè)染色體上的分布情況和在注釋的外顯子、內(nèi)含子、剪接接合區(qū)、基因間區(qū)的分布情況等.目前,已經(jīng)有一些用于測(cè)序數(shù)據(jù)注釋的生物信息學(xué)軟件,比如SAMtools、BEDtools等,但由于測(cè)序技術(shù)發(fā)展迅速,用戶需求因人而異,用戶經(jīng)常還需要根據(jù)需求編寫一定的程序或腳本完成或完善注釋分析的任務(wù).對(duì)于熟悉圖形用戶界面的研究人員,還可以利用UCSCTableBrowse和Galaxy來配合完成注釋分析.由于UCSCTableBrowser集成了大量基因組尺度上的注釋信息,而Galaxy又為用戶提供了書寫簡(jiǎn)單、接口明晰和直觀的數(shù)據(jù)處理流程,這種方法十分方便有效,也為很多學(xué)者在展示研究成果時(shí)所采用.以上描述了對(duì)于基因組已知的物種進(jìn)行RNA-seq數(shù)據(jù)處理基本流程,圖2a給出了其主要步驟.若研究對(duì)象尚未完成基因組測(cè)序,則需要采用讀段的從頭拼裝(denovoassembly)來代替讀段定位,后續(xù)流程也須做相應(yīng)的調(diào)整.若RNA-seq實(shí)驗(yàn)在文庫(kù)制備時(shí)保留了RNA的方向信息,則應(yīng)分別研究來自正鏈和反鏈的轉(zhuǎn)錄產(chǎn)物,并通過與基因注釋比較來檢測(cè)反義轉(zhuǎn)錄本.最近,RNA-MATE軟件在其分析流程中加入如此的處理策略.此外,通過分析定位到外顯子接合區(qū)的讀段,還可以獲取轉(zhuǎn)錄本結(jié)構(gòu),這為研究基因的剪接調(diào)控機(jī)理提供了重要信息.而利用RNA-seq數(shù)據(jù)提供的序列信息,通過與DNA序列的細(xì)致比較可分析轉(zhuǎn)錄組的序列差異(如SNP等),從而研究等位基因的表達(dá)模式及RNA編輯等.最后需要指出,由于miRNA在序列和結(jié)構(gòu)上具有一定的特點(diǎn),miRNA-seq數(shù)據(jù)的基本處理流程也與本節(jié)所述有所不同,感興趣的讀者可參考軟件工具miRDeep提供的處理策略.4剪接異構(gòu)體的差異檢驗(yàn)很多RNA-seq實(shí)驗(yàn)的目的是為了比較兩種或多種樣本中基因表達(dá)或整個(gè)轉(zhuǎn)錄組的差異,如比較癌癥組織和正常組織的轉(zhuǎn)錄組差異等.這些差異既包括通常意義下的差異表達(dá)基因,也主要包括選擇性剪接模式的差異、剪接異構(gòu)體表達(dá)的差異、非編碼轉(zhuǎn)錄本的差異等.這些差異一般可以用一些統(tǒng)計(jì)假設(shè)檢驗(yàn)方法檢測(cè),但這種檢驗(yàn)有時(shí)會(huì)受到測(cè)序深度、基因長(zhǎng)度等因素的影響,需要對(duì)結(jié)果進(jìn)行仔細(xì)分析,消除可能的混雜因素,必要時(shí)可以用讀段的絕對(duì)表達(dá)值倍數(shù)變化(fold-change)來作為補(bǔ)充.圖2b給出了兩類樣本數(shù)據(jù)分析的框架.4.1差異表達(dá)基因檢測(cè)雖然新一代測(cè)序相對(duì)第一代測(cè)序的單位成本大大降低,但是,利用RNA測(cè)序進(jìn)行基因表達(dá)研究的成本仍很高,因此,很多實(shí)驗(yàn)室沒有條件進(jìn)行樣本重復(fù).如果兩類樣本均沒有生物重復(fù),例如只對(duì)兩個(gè)細(xì)胞系各進(jìn)行一次mRNA樣本測(cè)序,則可以用隨機(jī)采樣模型通過假設(shè)檢驗(yàn)來分析差異表達(dá).對(duì)于某個(gè)基因,如果一個(gè)讀段來自于這個(gè)基因,我們稱事件A發(fā)生.對(duì)于一次RNA-seq實(shí)驗(yàn),事件A發(fā)生的概率可以用這個(gè)基因上的讀段數(shù)n除以所有基因上的讀段總數(shù)N來估計(jì),即RPM.事件A發(fā)生的概率反應(yīng)了這個(gè)基因的表達(dá)水平.如果要判斷某個(gè)基因在兩個(gè)樣本中的表達(dá)水平是否一致,就可以通過檢驗(yàn)事件A在兩種條件下發(fā)生的概率是否一致來實(shí)現(xiàn),采用似然比檢驗(yàn)、Fisher精確檢驗(yàn)以及基于MA圖的統(tǒng)計(jì)檢驗(yàn)方法等.同樣,也可用RPKM作為統(tǒng)計(jì)量來進(jìn)行假設(shè)檢驗(yàn)分析,由于是比較同一個(gè)基因在兩個(gè)樣本間的差異,基因長(zhǎng)度的影響被抵消,用RPKM和用RPM得到的結(jié)果相似.對(duì)無(wú)生物重復(fù)的RNA-seq數(shù)據(jù)進(jìn)行差異表達(dá)基因分析,已經(jīng)有幾個(gè)公開發(fā)表的軟件,包括DEGseq、Useq、Cufflinks中的Cuffdiff模塊等.圖4展示了我們開發(fā)的DEGseq軟件提供的多種差異表達(dá)基因識(shí)別方法的應(yīng)用例子.如果每一類樣本都包含了若干生物重復(fù),如病人和正常人對(duì)照研究,則可以沿用基因芯片數(shù)據(jù)分析中的很多方法.比如,可以用t檢驗(yàn)結(jié)合倍數(shù)變化的方法來分析差異表達(dá).如果兩類樣本具有配對(duì)的信息,也可以通過整合每對(duì)樣本分析結(jié)果來實(shí)現(xiàn).其步驟為,先在每對(duì)樣本中識(shí)別出差異表達(dá)的基因,再尋找這若干組差異表達(dá)基因之間的相同者,或用投票的方法來為基因的差異程度打分.針對(duì)某些RNA-seq數(shù)據(jù)生物樣本量小,R軟件包DEGseq和edgeR等還專門提供了基于改進(jìn)模型的統(tǒng)計(jì)方法.此外,一類將分類器與特征選擇包裹在一起的方法也同樣適用于此類問題(見4.3).4.2差異表達(dá)剪接異構(gòu)體的識(shí)別差異表達(dá)剪接異構(gòu)體的識(shí)別方法與差異表達(dá)基因的識(shí)別相似.如果把剪接異構(gòu)體看成是獨(dú)立的基因,那么前面討論的用于識(shí)別差異表達(dá)基因的方法對(duì)剪接異構(gòu)體完全適用.但是,注意到來自于同一個(gè)基因的剪接異構(gòu)體并不獨(dú)立,某些假設(shè)檢驗(yàn)的基本條件并不滿足,得到的結(jié)果就不一定正確.此外,由于現(xiàn)在剪接異構(gòu)體表達(dá)推斷的方法還不夠成熟,加之在基因結(jié)構(gòu)不可辨識(shí)的剪接異構(gòu)體上作表達(dá)推斷會(huì)出現(xiàn)病態(tài)結(jié)果,差異表達(dá)剪接異構(gòu)體的識(shí)別問題還處于探索的階段.目前,在剪接異構(gòu)體表達(dá)水平可辨識(shí)且讀段覆蓋度較高的基因上,BASIS方法通過貝葉斯模型來推斷差異表達(dá)的剪接異構(gòu)體.換一個(gè)角度,剪接異構(gòu)體由選擇性剪接造成,如果剪接異構(gòu)體的表達(dá)有差異,那么導(dǎo)致這些異構(gòu)體的選擇性剪接事件及異構(gòu)體特異的外顯子的表達(dá)也會(huì)有差異.因此,對(duì)差異表達(dá)剪接異構(gòu)體的識(shí)別可以轉(zhuǎn)變?yōu)榉治鲞x擇性剪接事件和外顯子表達(dá)的差異.外顯子表達(dá)差異的分析可以完全利用基因表達(dá)差異的分析方法.而剪接接合區(qū)也可以看成是一個(gè)較短的“外顯子”(長(zhǎng)度一般與測(cè)序長(zhǎng)度相當(dāng)).不過,由于外顯子長(zhǎng)度較基因的長(zhǎng)度短,對(duì)應(yīng)的讀段數(shù)量較少,差異識(shí)別的敏感度會(huì)有所下降.Solas方法就是根據(jù)類似的原理,采用統(tǒng)計(jì)學(xué)假設(shè)檢驗(yàn)的方法來識(shí)別差異表達(dá)的剪接異構(gòu)體.4.3核心評(píng)估方法通過統(tǒng)計(jì)方法識(shí)別出來的差異表達(dá)基因及剪接異構(gòu)體能否有效地區(qū)別兩類樣本,可以通過分類分析進(jìn)一步證實(shí).如果把每個(gè)基因(或剪接異構(gòu)體)的表達(dá)值作為特征,則差異表達(dá)基因(或剪接異構(gòu)體)的選取也就是特征篩選的過程.把前面用統(tǒng)計(jì)方法等檢測(cè)出來的差異表達(dá)基因(或剪接異構(gòu)體)用于分類分析,常被稱為過濾法.另一類基于分類器的包裹法,例如R-SVM、SVM-RFE等,可以根據(jù)每個(gè)特征在分類器中所占的權(quán)重來篩選特征,因此也可以用于差異表達(dá)基因(或剪接異構(gòu)體)的識(shí)別.分類的性能可以用交叉驗(yàn)證(cross-validation,CV)方法來評(píng)估.需要特別注意的是,交叉驗(yàn)證應(yīng)該包括對(duì)特征選擇步驟的交叉驗(yàn)證,防止發(fā)生信息泄露而導(dǎo)致評(píng)估結(jié)果過于樂觀.具體做法是:將樣本按一定的策略分成兩份,一份(通常是樣本數(shù)多的一份)用于特征選取和分類器訓(xùn)練,而用余下的樣本進(jìn)行分類器性能的估計(jì);重復(fù)以上步驟多次,就得到交叉驗(yàn)證錯(cuò)誤率.必要時(shí)還可以用隨機(jī)置換檢驗(yàn)(permutationtest)來推斷所得錯(cuò)誤率的統(tǒng)計(jì)顯著性.當(dāng)樣本數(shù)較小時(shí),可以采用留一法交叉驗(yàn)證(leave-one-outcross-validation,LOOCV).4.4goseq在生物數(shù)學(xué)過程分析中的應(yīng)用檢測(cè)差異表達(dá)的基因或差異表達(dá)異構(gòu)體是人們認(rèn)識(shí)所研究的生物問題機(jī)理的第一步,接下來需要從功能上研究這些差異轉(zhuǎn)錄現(xiàn)象的分子機(jī)理.這與在基因芯片應(yīng)用中所面臨的是同樣的生物學(xué)問題,對(duì)芯片數(shù)據(jù)分析結(jié)果的后續(xù)處理方法,都可以借鑒到測(cè)序數(shù)據(jù)上來.如何進(jìn)一步地從機(jī)理來解釋結(jié)果,還需結(jié)合已知生物學(xué)知識(shí)進(jìn)行后續(xù)分析.人們對(duì)基因芯片得到的基因表達(dá)數(shù)據(jù)進(jìn)行分析的很多方法都可以用到RNA-seq數(shù)據(jù)上來,比如利用機(jī)器學(xué)習(xí)方法進(jìn)行分類和特征選擇,對(duì)差異表達(dá)的基因進(jìn)行GO(geneontology)類別富集分析、信號(hào)通路富集分析等,一些常用的分析工具包括GoMiner、DAVID和VisNT等.需要說明的是,在各種以差異表達(dá)基因?yàn)榛A(chǔ)的分析中,由于基因表達(dá)水平都是通過讀段計(jì)數(shù)來估計(jì)的,表達(dá)水平較高或轉(zhuǎn)錄本較長(zhǎng)的基因擁有更多的讀段,更容易被多數(shù)統(tǒng)計(jì)方法識(shí)別為差異表達(dá)基因.這種偏好可能對(duì)后續(xù)分析帶來影響.以GO類別富集分析為例,這種偏好將導(dǎo)致長(zhǎng)基因占主導(dǎo)的功能類別更有可能被識(shí)別為富集的功能.這將對(duì)生物機(jī)理的研究帶來誤導(dǎo).最近,Young等發(fā)展了一種GOseq方法,針對(duì)這一偏好對(duì)GO類別富集分析做了改進(jìn).5實(shí)驗(yàn)數(shù)據(jù)處理高通量測(cè)序技術(shù)的發(fā)展十分迅速,這要求相應(yīng)的數(shù)據(jù)處理與分析方法快速跟進(jìn).正是這些方法,架起了高通量實(shí)驗(yàn)數(shù)據(jù)與科學(xué)問題之間的橋梁.這種橋梁作用正日趨重要,也為生物信息學(xué)帶來了挑戰(zhàn).這里,我們重點(diǎn)討論兩方面的挑戰(zhàn):a.如何實(shí)現(xiàn)剪接接合區(qū)讀段的準(zhǔn)確定位?b.在數(shù)據(jù)處理各階段中,如何對(duì)RNA-seq數(shù)據(jù)的系統(tǒng)誤差和固有偏好建?；蜓a(bǔ)償,以消除它們可能帶來的錯(cuò)誤推斷及結(jié)論?5.1剪接接合區(qū)讀段映射定位算法的突破測(cè)序技術(shù)的一個(gè)發(fā)展趨勢(shì)是測(cè)序長(zhǎng)度不斷增加.隨著讀長(zhǎng)的增加,RNA-seq中來自剪接接合區(qū)的讀段會(huì)越來越多.我們粗略估算,按照人類基因組refSeq基因注釋,一般情況下,如果測(cè)序讀長(zhǎng)為50個(gè)堿基,則約有10%的讀段來自剪接接合區(qū).而當(dāng)測(cè)序長(zhǎng)度達(dá)到100個(gè)堿基時(shí),這個(gè)比例將達(dá)到25%左右.對(duì)這些剪接接合區(qū)讀段的分析,將使我們能夠更準(zhǔn)確地檢測(cè)剪接事件和推斷剪接異構(gòu)體的表達(dá)水平,大大推進(jìn)人們對(duì)選擇性剪接的研究.在RNA-seq出現(xiàn)的早期,人們沒有意識(shí)到剪接接合區(qū)讀段的重要性.因?yàn)楫?dāng)時(shí)的讀長(zhǎng)只有20~30個(gè)堿基,來自剪接接合區(qū)的讀段所占比例甚小.當(dāng)時(shí)讀段定位的通常做法是,先將讀段與全基因組序列做映射定位,再考慮不能定位的讀段是否來自于剪接接合區(qū).這種做法雖然在一定程度上保證了讀段定位的比率,但由于基因組中重復(fù)序列和相似序列的存在,部分接合區(qū)讀段有可能在容許錯(cuò)配的情況下被定位到基因組上其他位置,從而失去了定位到正確的剪接接合區(qū)的機(jī)會(huì).在讀段定位時(shí),如果要同時(shí)考慮基因組序列和剪接接合區(qū)序列,就要利用已知的剪接事件注釋,這是目前軟件通用的方法.然而,包括人類在內(nèi)的各物種的基因注釋信息都還有待完善,也沒有較完