歐書華，劉學(xué)軍，張　禮

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

ISSN 1673-9418 CODEN JKYTA8

Journal of Frontiers of Computer Science and Technology 1673-9418/2016/10(03)-0381-08

?

基于平滑LDA的RNA-Seq數(shù)據(jù)表達(dá)分析研究*

歐書華+，劉學(xué)軍，張禮

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

ISSN 1673-9418 CODEN JKYTA8

Journal of Frontiers of Computer Science and Technology 1673-9418/2016/10(03)-0381-08

E-mail: fcst@vip.163.com

http://www.ceaj.org

Tel: +86-10-89056056

* The National Natural Science Foundation of China under Grant No. 61170152 (國家自然科學(xué)基金); the Qinglan Project of Jiangsu Province (江蘇省青藍(lán)工程) ; the Fundamental Research Funds for the Central Universities of China under Grant No. CXZZ11_0217(中央高?；究蒲袠I(yè)務(wù)費專項資金).

Received 2015-04,Accepted 2015-06.

CNKI網(wǎng)絡(luò)優(yōu)先出版: 2015-06-19, http://www.cnki.net/kcms/detail/11.5602.TP.20150619.1044.001.html

摘要：RNA-Seq是目前轉(zhuǎn)錄組研究的一種重要技術(shù)，針對RNA-Seq數(shù)據(jù)分析中讀段的多源映射，參考序列分布的不均勻性，一些轉(zhuǎn)錄本中外顯子分布稀疏以及跨結(jié)合區(qū)讀段處理問題，提出了一個新的轉(zhuǎn)錄組表達(dá)研究book=382,ebook=86模型sLDASeq。該模型根據(jù)基因中轉(zhuǎn)錄本注釋信息對模型參數(shù)進(jìn)行約束，對跨結(jié)合區(qū)的讀段按長度分配處理，解決了讀段非均勻分布和跨結(jié)合區(qū)問題；在模型中增加一個超參數(shù)，從而解決了外顯子的稀疏問題。將該模型應(yīng)用到3個真實的數(shù)據(jù)集上，并與其他主流方法進(jìn)行比較，結(jié)果表明該模型獲得了較為準(zhǔn)確的基因以及轉(zhuǎn)錄本表達(dá)水平計算結(jié)果。

關(guān)鍵詞：RNA-Seq；基因轉(zhuǎn)錄本表達(dá)水平；平滑LDA；結(jié)合區(qū)；多源映射；非均勻性

1　引言

RNA-Seq是近年來發(fā)展迅速的深度測序技術(shù)[1-3]，該技術(shù)通過記錄每個基因中讀段的表達(dá)頻率來量化轉(zhuǎn)錄組，相比其他方法具有信噪比高、靈敏度高、分辨率高、背景噪聲低、所需樣本少等特點[4]，并提供了尋找新的基因、轉(zhuǎn)錄本和病變的組織特異基因的方法[5-7]。RNA-Seq數(shù)據(jù)處理通常分為3個步驟:首先，將讀段映射到參考基因組或轉(zhuǎn)錄組上；其次，利用比對上的讀段重建參考基因組或轉(zhuǎn)錄組；最后，計算比對到基因組或轉(zhuǎn)錄組上的讀段數(shù)，獲得基因或轉(zhuǎn)錄本的表達(dá)水平。

通過序列對比將RNA-Seq原始讀段映射到基因組或轉(zhuǎn)錄組后，若以映射到基因上讀段數(shù)衡量基因的表達(dá)水平是不合理的，為此Mortazavi等人提出RPKM（reads per kilo bases per million reads）方法估計基因的表達(dá)水平[8]。RPKM方法是用讀段的覆蓋度來刻畫基因的表達(dá)水平，僅在讀段是均勻分布和不含異構(gòu)體基因上基因的表達(dá)水平估計較為準(zhǔn)確，而對含選擇性剪切異構(gòu)體的基因整體表達(dá)水平上低于真實值。針對讀段在多異構(gòu)體基因上多源映射問題，Bullard等人提出UI（union-intersection）基因來估計基因整體表達(dá)水平[9]。Trapnell等人提出采用基因的異構(gòu)體表達(dá)水平之和來估計基因表達(dá)水平的方法[10]。Jiang等人利用泊松分布對讀段在外顯子上分布建模，提出極大似然估計求得異構(gòu)體表達(dá)水平估計值[11]。針對讀段在參考序列不均勻分布問題，Li等人提出了通過利用經(jīng)驗分布和讀段質(zhì)量匹配信息來消除不均勻分布影響的RSEM（RNA-Seq by expectation maximization）方法[12-13]。Zhang等人提出了利用全局和局部偏差曲線描述讀段非均勻分布特性的N-URD （non-uniformreaddistribution）方法[14-15]。此外還有Glaus提出的BitSeq方法[16]、模擬讀段隨機采用特性的Cufflinks方法[17]以及Seqguio方法[8]均對讀段的非均勻分布采取了一些措施。

