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

1、本文檔下載自HYPERLINK /文檔下載網(wǎng)，內(nèi)容可能不完整，您可以點(diǎn)擊以下網(wǎng)址繼續(xù)閱讀或下載：HYPERLINK /doc/be6b0724c324a6024ee96d8b/doc/be6b0724c324a6024ee96d8b新一代高通量 RNA 測(cè)序數(shù)據(jù)的處理與分析 _新一代高通量 RNA 測(cè)序數(shù)據(jù)的處理與分析 _生物化學(xué)與生物物理進(jìn)展,37(8):834846此處圖片未下載成功此處圖片未下載成功新一代高通量RNA測(cè)序數(shù)據(jù)的處理與分析*王曦1)汪小我1)王立坤1,2)馮智星1)張學(xué)工1)*)(1)生物信息學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，清華信息科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室(籌)生物信息學(xué)研究部，清華大學(xué)自

2、動(dòng)化系，北京100084；吉林大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春130012)摘要隨著新一代高通量DNA測(cè)序技術(shù)的快速發(fā)展，RNA測(cè)序(RNA-seq)已成為基因表達(dá)和轉(zhuǎn)錄組分析新的重要手段seq技術(shù)產(chǎn)生的海量數(shù)據(jù)為生物信息學(xué)帶來(lái)了新的機(jī)遇和挑戰(zhàn)有效地對(duì)測(cè)序數(shù)據(jù)進(jìn)行針對(duì)性的生物信息學(xué)處理和分析，成為RNA-seq技術(shù)能否在科學(xué)探索中發(fā)揮重大作用的關(guān)鍵以新一代Illumina/Solexa測(cè)序平臺(tái)所產(chǎn)生的數(shù)據(jù)為例，在扼要介紹高通量RNA-seq測(cè)序流程的基礎(chǔ)上，對(duì)RNA-seq數(shù)據(jù)處理和分析的方法和現(xiàn)有軟件做一個(gè)較為全面的綜述，并對(duì)其中有待進(jìn)一步研究的問(wèn)題進(jìn)行展望關(guān)鍵詞高通量RNA測(cè)序，轉(zhuǎn)錄組，基因

3、表達(dá)，數(shù)據(jù)處理與分析，生物信息學(xué)，Q6，Q7:10.3724/SP.J.1206.2010.00151學(xué)科分類(lèi)號(hào)近年來(lái)，新一代高通量測(cè)序技術(shù)得到了突飛猛進(jìn)的發(fā)展，在此基礎(chǔ)上，高通量RNA測(cè)序即RNA-seq1-5也迅速發(fā)展與基因芯片技術(shù)相比，RNA-seq無(wú)需設(shè)計(jì)探針，能在全基因組范圍內(nèi)以單堿基分辨率檢測(cè)和量化轉(zhuǎn)錄片段，并能應(yīng)用于基因組圖譜尚未完成的物種6，具有信噪比高、分辨率高、應(yīng)用范圍廣等優(yōu)勢(shì)，正成為研究基因表達(dá)和轉(zhuǎn)錄組的重要實(shí)驗(yàn)手段seq為基因組學(xué)的研究帶來(lái)了高分辨率的海量數(shù)據(jù)，如何有效處理和分析這些海量數(shù)據(jù)成為這一新技術(shù)能否帶來(lái)新的科學(xué)發(fā)現(xiàn)的關(guān)鍵，一些/doc/be6b0724c324

4、a6024ee96d8b生物信息學(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)生過(guò)程及數(shù)據(jù)處理和分析的基本流程、關(guān)鍵方法和現(xiàn)有軟件進(jìn)行較全面的介紹，并討論RNA-seq數(shù)據(jù)分析中存在的挑戰(zhàn)格更便宜、自動(dòng)化程度更高的測(cè)序技術(shù)自2005年以來(lái)，以Roche公司的454技術(shù)、Illumina公司的Solexa技術(shù)和ABI公司的SOLiD技術(shù)為標(biāo)志的新一代測(cè)序技術(shù)相繼誕生8新一代測(cè)序技術(shù)又稱(chēng)作深度測(cè)序技術(shù)，主要特點(diǎn)是測(cè)序通量高、測(cè)序時(shí)間和成本顯著下降9把這種高通量測(cè)序技術(shù)應(yīng)用到由mRNA逆

5、轉(zhuǎn)錄生成的cDNA上，從而獲得來(lái)自不同基因的mRNA片段在特定樣本中的含量，這就是mRNA測(cè)序或mRNA-seq同樣原理，各種類(lèi)型的轉(zhuǎn)錄本都可以用深度測(cè)序技術(shù)進(jìn)行高通量定量檢測(cè)，統(tǒng)稱(chēng)作RNA-seq或RNA測(cè)序目前，在已經(jīng)推出的幾種新一代測(cè)序平臺(tái)中，Illumina/Solexa測(cè)序平臺(tái)上的RNA-seq應(yīng)用最廣，我們以此為例來(lái)綜述RNA-seq數(shù)據(jù)處理和分析的生物信息學(xué)問(wèn)題和方法.高通量測(cè)序技術(shù)簡(jiǎn)介*國(guó)家自然科學(xué)基金資助項(xiàng)目(60702002,60721003,30873464,60905013)和東南大學(xué)生物電子學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放研究基金資助項(xiàng)目.*通訊聯(lián)系人./p>

6、,E-mail:zhangxg收稿日期：2010-03-25，接受日期：2010-04-30誕生于20世紀(jì)70年代的Sanger法是最早被廣泛應(yīng)用的DNA測(cè)序技術(shù)7，也是完成人類(lèi)基因組計(jì)劃的基礎(chǔ)但是，由于它測(cè)序通量低，費(fèi)時(shí)費(fèi)力，科學(xué)家們一直在尋求通量更高、速度更快、價(jià);37(8)王曦等：新一代高通量RNA測(cè)序數(shù)據(jù)的處理與分析835/Solexa測(cè)序技術(shù)的基本原理是邊合成邊測(cè)序(sequencingb/doc/be6b0724c324a6024ee96d8bysynthesis，SBS)10-12，即測(cè)序過(guò)程是以DNA單鏈為模板，在生成互補(bǔ)鏈時(shí)，利用帶熒光標(biāo)記的dNTP發(fā)出不同顏色的熒光來(lái)確定不

7、同的堿基新加入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)13(/)初期的Illumina/Solexa測(cè)序技術(shù)只能在較短的測(cè)序讀長(zhǎng)上(2030堿基)保證較高的正確率隨著技術(shù)的改進(jìn)，目前的讀長(zhǎng)已經(jīng)增加到100堿基以上同時(shí)，隨著雙端測(cè)序(paired-end，PE)技術(shù)的成熟，測(cè)序長(zhǎng)度更可達(dá)到單端測(cè)序的2倍，測(cè)序通量也隨之增加這種測(cè)序技術(shù)是Solexa公司發(fā)展起來(lái)的，2007年被Il

8、lumina公司收購(gòu)，因此現(xiàn)在通常被稱(chēng)為Illumina/Solexa測(cè)序技術(shù)近兩年來(lái)，Illumina/Solexa測(cè)序平臺(tái)不斷升級(jí)，相繼推出了GA(GenomeAnalyzer)、GAIIx、HiSeq2000等測(cè)序儀更多關(guān)于高通量測(cè)序平臺(tái)的介紹，可以查閱相關(guān)文獻(xiàn)9,14-16RNA鄄seq測(cè)序文庫(kù)制備和測(cè)序平臺(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)單的介紹援1RNA鄄seq測(cè)序文庫(kù)制備對(duì)于mRNA-seq實(shí)驗(yàn)，從總RNA到最終的cDNA文庫(kù)制備完成主要包括以下步驟首先，用Poly(T)寡聚核苷酸從總RNA中抽取全部帶P

