石新新 劉學(xué)軍 張 禮

(南京航空航天大學(xué)計算機科學(xué)與技術(shù)學(xué)院，南京，210016)

?

改進的RNA-Seq數(shù)據(jù)轉(zhuǎn)錄組表達分析研究

石新新 劉學(xué)軍 張 禮

(南京航空航天大學(xué)計算機科學(xué)與技術(shù)學(xué)院，南京，210016)

基于高通量測序的RNA-Seq(RNA-sequencing)是用于轉(zhuǎn)錄組研究的一種新技術(shù)，針對該技術(shù)在轉(zhuǎn)錄組表達分析研究中存在的讀段多源映射和讀段非均勻分布等難點，提出一個改進的轉(zhuǎn)錄組表達研究方法LDASeqII(Improvement of latent Dirichlet allocation for sequencing data)。模型利用剪接異構(gòu)體結(jié)構(gòu)信息對參數(shù)進行約束并進行外顯子讀段數(shù)目歸一化處理，解決了讀段非均勻分布下的多源映射問題。通過引入“偽外顯子”和“偽轉(zhuǎn)錄本”分別處理接合區(qū)讀段和噪聲讀段。將模型應(yīng)用到真實數(shù)據(jù)集上，并與原LDASeq(Latent Dirichlet allocation for sequencing data)模型和目前流行的 Cufflinks與RSEM(RNA-Seq by expectation maximization)方法進行對比。結(jié)果顯示,改進方法獲得了更為準確的轉(zhuǎn)錄本及基因表達水平計算結(jié)果。

基因表達；RNA-Seq；轉(zhuǎn)錄組表達；多源映射；非均勻性

引 言

RNA-Seq(RNA sequencing)是基于高通量測序技術(shù)對轉(zhuǎn)錄組進行研究的一種新方法，具有信噪比高、分辨率高和所需樣本少等優(yōu)勢，可以用來檢測和量化任何物種的轉(zhuǎn)錄片段[1-2]，近年來在轉(zhuǎn)錄組研究中得到大量應(yīng)用。RNA-Seq數(shù)據(jù)處理過程主要分為3個方面：(1) 將讀段映射到參考基因組或轉(zhuǎn)錄組上；(2) 利用匹配到的讀段對基因或轉(zhuǎn)錄組進行構(gòu)建；(3) 通過計數(shù)匹配到基因或轉(zhuǎn)錄本上的讀段數(shù)目來計算基因或轉(zhuǎn)錄本表達值[3]。本文主要針對最后一個步驟進行研究，即給定轉(zhuǎn)錄組構(gòu)建、計算基因以及異構(gòu)體的表達水平。

RNA-Seq通過計數(shù)匹配到基因元件上讀段數(shù)目來表示相關(guān)基因的轉(zhuǎn)錄水平[4],為了保持對不同基因和不同實驗間估計的基因表達值的可比性，必須對基因長度和測序深度進行歸一化。最常用的歸一化方法是Mortazavi等人提出RPKM(Reads per kilobases per million reads)方法，即每百萬讀段中來自該基因的讀段數(shù)[5]?；蛟谵D(zhuǎn)錄過程中會形成一個或多個轉(zhuǎn)錄本，這些轉(zhuǎn)錄本共享基因的部分外顯子[6]，造成讀段到異構(gòu)體的多源映射，故只有一少部分比例的讀段能夠唯一確定其來源轉(zhuǎn)錄本[7]。因此，RPKM方法不能直接應(yīng)用于轉(zhuǎn)錄本表達值計算。一種解決思路為可以根據(jù)已知外顯子組成和各外顯子長度對轉(zhuǎn)錄本建立數(shù)學(xué)模型，利用基因外顯子上讀段數(shù)求解轉(zhuǎn)錄本表達值。文獻[8]提出的轉(zhuǎn)錄本表達值計算方法和文獻[9]設(shè)計的Cufflinks模型都采用了這種思路來實現(xiàn)剪接異構(gòu)體的表達推斷。此外，還有文獻[10-11]提出的RSEM(RNA-Seq by expectation maximization)方法，文獻[12]提出的BitSeq(Bayesian inference of transcripts from sequencing data)方法，都采取了一定手段來解決讀段的多源映射問題。RNA-Seq數(shù)據(jù)的另外一個特點是讀段在參考序列上的分布有著自身的一些模式，呈現(xiàn)出不均勻性[13-16]。Cufflinks方法把位置偏差和序列偏差融入到計算中，模擬了讀段非均勻采樣的隨機性質(zhì)[9]。RSEM方法通過一個產(chǎn)生式概率模型模擬讀段產(chǎn)生過程的經(jīng)驗分布等信息，在一定程度上消除了均勻分布的假設(shè)[10]。