針對RNA-Seq讀段數(shù)據(jù)在參考序列上的多源映射和分布的不均勻性，基于LDA（latent Dirichlet allocation）模型[18]處理的文本數(shù)據(jù)和RNA-Seq數(shù)據(jù)的相似性[19]，在先前的工作中提出了基于LDA模型的LDASeq方法用于計算基因和轉(zhuǎn)錄本的表達(dá)水平。LDASeq模型通過引入隱含變量模擬異構(gòu)體和對LDA模型中β矩陣的優(yōu)化，較好地解決了讀段的多源映射和非均勻分布問題，但忽略了外顯子長度對讀段的影響，沒有處理讀段落在外顯子之間的結(jié)合區(qū)的情況，同時也沒有考慮多外顯子基因中異構(gòu)體的外顯子稀疏問題。針對以上問題本文提出sLDASeq （smoothed LDASeq）模型，對外顯子長度歸一化和結(jié)合區(qū)的讀段進(jìn)行處理，增加了一個超參數(shù)來解決模型中的稀疏問題。本文通過3個真實數(shù)據(jù)集對提出的模型進(jìn)行了驗證，并與其他主流方法進(jìn)行了對比。

2　方法

2.1平滑LDA

平滑LDA模型[18]是文本數(shù)據(jù)分析中，針對LDA模型中不在訓(xùn)練集中的詞項，多項式參數(shù)β會將其概率設(shè)置為0而帶來的稀疏問題的改進(jìn)。如圖1所示，α為隱含主題概率的先驗概率分布參數(shù)，θd為文檔d下的主題概率分布, zn,d為第n個單詞在文檔d下生成的主題，wn,d為主題產(chǎn)生的單詞，β為主題和單詞的K×V概率矩陣，βij表示第i個隱含主題生成第j個單詞的概率，η為生成β的先驗概率參數(shù)。

2.2 sLDASeq模型

考慮到RNA-Seq數(shù)據(jù)和文本數(shù)據(jù)結(jié)構(gòu)上的相似性，LDASeq模型[20]是基于LDA提出的RNA-Seq數(shù)據(jù)分析方法，它是一個三層貝葉斯結(jié)構(gòu)模型。該模型將固定長度的探針對應(yīng)文檔的單詞，探針上的讀段數(shù)對應(yīng)單詞出現(xiàn)的頻率，將單個通道中一個基因上的探針及讀段個數(shù)看作一篇文檔，M篇文檔形成一個文檔集，基因的異構(gòu)體則對應(yīng)文檔的隱含主題，模型中的隱含θ變量反映基因中異構(gòu)體的比重。LDASeq模型中β矩陣初值是由注釋文件中每個基因的異構(gòu)體所含外顯子的信息得到的：若基因第i個異構(gòu)體含有第j個外顯子，那么矩陣第i行j列初值為1，反之為0，再將矩陣歸一化得到β矩陣的初值，并在計算過程中只優(yōu)化β矩陣中的非零元素。因此在基因的外顯子多而異構(gòu)體的外顯子少的情況下，按LDASeq模型處理較容易造成β矩陣稀疏問題。

Fig.1　 Model representation of smoothed LDA圖1　平滑LDA模型