9、oly(A)尾的RNA，其中的主要部分就是編碼基因所轉(zhuǎn)錄的mRNA將所得RNA隨機(jī)打斷成片段，再用隨機(jī)引物和逆轉(zhuǎn)錄酶從RNA片段合成cDNA片段然后，對(duì)cDNA片段進(jìn)行末端修復(fù)并連接測(cè)序接頭(adapter)，得到將用于測(cè)序的cDNA在以上過(guò)程，將RNA隨機(jī)片段化和采用隨機(jī)引物進(jìn)行反轉(zhuǎn)錄，都是為了使所得cDNA片段較均勻地取自各個(gè)轉(zhuǎn)錄本為提高測(cè)序效率，一般還需要用電泳切膠/doc/be6b0724c324a6024ee96d8b法獲取長(zhǎng)度范圍在200bp(依25bp)的cDNA片段，再通過(guò)RCR擴(kuò)增，得到最終的cDNA文庫(kù)在上述文庫(kù)制備過(guò)程中，如果不是只抽取帶Poly(A)尾的RNA，而是使用

10、全部的RNA，則RNA-seq測(cè)得的就是細(xì)胞中的全部轉(zhuǎn)錄本，如果把帶Poly(A)尾的RNA過(guò)濾掉，也可以得到非編碼的RNA轉(zhuǎn)錄本，如果從總RNA中只提取長(zhǎng)度為2123個(gè)堿基左右的RNA，則得到全部的miRNA(microRNA)轉(zhuǎn)錄本,相應(yīng)的方法也稱(chēng)作miRNA-seq.樣品制備最終得到的是雙鏈cDNA文庫(kù)在后續(xù)測(cè)序中，測(cè)得的每個(gè)讀段(read)隨機(jī)地來(lái)自雙鏈cDNA的某一條鏈，從讀段序列本身無(wú)法得知它是與RNA方向相同還是倒轉(zhuǎn)互補(bǔ)，在后續(xù)的讀段定位時(shí)需要兩個(gè)方向都考慮在新基因識(shí)別等應(yīng)用中，轉(zhuǎn)錄本的方向?qū)蜃⑨層葹橹匾?，需要在文?kù)制備和測(cè)序中保留RNA的方向信息最近有文獻(xiàn)報(bào)道了保留方向信息

11、的RNA-seq樣品制備方法17-20援2測(cè)序平臺(tái)數(shù)據(jù)輸出將RNA-seq測(cè)序文庫(kù)加入流動(dòng)槽(flowcell)中的各通道(lane)，在橋式PCR擴(kuò)增后，就可以進(jìn)行測(cè)序了測(cè)序過(guò)程中，計(jì)算機(jī)軟件同步地對(duì)熒光圖像數(shù)據(jù)進(jìn)行處理，通過(guò)分析熒光信號(hào)來(lái)確定被測(cè)堿基，并給出質(zhì)量評(píng)分按照?qǐng)D像上的位置坐標(biāo)，計(jì)算機(jī)程序?qū)⑼晃恢脺y(cè)得的堿基根據(jù)測(cè)序順序連成讀段(read)由于熒光圖像文件所占有的磁盤(pán)空間很大，通常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è)序后得到的

12、最原始的數(shù)據(jù)為了便于測(cè)序數(shù)據(jù)的發(fā)布和共享，高通量測(cè)序數(shù)據(jù)以FASTQ格式來(lái)記錄所測(cè)的堿基讀段和質(zhì)量分?jǐn)?shù)如圖1所示，F(xiàn)ASTQ格式以測(cè)序讀段為單位存儲(chǔ)，每條讀段占4行，其中第1行和第3行由文件識(shí)別標(biāo)志和讀段名(ID)組成(第1行以“”開(kāi)頭而第3行以“ ”開(kāi)頭；第3行中ID可以省略，但“ ”不能省略)，第2行為堿基序列，第4行為對(duì)應(yīng)的測(cè)序質(zhì)量分?jǐn)?shù)關(guān)于FASTQ格式更多地介紹可參考文獻(xiàn)21為方便保存和共享各實(shí)驗(yàn)室產(chǎn)生/doc/be6b0724c324a6024ee96d8b的高通量測(cè)序數(shù)據(jù)，NCBI、EBI、DDBJ等數(shù)據(jù)中心建立了大容量的數(shù)據(jù)庫(kù)SRA(SequenceReadArchive,/Tr

13、aces/sra)來(lái)存放共享的測(cè)序數(shù)據(jù)22-23836生物化學(xué)與生物物理進(jìn)展Prog.Biochem.Biophys.2010;37(8)每4行標(biāo)識(shí)為一個(gè)測(cè)序讀段讀段識(shí)別碼堿基序列此處圖片未下載成功此處圖片未下載成功此處圖片未下載成功 讀段識(shí)別碼測(cè)序質(zhì)量分?jǐn)?shù)鄄seq數(shù)據(jù)的基本處理seq的基本應(yīng)用是測(cè)量一個(gè)樣本中的基因表達(dá)或轉(zhuǎn)錄組有實(shí)驗(yàn)表明，新一代高通量測(cè)序技術(shù)重復(fù)數(shù)據(jù)之間的相關(guān)度較高(R2抑0.96)1-2，因此，如果對(duì)同一樣本在多個(gè)通道上進(jìn)行了RNA測(cè)序的技術(shù)重復(fù)，我們建議可以把幾個(gè)通道的數(shù)據(jù)進(jìn)行合并，這樣等效地增加了測(cè)序深度本節(jié)討論單個(gè)樣本RNA測(cè)序數(shù)據(jù)的基本處理流程，如圖2a所示圖1讀段

14、FASTQ數(shù)據(jù)格式示例Fig.1FASTQformatexamples(a)(b)圖2所示的流程,虛線(xiàn)箭頭表示可選輸入.鄄seq數(shù)據(jù)處理和分析流程圖鄄seqdataprocessingandanalysis(a)RNA-seq數(shù)據(jù)的基本處理,其方法介紹見(jiàn)正文第3節(jié).(b)兩類(lèi)樣本RNA-seq數(shù)據(jù)比較分析的框架,對(duì)應(yīng)于正文的第4節(jié).(b)中虛線(xiàn)框內(nèi)為(a)援1讀段定位獲得RNA-seq的原始數(shù)據(jù)后，首先需要將所有測(cè)序讀段通過(guò)序列映射(mapping)定位到參考基因組上，這是所有后續(xù)處理和分析的基礎(chǔ)在讀段定位之前，有時(shí)還需要根據(jù)測(cè)序數(shù)據(jù)情況對(duì)其做某些基本的預(yù)處理例如，過(guò)濾掉測(cè)序質(zhì)量較差的讀/do

15、c/be6b0724c324a6024ee96d8b段、對(duì)miRNA測(cè)序讀段數(shù)據(jù)去除接頭序列等高通量測(cè)序的海量數(shù)據(jù)對(duì)計(jì)算機(jī)算法的運(yùn)行時(shí)間提出了很高的要求針對(duì)諸如Illumina/Solexa等測(cè)序平臺(tái)得到的讀段一般較短、且插入刪除錯(cuò)誤較;37(8)王曦等：新一代高通量RNA測(cè)序數(shù)據(jù)的處理與分析837少等特點(diǎn)，人們開(kāi)發(fā)了一些短序列定位算法這些算法主要采用空位種子索引法(spaced-seedindexing)或Burrows-Wheeler轉(zhuǎn)換(Burrows-WheelerTransform，BWT)技術(shù)來(lái)實(shí)現(xiàn)24空位種子索引法首先將讀段切分，并選取其中一段或幾段作為種子建立搜索索引，再通過(guò)查