針對RNA-Seq讀段的異構(gòu)體多源映射與分布不均勻兩個難點，文獻[17]提出了LDASeq(Latent Dirichilet allocation for sequencing)模型來計算異構(gòu)體表達水平。LDASeq模型通過引入隱含變量，較好模擬了讀段的非均勻分布和多源映射現(xiàn)象。后續(xù)工作中發(fā)現(xiàn)，該模型存在以下缺點：首先忽視了兩類重要讀段，即接合區(qū)讀段和噪聲讀段對結(jié)果的重要影響。同時，對小外顯子沒有采取讀段數(shù)目歸一化且由于外顯子分割成探針而帶來外顯子尾部不足一個探針部分信息丟失，并且沒有考慮由于噪聲和未知轉(zhuǎn)錄本的存在對已知轉(zhuǎn)錄本表達值的推斷產(chǎn)生影響。針對LDASeq模型存在的這些問題，本文對LDASeq方法進行改進。改進的方法舍棄“偽探針”直接對外顯子進行長度歸一化，減小了探針引入所帶來的信息丟失，并且增加了對接合區(qū)讀段和噪聲讀段的處理。本文通過采用兩個真實數(shù)據(jù)集對所提出的改進方法在基因和異構(gòu)體表達水平計算精度上進行了驗證。

1 DASeqⅡ方法設(shè)計

1.1 LDASeq模型

基于RNA-Seq數(shù)據(jù)和文本數(shù)據(jù)在結(jié)構(gòu)上的相似性，文獻[17]提出了LDASeq模型。LDASeq是基于LDA(Latent Dirichlet allocation)模型[18]的一個3層貝葉斯網(wǎng)絡(luò)模型，如圖1所示。LDASeq模型引入固定長度的探針，將探針和單詞概念對應(yīng)，探針上讀段的匹配個數(shù)看作單詞在文檔中出現(xiàn)的頻數(shù)。單個通道下一個基因的N個探針上的讀段對應(yīng)一篇文檔，若干個異構(gòu)體看作文檔的隱含主題，α是一個超參數(shù)，模型中引入隱含變量θ，θ是一個向量，符合參數(shù)為α的狄利克雷分布，θ的各個分量反映了異構(gòu)體表達的相對強弱。模型中β是K×V矩陣，其中K表示該基因形成的剪切異構(gòu)體個數(shù)，V表示該基因含有的探針個數(shù)，第i(1≤i≤K)行第j(1≤j≤V)個分量表示第i個剪切異構(gòu)體是否含有第j個探針，因此每個基因?qū)?yīng)一個β矩陣。將β矩陣進行歸一化作為模型的輸入，文獻[17]中描述了β矩陣的產(chǎn)生過程。

LDASeq模型對每個基因模擬RNA-Seq數(shù)據(jù)的產(chǎn)生過程如下：(1) 根據(jù)狄利克雷分布為每個通道產(chǎn)生θl～Dirichlet(α)；(2) 對于通道l，根據(jù)多項分布產(chǎn)生異構(gòu)體isoformn。isoformn～Multinomial(θl)；(3) 在異構(gòu)體isoformn條件下，根據(jù)多項分布產(chǎn)生探針proben，proben～Multinomial(β)。

