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

南京航空航天大學(xué)計(jì)算機(jī)科學(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é)計(jì)算機(jī)科學(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 (國(guó)家自然科學(xué)基金); the Qinglan Project of Jiangsu Province (江蘇省青藍(lán)工程) ; the Fundamental Research Funds for the Central Universities of China under Grant No. CXZZ11_0217(中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金).

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

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

1　引言

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

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

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

2　方法

2.1平滑LDA

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

2.2 sLDASeq模型

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

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

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

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

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

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

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

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

θd~Dirichlet(α)

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

isoformn~multinomial(θd)

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

βi~Dirichlet(η)

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

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

p(exonn|isoformn,β)

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

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

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

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

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

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

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

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

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

2.3 sLDASeq處理流程

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

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

3　實(shí)驗(yàn)及討論

本文使用了3個(gè)真實(shí)的數(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]。實(shí)驗(yàn)分別采用以上3個(gè)數(shù)據(jù)集對(duì)sLDASeq模型的性能進(jìn)行驗(yàn)證，并將實(shí)驗(yàn)結(jié)果和LDASeq模型以及另外兩個(gè)主流方法Cufflinks[17]和RSEM進(jìn)行對(duì)比。

3.1基因水平上驗(yàn)證

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

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

本文實(shí)驗(yàn)主要采用的是UHR條件下雙末端數(shù)據(jù)集以及HBR和UHR條件下的單末端數(shù)據(jù)集。sLDASeq模型分別計(jì)算了經(jīng)過(guò)qRT-PCR驗(yàn)證的約1 000基因中，在UHR條件下雙末端數(shù)據(jù)集、單末端數(shù)據(jù)集和HBR條件下的單末端數(shù)據(jù)集中被檢測(cè)到并匹配到注釋文件中的多異構(gòu)體的711、728、730個(gè)基因的表達(dá)值與經(jīng)qRT-PCR驗(yàn)證的基因表達(dá)值的相關(guān)系數(shù)，即反映變量之間相關(guān)關(guān)系密切程度的統(tǒng)計(jì)指標(biāo)，其值越接近于1表明該方法與qRT-PCR驗(yàn)證的結(jié)果吻合度越高，則越準(zhǔn)確。并將結(jié)果和Cufflinks、RSEM及LDASeq方法得到的相關(guān)系數(shù)進(jìn)行比較，由表1可見(jiàn)，在MAQC雙末端和單末端數(shù)據(jù)集上sLDASeq獲得了比其他方法更高的相關(guān)系數(shù)。結(jié)果表明sLDASeq相比其他方法在基因表達(dá)水平計(jì)算上取得了較為準(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個(gè)條件下對(duì)樣本RNA數(shù)據(jù)進(jìn)行測(cè)序。本文在該數(shù)據(jù)集上分別用Cufflinks、RSEM、LDASeq、sLDASeq方法計(jì)算“S vs NS”和“C vs NC”條件分組下的7個(gè)基因LFC（logged foldchange）值，即將基因在兩個(gè)條件下求得的表達(dá)值相除再取以2為底的對(duì)數(shù)。將結(jié)果與qRT-PCR的LFC值進(jìn)行比較，實(shí)驗(yàn)結(jié)果如表2所示。

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

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

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

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

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

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

4　結(jié)論

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

本文所用的實(shí)驗(yàn)數(shù)據(jù)均是已知基因的轉(zhuǎn)錄本注釋信息，但生物基因的轉(zhuǎn)錄本注釋信息并不是全部已知的，在后續(xù)的工作中，可以構(gòu)造多種未知的轉(zhuǎn)錄本的組合，利用模型學(xué)習(xí)出新的轉(zhuǎn)錄本。另外對(duì)于結(jié)合區(qū)的讀段處理，本文是將這些跨外顯子的讀段按長(zhǎng)度比例分配到各個(gè)外顯子上，但是結(jié)合區(qū)的讀段含有選擇性剪切的信息，若保留這些讀段的選擇性剪切的信息對(duì)其單獨(dú)處理的效果，將在后續(xù)的工作中進(jìn)一步驗(yà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測(cè)序數(shù)據(jù)的處理與分析[J].生物化學(xué)與生物物理進(jìn)展, 2010, 37(8): 837-846.

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

[20]劉學(xué)軍,李蒙,張禮.一種針對(duì)RNA-Seq數(shù)據(jù)的基因異構(gòu)體表達(dá)水平計(jì)算方法[J].中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào), 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é)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院碩士研究生，主要研究領(lǐng)域?yàn)樯镄畔W(xué)，機(jī)器學(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年于英國(guó)曼徹斯特大學(xué)獲得博士學(xué)位，現(xiàn)為南京航空航天大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院教授，主要研究領(lǐng)域?yàn)樯镄畔W(xué)，機(jī)器學(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é)計(jì)算機(jī)應(yīng)用專業(yè)獲得碩士學(xué)位，現(xiàn)為南京航空航天大學(xué)博士研究生，主要研究領(lǐng)域?yàn)樯镄畔W(xué)，機(jī)器學(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

中圖分類號(hào)：TP391