16、找索引、延展匹配來(lái)實(shí)現(xiàn)讀段定位，通過(guò)輪換種子考慮允許出現(xiàn)錯(cuò)配(mismatch)的各種可能的位置組合BWT方法通過(guò)B-W轉(zhuǎn)換25將基因組序列按一定規(guī)則壓縮并建立索引，再通過(guò)查找和回溯來(lái)定位讀段，在查找時(shí)可通過(guò)堿基替代來(lái)實(shí)現(xiàn)允許的錯(cuò)配表1列出了目前可免費(fèi)下載使用的部分短序列定位軟件其中采用空位種子片段索引法的代表是Maq26，而采用Burrows-Wheeler轉(zhuǎn)換的代表是Bowtie27總的來(lái)說(shuō)，采用BWT的定位算法在時(shí)間效率上要優(yōu)于空位種子片段索引法24,28隨著讀長(zhǎng)的增加，允許讀段序列中存在插入刪除(indel)的定位變得可行而重要由于以上兩類(lèi)方法對(duì)序列中插入刪除的處理較為困難，近來(lái)人們開(kāi)

17、發(fā)了一些基于改進(jìn)的Smith-Waterman動(dòng)態(tài)規(guī)劃算法29的序列比對(duì)工具，如BFAST30、SHRiMP31、Mosaik(/marthlab/Mosaik)等，但算法速度較慢，大多需采用計(jì)算機(jī)并行編程技術(shù)來(lái)解決運(yùn)行時(shí)間的問(wèn)題表1Table1名稱(chēng)MAQ26Bowtie27BWA32ZOOM33ELANDSOAP234RazerS35NovoalignSHRiMP31BFAST30Mosaik)否是是否否/doc/be6b0724c324a6024ee96d8b否否是否質(zhì)量2)是是是否否否否是是適用于Illumina/Solexa測(cè)序平臺(tái)的讀段定位軟件Mappers/alignersforI

18、llumina/Solexasequencingdata主要采用技術(shù)空位種子BWTBWT空位種子空位種子BWTq-grams過(guò)濾Needleman-Wunsch算法空位種子q-grams過(guò)濾Smith-Waterman算法:/shrimp:/index.php/BFAST/marthlab/Mosaik網(wǎng)址:/:/index.shtml/bwa.shtml:/products/zoom:/software/genome_analyzer_software.ilmn:/www.seqan.de/projects/razers.html是是是是Waterman算法并行編程Smith-Waterma

19、n算法并行編程):是否能以SAM格式輸出;2)質(zhì)量:是否提供讀段定位質(zhì)量信息;BWT:Burrows-Wheeler轉(zhuǎn)換.在RNA測(cè)序數(shù)據(jù)的基因組定位中，一個(gè)特殊的問(wèn)題是跨越兩個(gè)外顯子接合區(qū)的讀段(junctionreads)定位在真/doc/be6b0724c324a6024ee96d8b核生物中，成熟的mRNA是經(jīng)過(guò)由mRNA前體中的外顯子經(jīng)過(guò)剪接形成的如果一個(gè)讀段跨越了兩個(gè)外顯子，那么就無(wú)法將這個(gè)讀段完整地定位到基因組序列上而同時(shí)，這種跨兩個(gè)外顯子的讀段在分析轉(zhuǎn)錄本的剪接形式和研究選擇性剪接中有重要的作用為了解決這一問(wèn)題，人們采取兩種典型的策略來(lái)進(jìn)行接合區(qū)讀段的定位：一是根據(jù)已知的基因外

20、顯子注釋?zhuān)瑯?gòu)建所有可能的外顯子接合區(qū)序列，與基因組序列一并作為定位的參考基因組；二是不依賴(lài)基因注釋?zhuān)窍壤媚芡暾ㄎ坏交蚪M的讀段得到粗略的外顯子區(qū)域，并結(jié)合剪接位點(diǎn)序列構(gòu)建出可能的剪接位點(diǎn)，然后將不能完整定位的讀段分段定位到兩個(gè)外顯子可能的結(jié)合區(qū)域Illumina/Solexa平臺(tái)提供的RNA-seq軟件分析包GApipeline采用了第一種策略采用第二種策略的軟件有Tophat36和G-Mo.R-Se37等，最新的Tophat軟件增加了利用已知外顯子邊界注釋信息的選項(xiàng)838生物化學(xué)與生物物理進(jìn)展不論是哪種測(cè)序平臺(tái)，測(cè)序中都不可避免地存在一定的錯(cuò)誤，基因組中又存在單核苷酸多態(tài)性等引起的序

21、列變化，所以在讀段定位時(shí)通常允許一定數(shù)量的錯(cuò)配，可以根據(jù)不同應(yīng)用調(diào)節(jié)允許錯(cuò)配的程度另一方面，由于基因組中重復(fù)序列和高相似度序列的影響，某些讀段會(huì)出現(xiàn)定位到基因組多個(gè)位置的情況這些因素影響了各個(gè)讀段到基因組的定位質(zhì)量，在一些新的讀段定位算法中，同時(shí)給出每個(gè)讀段與基因組匹配質(zhì)量通常在后續(xù)處理前，人們將多定位的讀段都過(guò)濾掉，也有人嘗試用適當(dāng)?shù)牟呗园讯喽ㄎ蛔x段“分配”到其中某些位置上2,38.讀段定位到基因組后通常采用SAM(SequenceAlignment/Map)格式或其二進(jìn)制版本BAM格式39來(lái)存儲(chǔ)二進(jìn)制版本可大大節(jié)省存儲(chǔ)空間，但不能直接用普通文本編輯工具顯示關(guān)于SAM格式的詳細(xì)介紹，可查閱(

22、/SAM1.pdf)援2基因表達(dá)水平估計(jì)在深度測(cè)序技術(shù)出現(xiàn)之前，高通量測(cè)量不同基因表達(dá)水平的主要手段是基因芯片，在此基礎(chǔ)上可以對(duì)不同組織或者不同發(fā)育階段的基因表達(dá)差異和模式進(jìn)行分析mRN/doc/be6b0724c324a6024ee96d8bA-seq數(shù)據(jù)最基本的應(yīng)用也是檢測(cè)基因的表達(dá)水平，與基因芯片數(shù)據(jù)相比，RNA測(cè)序得到的是數(shù)字化的表達(dá)信號(hào)，具有靈敏度高、分辨率高、無(wú)飽和區(qū)等優(yōu)勢(shì)40-42測(cè)序數(shù)據(jù)是對(duì)提取出的RNA轉(zhuǎn)錄本中隨機(jī)進(jìn)行的短片段測(cè)序，如果一個(gè)轉(zhuǎn)錄本的豐度高，則測(cè)序后定位到其對(duì)應(yīng)的基因組區(qū)域的讀段也就多，可以通過(guò)對(duì)定位到基因外顯子區(qū)的讀段計(jì)數(shù)來(lái)估計(jì)基因表達(dá)水平很顯然，讀段計(jì)數(shù)除了