圖1 LDASeq圖模型表示Fig.1 Graphic model representation of LDASeq

文獻[17]對該模型使用變分最大期望算法(Expectation maximization，EM)來計算。模型求解得到參數(shù)θ(θ～Dirichlet(α))，每個分量代表著一個異構(gòu)體表達比重，然后將基因每個探針上映射的讀段計數(shù)按這個比重分配給對應(yīng)的異構(gòu)體。采用文獻[5]中的RPKM異構(gòu)體表達式計算每個異構(gòu)體表達值，將一個基因所對應(yīng)的所有異構(gòu)體表達值求和即可得到該基因表達值[9]。

1.2 LDASeqⅡ模型

針對LDASeq模型存在的一些缺點，本文對其進行了改進, 將改進后的方法記為LDASeqⅡ(Improment of latent Dirichilet allocation for sequencing data)，改進后模型如圖2所示。與原有模型相比：(1) 舍棄原模型探針概念，將外顯子與單詞對應(yīng)。為了減小不同外顯子長度對讀段數(shù)目的偏好，改進方法中將讀段數(shù)目按照外顯子長度歸一化，即單位外顯子長度上讀段個數(shù)作為單詞出現(xiàn)次數(shù)。(2) 為了處理由未知轉(zhuǎn)錄本和測序錯誤等產(chǎn)生的噪聲讀段和帶有重要剪切信息的接合區(qū)讀段，對每個基因所屬異構(gòu)體結(jié)構(gòu)進行擴展，增加偽異構(gòu)體和偽外顯子。

圖2 LDASeqⅡ圖模型表示Fig.2 Graphic model representation of LDASeqⅡ

1.2.1 偽外顯子的設(shè)計

圖3以基因G為例說明偽外顯子構(gòu)造過程。首先，擴充該基因的外顯子集合。對于基因G，原有外顯子集合記為F={e1，e2，e3，e4}，由圖3結(jié)構(gòu)關(guān)系可知，該基因發(fā)生了選擇性剪接形成兩個轉(zhuǎn)錄本。由轉(zhuǎn)錄本1可知，基因外顯子e1與外顯子e3之間發(fā)生選擇性剪接，因此在e1與e3之間構(gòu)造一個接合區(qū)記為e1-e3，作為偽外顯子加入集合F，這樣基因G外顯子集合F={e1，e2，e3，e4，e1-e3}。由轉(zhuǎn)錄本2可知，基因在外顯子e1與外顯子e2之間、外顯子e2與外顯子e4之間發(fā)生選擇性剪接，因此分別在e1與e2之間、e2與e4之間構(gòu)造接合區(qū)e1-e2，e2-e4，集合F中不含有元素e1-e2和e2-e4，因此將這兩個元素作為偽外顯子加入集合F，這樣集合F變?yōu)閧e1，e2，e3，e4，e1-e3，e1-e2，e2-e4}。

圖3 外顯子和讀段之間關(guān)系Fig.3 Relationship between exons and reads

接下來構(gòu)造新的轉(zhuǎn)錄本與外顯子映射關(guān)系。擴增后的基因外顯子集合F={e1，e2，e3，e4，e1-e3，e1-e2，e2-e4}。根據(jù)注釋信息可知，轉(zhuǎn)錄本1包含外顯子子集記為f1={e1，e3，e1-e3}，轉(zhuǎn)錄本2包含外顯子子集記為f2={e1，e2，e4，e1-e2，e2-e4}。這樣，改造后的基因G和轉(zhuǎn)錄本參考結(jié)構(gòu)如圖4所示。當統(tǒng)計外顯子讀段時，比如讀段1，完全落在外顯子e3上，那么將e3上讀段數(shù)目加1；而對于讀段2，落在外顯子e1與外顯子e2接合處，因此，將偽外顯子e1-e2上讀段數(shù)加1,其他讀段作類似統(tǒng)計。