本文針對于LDASeq中一些缺點提出了基于平滑LDA的sLDASeq模型?，F(xiàn)用圖2所示的基因結(jié)構(gòu)圖說明RNA-Seq數(shù)據(jù)與文本數(shù)據(jù)的相似性，在sLDASeq中舍棄了LDASeq模型中的探針，將基因的外顯子和LDA中的單詞對應(yīng)，外顯子上讀段個數(shù)對應(yīng)單詞的頻數(shù)，單通道中一個基因的外顯子及讀段數(shù)對應(yīng)一篇文檔，基因的剪切異構(gòu)體則對應(yīng)LDA模型中的隱含主題。為消除外顯子長度對讀段數(shù)目的影響，sLDASeq將外顯子上讀段數(shù)按外顯子長度歸一化，即將每個外顯子上讀段數(shù)除以該外顯子的長度，并將單位外顯子長度上的讀段個數(shù)作為單詞出現(xiàn)的頻數(shù)。假設(shè)圖2中基因的外顯子1和4上讀段數(shù)分別為1和2，外顯子長度分別為10bp、20bp，那么歸一化后基因上外顯子1和4上單位長度讀段數(shù)均為0.1。

Fig.2　 Gene structure圖2　基因結(jié)構(gòu)圖

LDASeq模型對圖3跨結(jié)合區(qū)的讀段處理是將讀段歸為前一個或后面的外顯子，導(dǎo)致外顯子上的讀段數(shù)不太準(zhǔn)確，為了更為準(zhǔn)確地計算外顯子上的讀段數(shù)，sLDASeq模型將跨結(jié)合區(qū)的讀段按其在各個外顯子上的長度分配，然后歸一化。假設(shè)圖3中所示的讀段1長度為20bp，落在外顯子1、2上的長度分別為5bp、15bp，那么就將外顯子1上的讀段數(shù)計為0.25，外顯子2上的讀段數(shù)計為0.75，這樣在最終每個外顯子的讀段計數(shù)中保留了跨結(jié)合區(qū)讀段的統(tǒng)計特性，有助于模型優(yōu)化時獲得相應(yīng)異構(gòu)體相對表達(dá)值。

Fig.3　 reads in exon-junction圖3　結(jié)合區(qū)的讀段處理

為解決LDASeq模型中β矩陣稀疏問題，該模型增加了一個超參數(shù)η生成異構(gòu)體中外顯子（exon）上讀段出現(xiàn)的概率矩陣β，α為生成異構(gòu)體（isoform）概率的狄利克雷分布參數(shù)，θ表示各個異構(gòu)體分布強弱關(guān)系概率分布。該模型下RNA-Seq數(shù)據(jù)生成流程如下。

（1）根據(jù)狄利克雷分布（Dirichlet）生成每個通道的異構(gòu)體概率分布θd：

θd~Dirichlet(α)

（2）對于通道d生成的θd，按多項式分布（multinomial）選擇第n個讀段的異構(gòu)體isoformn:

isoformn~multinomial(θd)

（3）根據(jù)超參數(shù)η生成第i個異構(gòu)體的外顯子上讀段出現(xiàn)的概率向量βi：

βi~Dirichlet(η)

重復(fù)K次得到一個基因上K個異構(gòu)體和外顯子上讀段出現(xiàn)的概率矩陣β。

（4）最后以生成的異構(gòu)體和β為條件，通過多項式分布概率選擇外顯子：

p(exonn|isoformn,β)

在給定超參數(shù)α和η，模型關(guān)于隱含變量θ、isoform、exon、β的聯(lián)合概率為：

每個通道均相互獨立，含有M個通道基因的似然函數(shù)為:

由于β和θ的耦合性使得上式直接求解很困難，通過變分推理將p(θ,isoform,β|exon,α,η)近似等于q(θ,isoform,β|γ,?,λ)，再根據(jù)變分EM算法求解模型：

E步給定α和η，根據(jù)變分推理優(yōu)化變分參數(shù)γ,φ,λ。

M步利用E步得到的γ,φ,λ，計算極大化下界p(D|α,η)的超參數(shù)α和η的值。

重復(fù)上述步驟，直到對數(shù)似然函數(shù)下界收斂。

根據(jù)文獻(xiàn)[20]，模型優(yōu)化出的參數(shù)θ反映了基因中異構(gòu)體表達(dá)的強弱程度，且θ服從狄利克雷分布θ~Dirichlet(α)，根據(jù)狄利克雷分布性質(zhì)：