23、與基因真實(shí)表達(dá)水平成正比，還與基因長(zhǎng)度成正比，同時(shí)也與測(cè)序深度即測(cè)序?qū)嶒?yàn)中得到的總讀段數(shù)正相關(guān)為了保持對(duì)不同基因和不同實(shí)驗(yàn)間估計(jì)的基因表達(dá)值的可比性，人們提出了RPM和RPKM的概念2RPM(readspermillionreads)即每百萬(wàn)讀段中來(lái)自于某基因的讀段數(shù)，考慮了測(cè)序深度對(duì)讀段計(jì)數(shù)的影響RPKM(readsperkilobasespermillionreads)是每百萬(wàn)讀段中來(lái)自于某基因每千堿基長(zhǎng)度的讀段數(shù)，公式表示為：=基因區(qū)讀段計(jì)數(shù)伊伊109不僅對(duì)測(cè)序深度作了歸一化，而且對(duì)基.Biochem.Biophys.2010;37(8)因長(zhǎng)度也作了歸一化，使得不同長(zhǎng)度的基因在不同測(cè)序深

24、度下得到的基因表達(dá)水平估計(jì)值具有了可比性，是目前最常用的基因表達(dá)估計(jì)方法軟件rSeq43、DEGseq軟件包44和Cufflinks45等都提供了用上述方法進(jìn)行基因表達(dá)水平計(jì)算的功能根據(jù)RNA-seq文庫(kù)制備標(biāo)準(zhǔn)，在不考慮基因結(jié)構(gòu)的理想情況下，讀段會(huì)均勻地分布在基因上而實(shí)際上，通過(guò)對(duì)實(shí)際數(shù)據(jù)的可視化分析很容易發(fā)現(xiàn)，讀段在基因上的分布有著自身的一些模式，呈現(xiàn)出不均勻性(圖3)這一問(wèn)題已經(jīng)引起很多學(xué)者的關(guān)注46-48造成讀段分布出現(xiàn)偏好的原因可能有多個(gè)方面：在制備cDNA文庫(kù)時(shí)，反轉(zhuǎn)錄所采用的隨機(jī)引物對(duì)RNA序列具有一定的偏好性，使得cDNA片段不能夠完全均勻地取自各轉(zhuǎn)錄本；在PCR擴(kuò)增中，擴(kuò)增效

25、率與序列的GC含量等特征相關(guān)，可導(dǎo)致GC含量高的cDNA片段在文庫(kù)中拷貝數(shù)增加超過(guò)其他片段；舍棄多定位的讀段也可能導(dǎo)致讀段的非均勻分布；等等如果能對(duì)讀段分/doc/be6b0724c324a6024ee96d8b布的不均勻性進(jìn)行建模并加以校正，可以提高RNA-seq推斷基因表達(dá)量的準(zhǔn)確度但根據(jù)對(duì)實(shí)際數(shù)據(jù)的觀察，對(duì)于較長(zhǎng)轉(zhuǎn)錄本，讀段非均勻分布帶來(lái)的誤差很大程度上可相互抵消，用RPKM來(lái)估計(jì)基因的表達(dá)水平可以得到比較滿(mǎn)意的結(jié)果援3選擇性剪接事件識(shí)別和剪接異構(gòu)體表達(dá)水平推斷在真核生物中，選擇性剪接現(xiàn)象普遍存在基因轉(zhuǎn)錄形成的mRNA前體(pre-mRNA)在剪接過(guò)程中因去掉不同的內(nèi)含子區(qū)域或保留不同的

26、外顯子區(qū)域，可形成不同的剪接異構(gòu)體根據(jù)RNA-seq原理，只要測(cè)序深度足夠深，就能檢測(cè)到所有轉(zhuǎn)錄本的全部序列，包括來(lái)自剪接接合區(qū)的序列利用考慮到接合區(qū)的讀段定位方法，就有可能系統(tǒng)地研究某一組織或某一條件下的基因選擇性剪接事件前面已經(jīng)提到，Tophat等軟件定位剪接接合區(qū)讀段的策略能標(biāo)定出剪接事件中的兩個(gè)剪接位點(diǎn)：供體位點(diǎn)和受體位點(diǎn)通過(guò)比較供體位點(diǎn)和受體位點(diǎn)的組合，就能識(shí)別選擇性剪接事件4,49圖3中包含了選擇性剪接識(shí)別的一個(gè)例子進(jìn)一步，通過(guò)對(duì)供體和受體位點(diǎn)的讀段計(jì)數(shù)，結(jié)合外顯子其他區(qū)域的讀段數(shù)據(jù)，還能定量地計(jì)算選擇性剪接事件之間的比例50-51對(duì)于每一個(gè)剪接異構(gòu)體，RNA-seq數(shù)據(jù)能在一;3

27、7(8)王曦等：新一代高通量RNA測(cè)序數(shù)據(jù)的處理與分析839定程度上推斷其表達(dá)水平比如，可以根據(jù)已知外顯子組成和各外顯子長(zhǎng)度對(duì)剪接異構(gòu)體建立數(shù)學(xué)模型，在測(cè)序讀段轉(zhuǎn)錄本上均勻分布的假設(shè)下，利用各外顯子上的讀段數(shù)和接合區(qū)讀段數(shù)求解異構(gòu)體的表達(dá)值Jiang等43的方法及軟件IsoInfer52和cufflinks45都采用了這種思路來(lái)實(shí)現(xiàn)剪接異構(gòu)體的表達(dá)推斷需要指出的是，某些形式的剪接異構(gòu)體表達(dá)水平在這種方法框架下不可辨識(shí)533援4新基因的檢測(cè)在對(duì)RNA-seq數(shù)據(jù)的分析中，人們發(fā)現(xiàn)，往往不是所有讀段都能定位到已有注釋的基因區(qū)，說(shuō)明除了轉(zhuǎn)錄噪聲或測(cè)序錯(cuò)誤等的影響外，可能還存在尚未被注釋的基因這里，我

28、們把這種尚未注釋的基因稱(chēng)為新基因，包括新的蛋白質(zhì)編碼基因和非編碼/doc/be6b0724c324a6024ee96d8bRNA基因能檢測(cè)新基因，尤其是低表達(dá)基因是RNA-seq技術(shù)優(yōu)于基因芯片的特點(diǎn)之一，因?yàn)樗恍枰靡阎蜃⑨寔?lái)設(shè)計(jì)檢測(cè)探針seq技術(shù)靈敏度高，但樣品污染、測(cè)序錯(cuò)誤等仍可能帶來(lái)背景噪聲從基因組未注釋區(qū)域的RNA測(cè)序讀段信號(hào)中檢測(cè)新基因是典型的信號(hào)檢測(cè)問(wèn)題如何控制新基因識(shí)別的誤發(fā)現(xiàn)率(FDR)是檢測(cè)方法的關(guān)鍵Useq軟件包54將ChIP-seq數(shù)據(jù)分析的方法移植到RNA-seq數(shù)據(jù)上，用滑窗的方法來(lái)識(shí)別測(cè)序讀段定位富集的區(qū)域，給出反映滑窗所在區(qū)域讀段富集顯著程度的P值(P-

29、value)及新基因誤發(fā)現(xiàn)率，通過(guò)設(shè)定P值或誤發(fā)現(xiàn)率的閾值，可篩選出讀段富集的區(qū)域，再將相鄰區(qū)域合并或根據(jù)剪接接合區(qū)讀段將相應(yīng)區(qū)域連接，完成新基因的檢測(cè)援5讀段的可視化及注釋對(duì)于復(fù)雜的組學(xué)數(shù)據(jù)，能盡可能方便地直接觀察數(shù)據(jù)對(duì)于數(shù)據(jù)的分析和解釋都非常重要，對(duì)新一代測(cè)序數(shù)據(jù)的可視化和交互展示是一個(gè)非常重要但容易被人忽視的問(wèn)題不深入考查數(shù)據(jù)的細(xì)節(jié)，而是滿(mǎn)足于對(duì)數(shù)據(jù)的統(tǒng)計(jì)分析，是高通量數(shù)據(jù)應(yīng)用中經(jīng)常容易陷入的誤區(qū)，方便有效的可視化工具能夠幫助避免這樣的誤區(qū)表2列出了部分適用于RNA-seq數(shù)據(jù)的全基因組瀏覽器，其中比較具有代表性的有UCSCGenomeBrowser、CisGenomeBrowser和I