圖4 增加偽外顯子、偽異構(gòu)體后的基因結(jié)構(gòu)圖Fig.4 Gene structure of adding ″pseudo-exon″ and ″pseudo-transcript″

1.2.2 偽轉(zhuǎn)錄本的引入

由于現(xiàn)有基因注釋信息的不全面，沒有包含那些尚未發(fā)現(xiàn)但已經(jīng)存在的的轉(zhuǎn)錄本，另外，RNA-Seq數(shù)據(jù)中還存在由于測序錯誤等產(chǎn)生的其他噪聲讀段。考慮到這兩類噪聲讀段對于已知轉(zhuǎn)錄本表達水平的影響，試圖在建模中加以矯正。通過改變每個基因所對應(yīng)的β的初始值，即對于每個基因，嘗試增加一條特殊的轉(zhuǎn)錄本參與優(yōu)化，將那些噪聲讀段看作這條特殊轉(zhuǎn)錄本所產(chǎn)生。對于圖3中的基因G，模型中嘗試引入了一條偽異構(gòu)體帶表這條基因尚未發(fā)現(xiàn)的剪切異構(gòu)體，圖4中即為基因G引入了偽異構(gòu)體。這條異構(gòu)體含有基因外顯子集合F中所有的外顯子。文章[17]描述了β矩陣產(chǎn)生規(guī)則。在對基因結(jié)構(gòu)進行擴增即增加了偽外顯子和偽轉(zhuǎn)錄本后，β產(chǎn)生規(guī)則不變。圖4所示擴增結(jié)構(gòu)后基因?qū)?yīng)的β矩陣初始值的產(chǎn)生過程

(1)

1.2.3 歸一化的改進

對于圖4中擴增結(jié)構(gòu)后的基因G′，該基因含有的7個外顯子對應(yīng)7個單詞。為了減小由于外顯子長度不同導(dǎo)致外顯子上讀段數(shù)目的偏好，將每個外顯子上讀段數(shù)目除以該外顯子長度，取單位長度上的讀段個數(shù)作為單詞出現(xiàn)頻數(shù)。對于基因G′原有的真實外顯子，即{e1，e2，e3，e4}，歸一化長度取真實長度-讀段長度，對于個別讀段小于讀段長度的小外顯子，將上面的讀段統(tǒng)計到相鄰的接合區(qū)，無需歸一化；對于新構(gòu)造的偽外顯子{e1-e3，e1-e2，e2-e4}，歸一化長度取讀段長度。對于雙末端讀段數(shù)據(jù)，由于讀段片段長度不定，實驗中取其平均長度。

2 實驗過程及結(jié)果

2.1 數(shù)據(jù)處理流程

由于增加了對接合區(qū)讀段、噪聲讀段的處理，改進后方法的的數(shù)據(jù)處理過程與LDASeq方法[17]有所不同,整個處理過程如圖5所示。

圖5 LDASeqⅡ處理流程Fig.5 Processing of LDASeqⅡ

處理過程為：第1步，用序列比對軟件將讀段比對到轉(zhuǎn)錄組上；第2步， 根據(jù)注釋信息，擴充基因外顯子集合、轉(zhuǎn)錄本集合；第3步，根據(jù)擴充后的基因結(jié)構(gòu)，寫出基因?qū)?yīng)的β矩陣并歸一化；第4步，統(tǒng)計基因各個外顯子以及偽外顯子上面的讀段數(shù)；第5步，將第3步和第4步的數(shù)據(jù)作為模型輸入，計算轉(zhuǎn)錄本和基因表達值。

2.2 實驗數(shù)據(jù)集