模型優(yōu)化出的歸一化后超參數(shù)α代表了異構(gòu)體的比重，將映射到基因外顯子中的讀段數(shù)按這個比例分配給異構(gòu)體，采用下式中的FPKM（reads per kilobase of exon model per million mapped reads）公式計算異構(gòu)體的表達(dá)水平：

其中，n表示異構(gòu)體上外顯子的總讀段數(shù)；N表示所有基因上外顯子的總讀段數(shù)；l表示異構(gòu)體的長度?；虮磉_(dá)水平的計算則是將基因中所有異構(gòu)體的表達(dá)水平相加即可。

2.3 sLDASeq處理流程

sLDASeq采用bowtie2[21]將讀段比對到轉(zhuǎn)錄組參考序列。圖4給出了sLDASeq模型處理RNA-Seq數(shù)據(jù)流程。第一步，根據(jù)參考序列用bowtie2將讀段比對到轉(zhuǎn)錄組序列上，比對的結(jié)果包含讀段在轉(zhuǎn)錄組上的相對位置。第二步，根據(jù)讀段的位置和注釋信息得到其在基因上的位置，從而獲得基因上各個外顯子映射的讀段個數(shù)，其中對于跨結(jié)合區(qū)的讀段單獨處理。第三步，應(yīng)用sLDASeq模型將第二步得到的結(jié)果作為輸入計算轉(zhuǎn)錄本和基因的表達(dá)值FPKM。

Fig.4　 Workflow of sLDASeq approach圖4　 sLDASeq處理流程圖

3　實驗及討論

本文使用了3個真實的數(shù)據(jù)集，分別是基因芯片質(zhì)量控制（micro array quality control，MAQC）雙末端和單末端的數(shù)據(jù)集[22]、人類乳腺癌數(shù)據(jù)集（humanbreast cancer，HBC）[23]以及吸煙與肺癌關(guān)系（smoke exposure an lung cancer pathogenesis，SELC）數(shù)據(jù)集[24]。實驗分別采用以上3個數(shù)據(jù)集對sLDASeq模型的性能進(jìn)行驗證，并將實驗結(jié)果和LDASeq模型以及另外兩個主流方法Cufflinks[17]和RSEM進(jìn)行對比。

3.1基因水平上驗證

3.1.1 MAQC數(shù)據(jù)集

MAQC是美國食品藥品監(jiān)督管理局（Food and Drug Administration，F(xiàn)DA）領(lǐng)導(dǎo)并聯(lián)合多國研究機構(gòu)實施的項目，該項目得到生物芯片數(shù)據(jù)具有平臺內(nèi)和平臺間的重現(xiàn)性的結(jié)論。MAQC包含單末端和雙末端數(shù)據(jù)集，對其中約1 000個基因進(jìn)行qRT-PCR （quantitative real time PCR）[8,10-11]驗證用以作為判別其他方法優(yōu)劣的依據(jù)。該項目中包含兩個條件，分別是通用人類參考（universal human reference，UHR）和人類大腦參考（human brain reference，HBR）。

本文實驗主要采用的是UHR條件下雙末端數(shù)據(jù)集以及HBR和UHR條件下的單末端數(shù)據(jù)集。sLDASeq模型分別計算了經(jīng)過qRT-PCR驗證的約1 000基因中，在UHR條件下雙末端數(shù)據(jù)集、單末端數(shù)據(jù)集和HBR條件下的單末端數(shù)據(jù)集中被檢測到并匹配到注釋文件中的多異構(gòu)體的711、728、730個基因的表達(dá)值與經(jīng)qRT-PCR驗證的基因表達(dá)值的相關(guān)系數(shù)，即反映變量之間相關(guān)關(guān)系密切程度的統(tǒng)計指標(biāo)，其值越接近于1表明該方法與qRT-PCR驗證的結(jié)果吻合度越高，則越準(zhǔn)確。并將結(jié)果和Cufflinks、RSEM及LDASeq方法得到的相關(guān)系數(shù)進(jìn)行比較，由表1可見，在MAQC雙末端和單末端數(shù)據(jù)集上sLDASeq獲得了比其他方法更高的相關(guān)系數(shù)。結(jié)果表明sLDASeq相比其他方法在基因表達(dá)水平計算上取得了較為準(zhǔn)確的結(jié)果。