30、GV(IntegrativeGenomicsViewer)等這些瀏覽器具有如下特點(diǎn)：a能在不同尺度下顯現(xiàn)單個(gè)或多個(gè)讀段在基因組上的位置，包括來(lái)源于剪接接合區(qū)的讀段；b能在不同尺度下顯示不同區(qū)域的讀段豐度，以反映不同區(qū)域的轉(zhuǎn)錄水平或測(cè)序效率；c能顯示基因及其剪接異構(gòu)體的注釋信息；d能顯示其他注釋信息，例如物種間基因組序列保守性、序列GC含量等；e能直接或間接支持SAM/BAM讀段定位數(shù)據(jù)存儲(chǔ)格式.UCSCGenomeBrowser55屬于基于網(wǎng)絡(luò)模式的全基因組瀏覽器，所有數(shù)據(jù)都需要上載到遠(yuǎn)程服務(wù)器，經(jīng)過(guò)處理后將圖形返回客戶(hù)端顯示圖3中的例子就是從UCSCGenomeBrowser的顯示截取的.C

31、isGenomeBrowser56是典型的本地版基因組瀏覽器，所有讀段數(shù)據(jù)、注釋信息都存于本地文件，因此不/doc/be6b0724c324a6024ee96d8b需要網(wǎng)絡(luò)連接，方便內(nèi)部考查數(shù)據(jù)用IGV(/igv)可以說(shuō)是以上兩種模式的融合，既可以從遠(yuǎn)程服務(wù)器端下載各種注釋信息，又可以從本地加載注釋信息表2名稱(chēng)IGV適用于mRNA鄄seq數(shù)據(jù)的全基因組瀏覽器/viewersapplicabletomRNA鄄seqdataviewing支持的數(shù)據(jù)格式網(wǎng)址:/igv/GFF3,BED,SAM/BAM,WIG,55BED,bigBed,BEDGRAPH,GFF,GTF,bigWig,BAM,5657

32、,BED,refFlat,FA,Wig,BED,GFF,FASTA,ELAND,GFF,BED,MAQ,SAMBAM,BED,GFF2,GFF3,FASTA,VCFExpressiondata,Annotationtracks:/jiangh/browser/sj/mochiview-start:/www.bioinformatics.bbsrc.ac.uk/projects/seqmonk/p/gambit-viewer:/packages/release/bioc/html/GenomeGraphs.html以上列出的全基因組瀏覽器均可/doc/be6b0724c324a6024ee96d

33、8b在Windows、Linux和蘋(píng)果公司的MacOS等計(jì)算機(jī)平臺(tái)下運(yùn)行.840生物化學(xué)與生物物理進(jìn)展Prog.Biochem.Biophys.2010;37(8)剪接接合區(qū)樣本A測(cè)序標(biāo)簽分布該基因在A、B樣本中差異表達(dá)內(nèi)含子區(qū)的測(cè)序標(biāo)簽接測(cè)序標(biāo)布剪接接合區(qū)樣本B測(cè)序標(biāo)簽分布基因注釋圖3鄄seq數(shù)據(jù)可視化示例鄄seqdatavisualization圖示區(qū)域?yàn)槿祟?lèi)基因CBX7.圖中紅色表示樣本A的數(shù)據(jù)，藍(lán)色表示樣本B.各軌道(track)依次為:基因組坐標(biāo)、樣本A的剪接接合區(qū)、樣本A的讀段分布、樣本B的剪接接合區(qū)、樣本B的讀段分布、UCSC基因注釋.圖中還標(biāo)識(shí)了：因受體位點(diǎn)不同而形成的選擇性剪

34、接;基因的5憶端出現(xiàn)讀段的非均勻分布;在兩個(gè)樣本中,差異表達(dá)基因的讀段信號(hào)強(qiáng)度不同;在基因標(biāo)注的內(nèi)含子(intron)區(qū)域存在少量不連續(xù)的讀段.除對(duì)讀段的可視化外，用描述統(tǒng)計(jì)學(xué)方法對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分類(lèi)別統(tǒng)計(jì)也十分重要例如，統(tǒng)計(jì)讀段在各個(gè)染色體上的分布情況和在注釋的外顯子、內(nèi)含子、剪接接合區(qū)、基因間區(qū)的分布情況等目前，已經(jīng)有一些用于測(cè)序數(shù)據(jù)注釋的生物信息學(xué)軟件，比如SAMtools39、BEDtools58等，但由于測(cè)序技術(shù)發(fā)展迅速，用戶(hù)需求因人而異，用戶(hù)經(jīng)常還需要根據(jù)需求編寫(xiě)一定的程序或腳本完成或完善注釋分析的任務(wù)對(duì)于熟悉圖形用戶(hù)界面的研究人員，還可以利用UCSCTableBrowse56和Ga

35、laxy59-60來(lái)配合完成注釋分析由于UCSCTableBrowser集成了大量基因組尺度上的注釋信息，而Galaxy又為用戶(hù)提供了書(shū)寫(xiě)簡(jiǎn)單、接口明晰和直觀的數(shù)據(jù)處理流程，這種方法十分方便有效，也為很多學(xué)者在展示研究成/doc/be6b0724c324a6024ee96d8b果時(shí)所采用以上描述了對(duì)于基因組已知的物種進(jìn)行RNA-seq數(shù)據(jù)處理基本流程，圖2a給出了其主要步驟若研究對(duì)象尚未完成基因組測(cè)序，則需要采用讀段的從頭拼裝(denovoassembly)6,61-62來(lái)代替讀段定位，后續(xù)流程也須做相應(yīng)的調(diào)整若RNA-seq實(shí)驗(yàn)在文庫(kù)制備時(shí)保留了RNA的方向信息，則應(yīng)分別研究來(lái)自正鏈和反鏈的

36、轉(zhuǎn)錄產(chǎn)物，并通過(guò)與基因注釋比較來(lái)檢測(cè)反義轉(zhuǎn)錄本63最近，RNA-MATE64軟件在其分析流程中加入如此的處理策略此外，通過(guò)分析定位到外顯子接合區(qū)的讀段，還可以獲取轉(zhuǎn)錄本結(jié)構(gòu)，這為研究基因的剪接調(diào)控機(jī)理提供了重要信息5而利用RNA-seq數(shù)據(jù)提供的序列信息，通過(guò)與DNA序列的細(xì)致比較可分析轉(zhuǎn)錄組的序列差異(如SNP等)65，從而研究等位基因的表達(dá)模式66-67及RNA編輯68等最后需要指出，由于miRNA在序列和結(jié)構(gòu)上具有一定的特點(diǎn)，miRNA-seq數(shù)據(jù)的基本處理流程也與本節(jié)所述有所不同，感興趣的讀者可參考軟件工具miRDeep69提供的處理策略多類(lèi)樣本mRNA鄄seq數(shù)據(jù)間的比較分析很多RN