本文使用兩個真實數(shù)據(jù)集對LDASeqII的計算性能進行驗證，分別為人類大腦雙末端數(shù)據(jù)集[19]和人類乳腺癌數(shù)據(jù)集(Human breast cancer, HBC)[20]。人類大腦數(shù)據(jù)集來自美國食品藥品管理局聯(lián)合全球多家研究機構(gòu)舉辦的基因芯片質(zhì)量控制(Microarray quality control， MAQC)項目[19]。該數(shù)據(jù)集包含1 000個經(jīng)過定量反轉(zhuǎn)錄聚合酶鏈式反應(yīng)(Quantificational real-time polymerase chain reaction, qRT-PCR)實驗驗證的基因，這些結(jié)果可以作為基因表達水平計算結(jié)果的判別標準。人類乳腺癌數(shù)據(jù)集[20]包含兩個實驗條件，即乳腺癌細胞(Human breast cancer cell line, MCF-7)和正常乳腺細胞(Normal cell line, HME)。文獻[21]也采用了該數(shù)據(jù)集，并提供了經(jīng)過驗證的4個多異構(gòu)體基因共8個轉(zhuǎn)錄本在不同條件下的表達調(diào)控值。這些結(jié)果可以作為驗證轉(zhuǎn)錄本表達計算性能的一個標準。

2.3 基因水平上的驗證結(jié)果

在人類大腦雙末端數(shù)據(jù)集上，LDASeqⅡ共計算出740個多異構(gòu)體基因的表達值，并求出了表達值與qRT-PCR結(jié)果的相關(guān)系數(shù)，將結(jié)果與改進前的LDASeq， Cufflinks和RSEM這3種方法求出的相關(guān)系數(shù)作比較，各種方法的相關(guān)性如表1所示。從表中可以看出，LDASeqⅡ相比其他3種方法獲得了最高相關(guān)系數(shù)，基因表達水平獲得了較準確的計算結(jié)果。

表1 不同方法在MAQC數(shù)據(jù)集上的處理結(jié)果對比

2.4 轉(zhuǎn)錄本水平上的驗證結(jié)果

在人類乳腺癌數(shù)據(jù)集上，LDAseqⅡ, Cufflinks, RSEM和LDASeq這4種方法分別計算,經(jīng)過qRT-PCR實驗驗證的8個轉(zhuǎn)錄本在HME和MCF-7兩個條件下的表達水平和表達值變化方向調(diào)控值，將結(jié)果與qRT-PCR結(jié)果進行對比,如表2所示。

表2 人類乳腺癌數(shù)據(jù)集各方法計算結(jié)果

表2顯示了同一轉(zhuǎn)錄本在不同實驗條件、同一基因的不同轉(zhuǎn)錄本在同一實驗條件下表達值的變化方向，共16種表達調(diào)控關(guān)系。“+”表示在相應(yīng)的比較條件中轉(zhuǎn)錄本表達上調(diào)，“-”表示下調(diào)，括號中的數(shù)字表示調(diào)控的對數(shù)倍數(shù)。表2的最后兩行顯示了各種計算方法得到的調(diào)控方向與qRT-PCR結(jié)果不一致的比較次數(shù)，以及16個對數(shù)倍數(shù)與qRT-PCR結(jié)果的相關(guān)系數(shù)。調(diào)控方向反映模型在定性計算上與qRT-PCR結(jié)果的一致性，從表中可以看出，只有LDASeqⅡ計算出的16個調(diào)控方向和qRT-PCR計算出的方向完全一致，而Cufflinks， RSEM和LDASeq方法的調(diào)控方向錯誤個數(shù)分別為4，5和1。LDASeqⅡ方法在保持調(diào)控方向與qRT-PCR結(jié)果一致的情況下，相關(guān)系數(shù)達到0.733 8。由此可以看出，改進后的方法即LDAseqⅡ，在人類乳腺癌數(shù)據(jù)集上獲得了相比其他方法更為準確的轉(zhuǎn)錄本表達水平計算結(jié)果。

3 結(jié)束語