Table 1　 Results of various methods for MAQC dataset表1　不同方法在MAQC數(shù)據(jù)集上結(jié)果

3.1.2 SELC數(shù)據(jù)集

SELC數(shù)據(jù)集分別在NS（healthy never smoker）、S（current smoker）、C（smokers with lung cancer）、NC （smokers without lung cancer）4個條件下對樣本RNA數(shù)據(jù)進(jìn)行測序。本文在該數(shù)據(jù)集上分別用Cufflinks、RSEM、LDASeq、sLDASeq方法計算“S vs NS”和“C vs NC”條件分組下的7個基因LFC（logged foldchange）值，即將基因在兩個條件下求得的表達(dá)值相除再取以2為底的對數(shù)。將結(jié)果與qRT-PCR的LFC值進(jìn)行比較，實驗結(jié)果如表2所示。

Table 2　 Results of various methods for SELC dataset表2　 SELC數(shù)據(jù)集上各種方法結(jié)果

表2中第3列表示qRT-PCR在第1列兩個條件分組下獲得的LFC值，后5列分別是基因在條件“S vs NS”和“C vs NC”下用Cufflinks、RSEM、LDASeq、sLDASeq方法獲得的LFC值。括號里的數(shù)值是4個方法得到的LFC值與qRT-PCR獲得的LFC值的誤差率，即將每種方法獲得的LFC值與qRT-PCR的LFC值相減取絕對值，再除以qRT-PCR的LFC值。表中最后一行是每個方法中7個基因的誤差率取均值得到的平均誤差率。結(jié)果表明sLDASeq在該數(shù)據(jù)集上相比其他方法平均誤差率最低，在基因表達(dá)水平計算中取得了較為準(zhǔn)確的結(jié)果。

3.2轉(zhuǎn)錄本水平上驗證

實驗采用了人類乳腺癌數(shù)據(jù)集進(jìn)行轉(zhuǎn)錄本水平上的驗證。本文分別用Cufflinks、RSEM、LDASeq、sLDASeq方法計算乳腺癌細(xì)胞（cancer cell line，MCF-7）和正常乳腺細(xì)胞（breast normal cell line，HME）兩個條件下的4個基因8個轉(zhuǎn)錄本LFC值和調(diào)控方向變化，并將結(jié)果和qRT-PCR實驗LFC值比較得到誤差率，實驗結(jié)果如表3所示。

Table 3　 Results of various methods for human breast cancer dataset表3　人類乳腺癌數(shù)據(jù)集上各種方法計算結(jié)果

表3顯示了同一轉(zhuǎn)錄本在兩個不同條件中以及同一基因在一個條件下兩個不同轉(zhuǎn)錄本的調(diào)控方向和誤差率。表中的“+”表示轉(zhuǎn)錄本在不同實驗下方向上調(diào)（up-regulation），“-”表示下調(diào)（down-regulation）。括號中的數(shù)字表示其他幾種方法獲得的LFC值與qRT-PCR的LFC值比較獲得的誤差率。表中顯示Cufflinks、RSEM、LDASeq、sLDASeq調(diào)控方向錯的個數(shù)分別為4、5、1、1，且sLDASeq平均誤差率最小。結(jié)果表明相比其他幾個方法，sLDASeq在轉(zhuǎn)錄本表達(dá)水平計算中取得了較為準(zhǔn)確的結(jié)果。

4　結(jié)論