37、A-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)度等因素的影響70-71，需要對(duì)結(jié)果進(jìn)行仔細(xì)分析，消除可能的混雜因素，必要時(shí)可以用讀段的絕對(duì)表達(dá)值倍數(shù)變化(fold-change)來(lái)作為補(bǔ)充圖2b給出了兩類(lèi)樣本數(shù)據(jù)分析的框架;37(8)王曦等：新一代高通量RNA測(cè)序數(shù)據(jù)的處理與分析841援1差異表達(dá)基因的識(shí)別雖然新一代測(cè)序相對(duì)第一代測(cè)序的單

38、位成本大大降低，但是，利用RNA測(cè)序進(jìn)行基因表達(dá)研究的成本仍很高，因此，很多實(shí)驗(yàn)室沒(méi)有條件進(jìn)行樣本重復(fù)如果兩類(lèi)樣本均沒(méi)有生物重復(fù)，例如只/doc/be6b0724c324a6024ee96d8b對(duì)兩個(gè)細(xì)胞系各進(jìn)行一次mRNA樣本測(cè)序，則可以用隨機(jī)采樣模型通過(guò)假設(shè)檢驗(yàn)來(lái)分析差異表達(dá)對(duì)于某個(gè)基因，如果一個(gè)讀段來(lái)自于這個(gè)基因，我們稱(chēng)事件A發(fā)生對(duì)于一次RNA-seq實(shí)驗(yàn)，事件A發(fā)生的概率可以用這個(gè)基因上的讀段數(shù)n除以所有基因上的讀段總數(shù)N來(lái)估計(jì)，即RPM事件A發(fā)生的概率反應(yīng)了這個(gè)基因的表達(dá)水平如果要判斷(a)某個(gè)基因在兩個(gè)樣本中的表達(dá)水平是否一致，就可以通過(guò)檢驗(yàn)事件A在兩種條件下發(fā)生的概率是否一致來(lái)實(shí)

39、現(xiàn)，采用似然比檢驗(yàn)1、Fisher精確檢驗(yàn)72以及基于MA圖的統(tǒng)計(jì)檢驗(yàn)方法44等同樣，也可用RPKM作為統(tǒng)計(jì)量來(lái)進(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è)公開(kāi)發(fā)表的軟件，包括DEGseq44、Useq54、Cufflinks45中的Cuffdiff模塊等圖4展示了我們開(kāi)發(fā)的DEGseq軟件提供的多種差異表達(dá)基因識(shí)別方法的應(yīng)用例子(b)-A(d)A-0A8log2(readcountsforeachgene)inB(c)-AvsBAvsB圖4用DEGseq軟件包

40、識(shí)別差異表達(dá)基因的結(jié)果(a)各基因在樣本A和樣本B中表達(dá)水平的散點(diǎn)圖.(b),(c),(d)圖中紅點(diǎn)表示分別用FET、LRT和MARS方法得到的差異表達(dá)基因.FET:FishersExactTest,Fisher精確檢驗(yàn).LRT:LikelihoodRatioTest,似然比檢驗(yàn).MARS:MA-plot-basedmethodwit/doc/be6b0724c324a6024ee96d8bhRandomSamplingmodel,基于MA圖的隨機(jī)采樣模型.如果每一類(lèi)樣本都包含了若干生物重復(fù)，如病人和正常人對(duì)照研究，則可以沿用基因芯片數(shù)據(jù)分析中的很多方法比如，可以用t檢驗(yàn)結(jié)合倍數(shù)變化的方法來(lái)分

41、析差異表達(dá)如果兩類(lèi)樣本具有配對(duì)的信息，也可以通過(guò)整合每對(duì)樣本分析結(jié)果來(lái)實(shí)現(xiàn)其步驟為，先在每對(duì)樣本中識(shí)別出差異表達(dá)的基因，再尋找這若干組差異表達(dá)基因之間的相同者，或用投票的方法來(lái)為基因的差異程度打分針對(duì)某些RNA-seq數(shù)據(jù)生物樣本量小，R軟件包DEGseq44和edgeR73等還專(zhuān)門(mén)提供了基于改進(jìn)模型的統(tǒng)計(jì)方法此外，一類(lèi)將分類(lèi)器與特征選擇包裹在一起的方法也同樣適用于此類(lèi)問(wèn)題(見(jiàn)4.3)842生物化學(xué)與生物物理進(jìn)展4援2差異表達(dá)剪接異構(gòu)體的識(shí)別差異表達(dá)剪接異構(gòu)體的識(shí)別方法與差異表達(dá)基因的識(shí)別相似如果把剪接異構(gòu)體看成是獨(dú)立的基因，那么前面討論的用于識(shí)別差異表達(dá)基因的方法對(duì)剪接異構(gòu)體完全適用但是，注

42、意到來(lái)自于同一個(gè)基因的剪接異構(gòu)體并不獨(dú)立，某些假設(shè)檢驗(yàn)的基本條件并不滿(mǎn)足，得到的結(jié)果就不一定正確此外，由于現(xiàn)在剪接異構(gòu)體表達(dá)推斷的方法還不夠成熟，加之在基因結(jié)構(gòu)不可辨識(shí)的剪接異構(gòu)體上作表達(dá)推斷會(huì)出現(xiàn)病態(tài)結(jié)果，差異表達(dá)剪接異構(gòu)體的識(shí)別問(wèn)題還處于探索的階段目前，在剪接異構(gòu)體表達(dá)水平可辨識(shí)且讀段覆蓋度較高的基因上，BASIS74方法通過(guò)貝葉斯模型來(lái)推斷差異表達(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á)的差異75外顯子表達(dá)差

43、異的分析可以完全利用基因表達(dá)差異的分析方法而剪接接合區(qū)也可以看成是一個(gè)較短的“外顯子”(長(zhǎng)度一般與測(cè)序長(zhǎng)度相當(dāng))不過(guò)，由于外顯子長(zhǎng)度較基因的長(zhǎng)度短，對(duì)應(yīng)的讀段數(shù)量較少，差異識(shí)別的敏感度會(huì)有所下降Solas方法75就是根據(jù)類(lèi)似的原理，采用統(tǒng)計(jì)學(xué)假設(shè)檢驗(yàn)的方法來(lái)識(shí)別差異表達(dá)的剪接異構(gòu)體4援3對(duì)樣本的分類(lèi)分析:/doc/be6b0724c324a6024ee96d8b通過(guò)統(tǒng)計(jì)方法識(shí)別出來(lái)的差異表達(dá)基因及剪接異構(gòu)體能否有效地區(qū)別兩類(lèi)樣本，可以通過(guò)分類(lèi)分析進(jìn)一步證實(shí)如果把每個(gè)基因(或剪接異構(gòu)體)的表達(dá)值作為特征，則差異表達(dá)基因(或剪接異構(gòu)體)的選取也就是特征篩選的過(guò)程把前面用統(tǒng)計(jì)方法等檢測(cè)出來(lái)的差異表達(dá)