本文針對已有的LDASeq模型存在的一些缺點進行改進，提出了LDASeqⅡ方法來計算基因及轉(zhuǎn)錄本表達值。模型打破了讀段在參考序列上均勻分布這一常用假設(shè)，通過統(tǒng)計基因外顯子和接合區(qū)上讀段數(shù)目，利用轉(zhuǎn)錄本的結(jié)構(gòu)信息構(gòu)造初始化β矩陣對模型參數(shù)進行約束，優(yōu)化不同轉(zhuǎn)錄本下各個外顯子表達強弱，最終估計每個轉(zhuǎn)錄本的表達比重，較好地解決了轉(zhuǎn)錄本表達的計算問題。LDASeqⅡ方法在LDASeq的基礎(chǔ)上，改進了讀段數(shù)目歸一化方法，增加了對接合區(qū)讀段、噪聲讀段的處理，實驗表明，改進后的模型獲得了較優(yōu)的計算性能。相比其他大多數(shù)計算轉(zhuǎn)錄本表達值方法的創(chuàng)新之處在于模型中通過嘗試引入未知異構(gòu)體，減小了噪聲讀段對已知異構(gòu)體表達結(jié)果的影響，可以作為提高轉(zhuǎn)錄本表達研究性能的一個借鑒。本文在構(gòu)造特殊外顯子時，假設(shè)它含有基因的所有外顯子，在后續(xù)工作中，可以嘗試隨機分配外顯子或者根據(jù)基因已知的選擇性剪接信息，讓模型自學(xué)習(xí)未知的轉(zhuǎn)錄本結(jié)構(gòu)。另外，接合區(qū)讀段包含選擇性剪切的重要信息，通過對該區(qū)域讀段的處理，可以提高異構(gòu)體表達水平的計算準確性，本文方法對接合區(qū)讀段單獨處理，充分保留了該區(qū)域讀段帶有的選擇性剪接信息。對于接合區(qū)長度，實驗中采用一個讀段長度，也可以根據(jù)接合區(qū)兩端的外顯子長度按比例設(shè)定接合區(qū)的長度，在后續(xù)工作中會進一步驗證該方法的實際效果。

[1] Wang Z, Gerstein M, Snyder M.RNA-Seq: A revolutionary tool for transcriptomics [J].Nature Reviews Genetics, 2009, 10(1): 57-63.

[2] Denoeud F, Aury J M, Da Silva C, et al.Annotating genomes with massive scale RNA sequencing[J]. Genome Biol, 2008, 9(12):R175.

[3] Garber M, Grabherr M G, Guttman M, et al.Computational methods for transcriptome annotation and quantification using RNA-Seq[J].Nature Methods, 2011, 8(6): 469-477.

[4] Marguerat S, Bahler J.RNA-seq: From technology to biology[J].Cell Mol Life Sci, 2010, 67: 569-579.

[5] Mortazavi A, Williams B A, Mccue K, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq[J]. Nature Methods, 2008, 5(7): 621-628.

[6] Pan Q, Shai O, Lee L J, et al.Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing [J].Nature Genetices, 2008, 40(12): 1413-1415.

[7] Turro E, Su S Y, Goncalves, et al. Haplotype and isoform specific expression estimation using multi-mapping RNA-Seq reads [J]. Genome Biology, 2011, 12: R13.

[8] Jiang Hui, Wong Winghung. Statistical inferences for isoform expression in RNA-Seq [J]. Bioinformatics, 2009, 25(8): 1026-1032.

[9] Trapnell C, Williams B A, Pertea G, Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[J]. Nat Biotechnol, 2010(5): 511-515.

[10]Li B, Ruotti V, Stewart R M, et al.RNA-Seq gene expression estimation with read mapping uncertainty [J]. Bioinformatics, 2010, 26(4): 493-500.

[11]Li B, Dewey C N.RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome[J].BMC Bioinformatics, 2011, 12: 323.

[12]Glaus P, Honkela A, Rattray M.Identifying differentially expressed transcripts from RNA-Seq data with biological variation[J].Bioinformatics, 2012, 28(3): 1721-1728.