本文提出了基于平滑LDA的sLDASeq模型來計算基因和轉(zhuǎn)錄本的表達(dá)水平。sLDASeq模型利用已知的注釋文件對模型參數(shù)進(jìn)行約束，增加了對于跨外顯子讀段的處理，并在原LDA模型的基礎(chǔ)上增加了一個超參數(shù)生成轉(zhuǎn)錄本和外顯子上讀段出現(xiàn)的概率矩陣，解決了原模型中的稀疏問題。將所設(shè)計的sLDASeq模型應(yīng)用于基因芯片質(zhì)量控制（MAQC）、人類乳腺癌（HBC）、吸煙與肺癌關(guān)系（SELC）3個數(shù)據(jù)集上進(jìn)行驗證，并與先前提出的LDASeq模型、Cufflinks、RSEM的結(jié)果進(jìn)行對比。在MAQC數(shù)據(jù)集上本文模型獲得的基因表達(dá)值與qRT-PCR實驗結(jié)果的相關(guān)系數(shù)均比其他模型要高；在經(jīng)過qRT-PCR驗證的SELC和HBC數(shù)據(jù)集上，sLDASeq模型與其他方法相比，與qRT-PCR結(jié)果更為接近。由此得出，sLDASeq模型獲得了較為準(zhǔn)確的基因以及異構(gòu)體上表達(dá)水平。

本文所用的實驗數(shù)據(jù)均是已知基因的轉(zhuǎn)錄本注釋信息，但生物基因的轉(zhuǎn)錄本注釋信息并不是全部已知的，在后續(xù)的工作中，可以構(gòu)造多種未知的轉(zhuǎn)錄本的組合，利用模型學(xué)習(xí)出新的轉(zhuǎn)錄本。另外對于結(jié)合區(qū)的讀段處理，本文是將這些跨外顯子的讀段按長度比例分配到各個外顯子上，但是結(jié)合區(qū)的讀段含有選擇性剪切的信息，若保留這些讀段的選擇性剪切的信息對其單獨處理的效果，將在后續(xù)的工作中進(jìn)一步驗證。

References:

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

[2] Sultan M, Amstislavskiy V, Risch T. Influence of RNA extraction methods and library selection schemes on RNA-seq data[J]. BMC Genomics, 2014, 15: 675-688.

[3] Robert A W, Philippa A S, Catherine M M. RNA Seq analysis of the Eimeria tenella gametocyte transcriptome reveals clues about the molecular basis for sexual reproduction and oocyst biogenesis[J]. BMC Genomics, 2015, 16: 94-114.

[4] Wang Xi, Wang Xiaowo, Wang Likun. A review on the processing and analysis of next-generation RNA-seq data[J]. Progress in Biochemistry and Biophysics, 2010, 37(8): 837-846.

[5] Xiao Shengjian, Zhang Chi, Zou Quan, et al. TiSGeD: a database for tissue-specific genes[J]. Bioinformatics, 2010, 26 (9): 1273-1275.

[6] Pan Jianbo, Hu Shichang, Shi Dan, et al. PaGenBase: a pattern gene database for the global and dynamic understanding of gene function[J]. PLoS ONE, 2013, 8(12): e80747.

[7] Pan Jianbo, Hu Shichang, Wang Hao, et al. PaGeFinder: quantitative identification of spatiotemporal pattern genes[J]. Bioinformatics, 2012, 28(11): 1544-1545.

[8] 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.

[9] Bullard J H, Purdom E, Hansen K D, et al. Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments[J]. BMC Bioinformatics, 2010, 11: 94-107.

[10] Trapnell C, Williams B A, Pertea G. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[J]. Nature Biotechnology, 2011, 28(5): 511-515.

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

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

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

[14] Wu Zhengpeng, Wang Xi, Zhang Xuegong. Using non-uniform read distribution models to improve isoform expression inference in RNA-Seq[J]. Bioinformatics, 2011, 27(4): 502-508.

[15] Ma Xinyun, Zhang Xuegong. NURD an implementation of a new method to estimate isoform expression from non-uniform RNA-seq data[J]. BMC Bioinformatics, 2013, 14: 220-227.

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

[17] Roberts A, Trapnell C, Donaghey J, et al. Improving RNASeq expression estimates by correcting for fragment bias[J]. Genome Biology, 2011, 12(3): 22-36.

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

[19] Zhao Zhilan, Liu Xuejun, Zhang Li. RNA-Seq data analysis based on probabilistic model[C]//The Monograph of Academic Annual Conference Held by Chinese Society of Biomedical Engineering (CD). Wuhan: Chinese Society of Biomedical Engineering, 2011.

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

[21] Langmead B, Salzberg S L. Fast gapped-read alignment with Bowtie 2[J]. Nature Methods, 2012, 9(4): 357-359.