44、基因(或剪接異構(gòu)體)用于分類(lèi)分析，常被稱(chēng)為過(guò)濾法另一類(lèi)基于分類(lèi)器的包裹法，例如R-SVM76、SVM-RFE77等，可以根據(jù)每個(gè)特征在分類(lèi)器中所占的權(quán)重來(lái)篩選特征，因此也可以用于差異表達(dá)基因(或剪接異構(gòu)體)的識(shí)別分類(lèi)的性能可以用交叉驗(yàn)證(cross-validation，CV)方法來(lái)評(píng)估需要特別注意的是，交叉驗(yàn)證應(yīng)該包括對(duì)特征選擇步驟的交叉驗(yàn)證，防止發(fā)生信息泄露而導(dǎo)致評(píng)估結(jié)果過(guò)于樂(lè)觀具體做法是：將樣本按一定的策略分成兩份，一份(通常是樣本數(shù)多的一份)用于特征選取和分類(lèi)器訓(xùn)練，而用余下的樣本進(jìn)行分類(lèi)器性能的估計(jì)；重復(fù)以上步驟多次，就得.Biochem.Biophys.2010;37(8)到交叉驗(yàn)

45、證錯(cuò)誤率必要時(shí)還可以用隨機(jī)置換檢驗(yàn)(permutationtest)來(lái)推斷所得錯(cuò)誤率的統(tǒng)計(jì)顯著性76當(dāng)樣本數(shù)較小時(shí)，可以采用留一法交叉驗(yàn)證(leave-one-outcross-validation，LOOCV)4援4其他高層分析方法檢測(cè)差異表達(dá)的基因或差異表達(dá)異構(gòu)體是人們認(rèn)識(shí)所研究的生物問(wèn)題機(jī)理的第一步，接下來(lái)需要從功能上研究這些差異轉(zhuǎn)錄現(xiàn)象的分子機(jī)理這與在基因芯片應(yīng)用中所面臨的是同樣的生物學(xué)問(wèn)題，對(duì)芯片數(shù)據(jù)分析結(jié)果的后續(xù)處理方法，都可以借鑒到測(cè)序數(shù)據(jù)上來(lái)如何進(jìn)一步地從機(jī)理來(lái)解釋結(jié)果，還需結(jié)合已知生物學(xué)知識(shí)進(jìn)行后續(xù)分析人們對(duì)基因芯片得到的基因表達(dá)數(shù)據(jù)進(jìn)行分析的很多方法都可以用到RNA-seq

46、數(shù)據(jù)上來(lái)，比如利用機(jī)器學(xué)習(xí)方法進(jìn)行分類(lèi)和特征選擇，對(duì)差異表達(dá)的基因進(jìn)行GO(geneontology)78類(lèi)別富集分析、信號(hào)通路富集分析等，一些常用的分析工具包括GoMiner79、DAVID80和VisANT81等需要說(shuō)明的是，在各種以差異表達(dá)基因?yàn)榛A(chǔ)的分析中，由于基因表達(dá)水平都是通過(guò)讀段計(jì)數(shù)來(lái)估計(jì)的，表達(dá)水平較高或/doc/be6b0724c324a6024ee96d8b轉(zhuǎn)錄本較長(zhǎng)的基因擁有更多的讀段，更容易被多數(shù)統(tǒng)計(jì)方法識(shí)別為差異表達(dá)基因70這種偏好可能對(duì)后續(xù)分析帶來(lái)影響以GO類(lèi)別富集分析為例，這種偏好將導(dǎo)致長(zhǎng)基因占主導(dǎo)的功能類(lèi)別更有可能被識(shí)別為富集的功能這將對(duì)生物機(jī)理的研究帶來(lái)誤導(dǎo)最

47、近，Young等71發(fā)展了一種GOseq方法，針對(duì)這一偏好對(duì)GO類(lèi)別富集分析做了改進(jìn)RNA鄄seq數(shù)據(jù)處理中的生物信息學(xué)挑戰(zhàn)高通量測(cè)序技術(shù)的發(fā)展十分迅速，這要求相應(yīng)的數(shù)據(jù)處理與分析方法快速跟進(jìn)正是這些方法，架起了高通量實(shí)驗(yàn)數(shù)據(jù)與科學(xué)問(wèn)題之間的橋梁這種橋梁作用正日趨重要，也為生物信息學(xué)帶來(lái)了挑戰(zhàn)7,71這里，我們重點(diǎn)討論兩方面的挑戰(zhàn)：a如何實(shí)現(xiàn)剪接接合區(qū)讀段的準(zhǔn)確定位？b在數(shù)據(jù)處理各階段中，如何對(duì)RNA-seq數(shù)據(jù)的系統(tǒng)誤差和固有偏好建?；蜓a(bǔ)償，以消除它們可能帶來(lái)的錯(cuò)誤推斷及結(jié)論？援1剪接接合區(qū)讀段的定位測(cè)序技術(shù)的一個(gè)發(fā)展趨勢(shì)是測(cè)序長(zhǎng)度不斷增加隨著讀長(zhǎng)的增加，RNA-seq中來(lái)自剪接接合區(qū)的讀段

48、會(huì)越來(lái)越多我們粗略估算，按照人類(lèi)基因組refSeq基因注釋?zhuān)话闱闆r下，如果測(cè)序讀長(zhǎng);37(8)王曦等：新一代高通量RNA測(cè)序數(shù)據(jù)的處理與分析843為50個(gè)堿基，則約有10%的讀段來(lái)自剪接接合區(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)的早期，人們沒(méi)有意識(shí)到剪接接合區(qū)讀段的重要性因?yàn)楫?dāng)時(shí)的讀長(zhǎng)只有2030個(gè)堿基，來(lái)自剪接接合區(qū)的讀段所占比例甚小當(dāng)時(shí)讀段定位的通常做法是，先將讀段與全基因組序列做映射定位，再考慮不能定位的讀段是否來(lái)自于剪接接合區(qū)2

49、這種做法雖然在一定程度上保證了讀段定位的比率，但由于基因組中重復(fù)序列和相似序列的存在，部分接合區(qū)讀段有可能在容許錯(cuò)配的情況下被定位到基因組上其他位置，從而失去了定位到正確的剪接接合區(qū)的機(jī)會(huì)在讀段定位時(shí)，如果要同時(shí)/doc/be6b0724c324a6024ee96d8b考慮基因組序列和剪接接合區(qū)序列，就要利用已知的剪接事件注釋?zhuān)@是目前軟件通用的方法然而，包括人類(lèi)在內(nèi)的各物種的基因注釋信息都還有待完善，也沒(méi)有較完整的剪接組(splicome)數(shù)據(jù)庫(kù)，能夠不依賴(lài)注釋信息和對(duì)剪接機(jī)理的現(xiàn)有認(rèn)識(shí)，高效、準(zhǔn)確地定位所有已知和未知的接合區(qū)讀段，仍然是對(duì)讀段映射定位算法的一個(gè)挑戰(zhàn)援2系統(tǒng)噪聲和偏好的分析雖

50、然深度測(cè)序技術(shù)的準(zhǔn)確性較以前的技術(shù)有了很大提高，但仍然存在錯(cuò)誤和噪聲比如從圖3中可以看到，內(nèi)含子區(qū)內(nèi)有一些不連續(xù)的讀段，很可能由系統(tǒng)噪聲造成，如樣品污染、測(cè)序錯(cuò)誤和不恰當(dāng)?shù)淖x段定位策略等從圖3還能看出，外顯子區(qū)域內(nèi)的讀段信號(hào)分布也很不均勻有文獻(xiàn)報(bào)道，序列組成尤其是GC含量46、RNA二級(jí)結(jié)構(gòu)2等也有可能是導(dǎo)致讀段不均勻分布的原因這些噪聲和分布偏好將影響新基因的識(shí)別和對(duì)剪接異構(gòu)體形式和表達(dá)水平的推斷合理地建模RNA-seq數(shù)據(jù)中的系統(tǒng)噪聲和偏好是解決上述問(wèn)題最有效的辦法基本的思路可以是：首先根據(jù)實(shí)驗(yàn)原理尋找可能產(chǎn)生系統(tǒng)噪聲或偏差的因素，并盡可能將這些因素轉(zhuǎn)化成可量化的特征，如序列特征、二級(jí)結(jié)構(gòu)等