[13]Pepke S, Wold B, Mortazavi A.Computation for ChIP-Seq and RNA-Seq studies[J].Nature Methods Supplement, 2009, 6: S22-S32.

[14]Dohm J C, Lottaz C, Borodina T, Substantial biases in ultra-short read data sets from high-throughput DNA sequencing[J].Nucleic Acids Res, 2008(16): e105.

[15]Li J, Jiang H, Wong W H.Modeling non-uniformity inshort-read rates in RNA-Seq data[J].Genome Biol, 2010(5): R50.

[16]Hansen K D, Brenner S E, Dudoit S.Biases in illumina transcriptome sequencing caused by random hexamer priming[J]. Nucleic Acids Research, 2010, 38(12): e131.

[17]劉學(xué)軍，李蒙，張禮.一種針對RNA-Seq數(shù)據(jù)的基因異構(gòu)體表達水平計算方法[J].中國生物醫(yī)學(xué)工程學(xué)報,2013,7(4): 454-463.

Liu Xuejun， Li Meng， Zhang Li.A method of isoform expression calculation for RNA-Seq data[J].Chinese Journal of Biomedical Engineering,2013,7(4): 454-463.

[18]Blei D M, Ng A Y, Jordan M I.Latent Dirichlet allocation[J].Journal of Machine Learning Research, 2003, 3: 993-1022.

[19]Shi L, Reid L H, Jones W D, et al. The microarry quality control(MAQC) project shows inter- and intraplatform reproducibility of gene expression measurements[J].Nat Biotechnol, 2006, 24(9): 1151-1161.

[20]Wang E T, Sandberg R, Luo S J, et al.Alternative isoform regulation in human tissue transcriptomes[J].Nature, 2008, 456(7221): 470-476.

[21]Kim H, Bi Y T, Pal S, et al.IsoformEx: Isoform level gene expression estimation using weighted non-negative least squares from mRNA-Seq data[J].BMC Bioinformatics, 2011, 12: 305.

Improved Trancriptome Expression Analysis for RNA-Seq Data

Shi Xinxin, Liu Xuejun, Zhang Li

(College of Computer Science and Technology， Nanjing University of Aeronautics and Astronautics， Nanjing, 210016， China)

RNA-Seq(RNA-sequencing), based on high-throughput sequencing, is a new technique for transcriptome research.Considering the difficulties in the analysis of transcript expression using RNA-Seq data, an improved method, improvement of latent dirichlet allocation for sequencing data(LDASeqⅡ) is proposed to calculate the transcript expression.To deal with multi-mappings between reads and isoforms and non-uniform distribution of reads along reference, LDASeqⅡ utilizes the known gene-isoform annotation to constrain the hyperparameters and normalizes the read counts by exon length for each individual exon.By introducing ″pseudo-exon″ and ″pseudo-transcript″, the conjunction reads and noise reads gain proper treatments.LDASeqⅡ is validated using two real datasets on gene and transcript expression calculation and compared with latent dirichlet allocation for sequencing data(LDASeq) and other two popular methods Cufflinks and RNA-Seq by expectation maximization(RSEM). The results show that LDASeqⅡ obtains more accurate transcript and gene expression measurements than other approaches.

gene expression; RNA-Seq; transcript expression; multi-mapping; non-uniformity

國家自然科學(xué)基金(61170152)資助項目；中央高?；究蒲袠I(yè)務(wù)費專項(CXZZ11_0217)資助項目。

2014-05-30；

2014-06-25

TP391.9

A

石新新(1989-)，女，碩士研究生，研究方向：生物信息學(xué)，E-mail: shixinxin61@126.com。

劉學(xué)軍(1976-)，女，博士,教授，研究方向：機器學(xué)習(xí)與生物信息學(xué)。

張禮(1985-)，男，博士，研究方向：生物信息學(xué)。