[22] Consortium M. The microarray quality control (MAQC) project shows inter- and intraplatform reproducibility ofgene expression measurements[J]. Nature Biotechnology, 2006, 24(9): 1151-1161.

[23] Kim H, Bi Yingtao, 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-314.

[24] Beane J, Vick J, Schembri F, et al. Characterizing the impact of smoking and lung cancer on the airway transcriptome using RNA-Seq[J]. Cancer Prevention, 2011, 4(6): 803-817.

附中文參考文獻(xiàn)：

[4]王曦,汪小我,王立坤.新一代高通量RNA測序數(shù)據(jù)的處理與分析[J].生物化學(xué)與生物物理進(jìn)展, 2010, 37(8): 837-846.

[19]趙志蘭,劉學(xué)軍,張禮.一種基于概率模型的RNA-Seq數(shù)據(jù)分析方法[C]//中國生物醫(yī)學(xué)工程聯(lián)合學(xué)術(shù)年會論文集(光盤版).武漢:中國生物醫(yī)學(xué)工程學(xué)會, 2011.

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

OU Shuhua was born in 1991. She is an M.S. candidate at College of Computer Science and Technology, Nanjing University of Aeronautics and Astronautics. Her research interests include bioinformatics and machine learning, etc.歐書華（1991—），女，南京航空航天大學(xué)計算機科學(xué)與技術(shù)學(xué)院碩士研究生，主要研究領(lǐng)域為生物信息學(xué)，機器學(xué)習(xí)等。

LIU Xuejun was born in 1976. She received the Ph.D. degree in computer science from the University of Manchester in 2006. Now she is a professor at College of Computer Science and Technology, Nanjing University of Aeronautics and Astronautics. Her research interests include bioinformatics and machine learning, etc.劉學(xué)軍（1976—），女，2006年于英國曼徹斯特大學(xué)獲得博士學(xué)位，現(xiàn)為南京航空航天大學(xué)計算機科學(xué)與技術(shù)學(xué)院教授，主要研究領(lǐng)域為生物信息學(xué)，機器學(xué)習(xí)等。

ZHANG Li was born in 1985. He received the M.S. degree in computer applications from Nanjing University of Aeronautics and Astronautics in 2010. Now he is a Ph.D. candidate at Nanjing University of Aeronautics and Astronautics. His research interests include bioinformatics and machine learning, etc.張禮（1985—），男，2010年于南京航空航天大學(xué)計算機應(yīng)用專業(yè)獲得碩士學(xué)位，現(xiàn)為南京航空航天大學(xué)博士研究生，主要研究領(lǐng)域為生物信息學(xué)，機器學(xué)習(xí)等。

RNA-Seq Data Expression Analysis Based on Smoothed LDA?

OU Shuhua+, LIU Xuejun, ZHANG Li
College of Computer Science and Technology, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
+ Corresponding author: E-mail: oshuhua@163.com

OU Shuhua, LIU Xuejun, ZHANG Li. RNA-Seq data expression analysis based on smoothed LDA. Journal of Frontiers of Computer Science and Technology, 2016, 10(3): 381-388.

Abstract:RNA-Seq is an important technique for transcriptome research. Considering the multi-mappings between reads and isoforms, non-uniform distribution of reads along the reference sequence, conjunction reads and the sparsity caused by the large exon size, this paper proposes a new method, sLDASeq, to calculate the gene and transcript expression. To solve the problems of multi-mappings, non-uniform distribution of reads and conjunction reads, the model utilizes the known gene-isoform annotation to constrain the hyper-parameters and allocate the read counts according to exon length. By adding a hyper-parameter, the model solves the problem of sparsity in the exons. sLDASeq is validated by using three real datasets on the gene and transcript expression calculation and compared with LDASeq and other popular methods. Results show that sLDASeq obtains more accurate transcript and gene expression measurements than other methods.

Key words:RNA-Seq; gene and transcript expression; smoothed LDA; exon-junction; multi-mapping; non-uniformity

doi:10.3778/j.issn.1673-9418.1505048

文獻(xiàn)標(biāo)志碼：A

中圖分類號：TP391