51、；然后，將用實(shí)驗(yàn)數(shù)據(jù)對(duì)這些特征做統(tǒng)計(jì)分析，構(gòu)造和訓(xùn)練模型，用模型來(lái)對(duì)數(shù)據(jù)進(jìn)行校正需要注意的是，某些偏好是由當(dāng)前的測(cè)序技術(shù)和分析方法共同造成的，難以完全消除71在這種情況下，后續(xù)處理和解釋時(shí)需要充分意識(shí)到這種偏好可能對(duì)生物學(xué)結(jié)論帶來(lái)的影響，必要時(shí)通過(guò)補(bǔ)充其他實(shí)驗(yàn)來(lái)驗(yàn)證和修正通過(guò)高通量測(cè)序得到的生物學(xué)結(jié)論總結(jié)與展望本文以Illumina/Solexa測(cè)序平臺(tái)為例，嘗試對(duì)新一代測(cè)序技術(shù)的RNA-seq數(shù)據(jù)處理和分析方法做了較為全面的梳理，并對(duì)各個(gè)環(huán)節(jié)上可用的軟件進(jìn)行了匯總高通量測(cè)序是正在飛速發(fā)展的技術(shù)，相應(yīng)的生物信息學(xué)方法也在快速發(fā)展，這里討論的是RNA-seq中一些代表性的方法和問(wèn)題，希望能對(duì)正在

52、或即將采用RNA-seq實(shí)驗(yàn)進(jìn)行科學(xué)研究的學(xué)者和進(jìn)行RNA測(cè)序數(shù)據(jù)處理的同行提供參考.測(cè)序和基因芯片有很多共同的應(yīng)用領(lǐng)域，盡管相對(duì)還不是很成熟，RNA-seq技術(shù)在很多方面已經(jīng)表現(xiàn)/doc/be6b0724c324a6024ee96d8b出了優(yōu)勢(shì)，有人甚至預(yù)言基因芯片時(shí)代即將結(jié)束36但也有報(bào)道認(rèn)為，RNA-seq數(shù)據(jù)在基因表達(dá)水平的估計(jì)上和基因芯片相比沒(méi)有明顯的優(yōu)勢(shì)60，加上測(cè)序的成本目前還遠(yuǎn)高于芯片實(shí)驗(yàn)的成本，所以更多人認(rèn)為測(cè)序和基因芯片將長(zhǎng)期共存，以各自不同的特點(diǎn)在現(xiàn)代組學(xué)研究中發(fā)揮作用新一代高通量測(cè)序技術(shù)的應(yīng)用面非常廣82，RNA-seq只是其中一個(gè)方面，除此之外，基因組的從頭測(cè)序和重測(cè)

53、序83-84、染色質(zhì)免疫沉淀測(cè)序(ChIP-seq)85-86、甲基化測(cè)序(Methyl-seq)87-88等技術(shù)都同樣有著廣泛的應(yīng)用尤其是，用ChIP-seq研究蛋白質(zhì)與DNA的相互作用，能夠得到高分辨率的轉(zhuǎn)錄因子結(jié)合數(shù)據(jù)和組蛋白修飾等表觀遺傳學(xué)數(shù)據(jù)發(fā)展有效的生物信息學(xué)方法，將ChIP-seq數(shù)據(jù)與RNA-seq得到的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行綜合分析，將大大推進(jìn)人們對(duì)復(fù)雜的基因轉(zhuǎn)錄調(diào)控系統(tǒng)的認(rèn)識(shí)致謝感謝本實(shí)驗(yàn)室劉霖曦、謝芃、孟璐等同學(xué)對(duì)本工作有意義的討論，感謝斯坦福大學(xué)WingHWong教授、HuiJiang博士和JunLi同學(xué)等的討論和幫助參考文獻(xiàn)1MarioniJC,MasonCE,ManeSM,

54、etal.RNA-seq:anassessment.GenomeRes,2008,18(9):1509-15172MortazaviA,WilliamsBA,McCueK,etal.MappingandSeq.NatMethods,2008,5(7):621-628844生物化學(xué)與生物物理進(jìn)展3NagalakshmiU,WangZ,WaernK,etal.Thetranscriptional:/doc/be6b0724c324a6024ee96d8bing.Science,2008,320(5881):1344-13494SultanM,SchulzMH,RichardH,etal.Aglob

55、alviewofgene.Science,2008,321(5891):956-9605WangET,SandbergR,LuoS,etal.Alternativeisoformregulation.Nature,2008,456(7221):470-4766BirzeleF,SchaubJ,RustW,etal.Intotheunknown:expression.NucleicAcidsRes,2010,doi:10.1093/nar/gkq1167SangerF,NicklenS,CoulsonAR.DNAsequencingwithchain-terminatinginhibitors.

56、ProcNatlAcadSciUSA,1977,74(12):5463-54678MarguliesM,EgholmM,AltmanWE,etal.Genomesequencingdensitypicolitrereactors.Nature,2005,437(7057):376-3809ShendureJ,JiH.Next-generationDNAsequencing.NatBiotechnol,26(10):1135-114510RuparelH,BiL,LiZ,etal.Designandsynthesisofa3憶-O-ally/doc/be6b0724c324a6024ee96d8

57、bl.ProcNatlAcadSciUSA,2005,102(17):5932-593711SeoTS,BaiX,KimDH,etal.Four-colorDNAsequencingby.ProcNatlAcadSciUSA,2005,102(17):5926-593112JuJ,KimDH,BiL,etal.Four-colorDNAsequencingbysynthesis.ProcNatlAcadSciUSA,2006,103(52):19635-1964013FedurcoM,RomieuA,WilliamsS,etal.BTA,anovelreagentforphaseamplifi

58、edDNAcolonies.NucleicAcidsRes,2006,34(3):e2214ShendureJA,PorrecaGJ,ChurchGM.OverviewofDNA/AusubelFM,BrentR,KingstonRE,A:JohnWileyandSons,Inc.,2008:Unit7.115MardisER.Next-generationDNAsequencingmethods.AnnuRev,2008,9:387-40216Full/doc/be6b0724c324a6024ee96d8berCW,MiddendorfLR,BennerSA,etal.Thechallen

59、gesof.NatBiotechnol,2009,27(11):1013-102317CroucherNJ,FookesMC,PerkinsTT,etal.Asimplemethodfor.NucleicAcidsRes,2009,37(22):e14818ParkhomchukD,BorodinaT,AmstislavskiyV,etal.TranscriptomespecificsequencingofcomplementaryDNA.NucleicAcidsRes,2009,37(18):e12319PerkinsTT,KingsleyRA,FookesMC,etal.Astrand-s

60、pecificSeqanalysisofthetranscriptomeofthetyphoidbacillusSalmonellatyphi.PLoSGenet,2009,5(7):e100056920MamanovaL,AndrewsRM,JamesKD,etal.FRT-seq:free,strand-specifictranscriptomesequencing.Nat.Biochem.Biophys.;37(8),2010,7(2):130-13221CockPJ,FieldsCJ,GotoN,etal.TheSangerFASTQfileformat,andtheSolexa/Il