基于多樣本RNA-Seq數(shù)據(jù)的表達(dá)水平估計方法*

張　禮，劉學(xué)軍，陳松燦

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

基于多樣本RNA-Seq數(shù)據(jù)的表達(dá)水平估計方法*

張禮+，劉學(xué)軍，陳松燦

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

ZHANG Li,LIU Xuejun,CHEN Songcan.Novel method to estimate expression level based on multi-sample RNA-Seq data.Journal of Frontiers of Computer Science and Technology,2016,10(2)：210-219.

隨著下一代高通量DNA測序的快速發(fā)展，RNA-Seq測序已成為轉(zhuǎn)錄組學(xué)分析的標(biāo)準(zhǔn)技術(shù)。在處理多樣本RNA-Seq數(shù)據(jù)時，現(xiàn)有表達(dá)水平估計方法通?；趩蝹€樣本逐個處理，忽略了基因讀段分布在樣本間高度相似的特點(diǎn)。因此，提出了一個基于多樣本RNA-Seq數(shù)據(jù)的表達(dá)水平估計方法，稱為MRSeq。其關(guān)鍵是通過建立偏差曲線估計模型獲得基因讀段分布在樣本之間的共享特征，通過偏差權(quán)重將共享特征嵌入到模型中，用來修正讀段數(shù)據(jù)，同時通過增加稀疏約束來表現(xiàn)基因和異構(gòu)體表達(dá)水平之間的稀疏性。進(jìn)而將該模型應(yīng)用到多個真實(shí)數(shù)據(jù)集進(jìn)行評測，與目前主流方法的比較結(jié)果表明：MRSeq不僅能得到準(zhǔn)確的基因和異構(gòu)體表達(dá)水平，同時也獲得了更有意義的生物解釋。

RNA-Seq；多樣本；偏差曲線；稀疏；基因和異構(gòu)體表達(dá)水平

1　引言

近幾年來，下一代高通量測序技術(shù)得到快速發(fā)展，RNA-Seq（RNA sequencing）測序技術(shù)被廣泛應(yīng)用到轉(zhuǎn)錄組學(xué)的研究上[1-2]。與傳統(tǒng)的基因芯片技術(shù)相比，RNA-Seq可在全基因組范圍內(nèi)進(jìn)行測序，具有高通量，高靈敏度，可重復(fù)性好，樣本需求低等特點(diǎn)，具有廣泛的應(yīng)用場景，比如估計基因或異構(gòu)體表達(dá)水平，尋找差異表達(dá)基因或異構(gòu)體，發(fā)現(xiàn)未知的異構(gòu)體等，正快速成為研究轉(zhuǎn)錄組學(xué)的基本實(shí)驗(yàn)手段[3-4]。

當(dāng)從RNA-Seq測序?qū)嶒?yàn)獲得海量讀段數(shù)據(jù)后，計算基因和異構(gòu)體的表達(dá)水平是研究轉(zhuǎn)錄組學(xué)中最為基礎(chǔ)的實(shí)驗(yàn)?zāi)康腫5-6]。rSeq（RNA-Seq analyzer）方法假設(shè)基因的每個外顯子上的讀段數(shù)目服從泊松分布，其參數(shù)是基因所包含的每個異構(gòu)體表達(dá)水平的線性加權(quán)和。rSeq方法解決了因選擇性剪接而造成的讀段多源映射問題，可同時估計基因和異構(gòu)體表達(dá)水平[7]。但rSeq方法是基于基因讀段分布是均勻的假設(shè)，而實(shí)際情況中，讀段分布具有明顯的非均勻特征，該特征通常是由兩類原因造成的：其一是技術(shù)性偏差，在RNA-Seq測序技術(shù)過程中，研究人員經(jīng)常使用不同的方法進(jìn)行RNA提取、大小選擇、片段化、轉(zhuǎn)換為cDNA（complementary DNA）、擴(kuò)增和最后測序，這導(dǎo)致同一個樣本的多次測序結(jié)果存在偏差[8]。其二是生物性偏差，如制備cDNA文庫時，RNA序列對反轉(zhuǎn)錄所采用的隨機(jī)引物具有一定偏好，在PCR擴(kuò)增階段，cDNA片段傾向GC堿基含量高的片段、ploy（A）端和低復(fù)雜度的序列[9-10]。技術(shù)上和生物上的偏差導(dǎo)致讀段分布具有非均勻分布的特征。NURD（non-uniform read distribution）方法改進(jìn)了rSeq模型，通過加入預(yù)先計算的基于基因的全局偏差權(quán)重和基于異構(gòu)體的局部偏差權(quán)重，來模擬讀段分布的非均勻分布特征[11-12]。POM（Poisson mixed-effects model）方法采用基于堿基層面的泊松分布，構(gòu)建一個統(tǒng)計模型來獲得序列中每個堿基的位置偏差和堿基之間的相互作用[13]。GP（generalized Poisson）和WemIQ方法采用廣義泊松分布，通過增加額外參數(shù)來解決因讀段非均勻分布而造成的過散布問題，相比傳統(tǒng)的泊松分布能更好地擬合數(shù)據(jù)[14-15]。Cufflinks方法把偏差特征分為位置偏差和序列偏差，通過分別建模來計算不同偏差的影響，進(jìn)一步量化偏差的構(gòu)成和更準(zhǔn)確地模擬讀段采樣的隨機(jī)特性[16-17]。RSEM（RNASeq by expectation-maximization）和BitSeq（Bayesian inference of transcripts from sequencing Data）方法都采用了產(chǎn)生式概率圖模型，模擬讀段產(chǎn)生過程中的經(jīng)驗(yàn)分布，以及考慮額外的數(shù)據(jù)匹配信息，如匹配質(zhì)量、映射方向等[18-20]。但兩者采用不同的方式來考慮讀段分布的非均勻特征。RSEM利用讀段起始位置的經(jīng)驗(yàn)分布來表示非均勻特征，而BitSeq采用了與Cufflinks同樣的偏差模型[20]。采用不同的偏差計算方式來模擬讀段非均勻分布，上述方法都能提高基因和異構(gòu)體表達(dá)水平的準(zhǔn)確程度。

在RNA-Seq測序?qū)嶒?yàn)中，為了避免實(shí)驗(yàn)中的技術(shù)性誤差，同一個測序樣本會進(jìn)行多次技術(shù)性重復(fù)實(shí)驗(yàn)，從而獲得多樣本RNA-Seq數(shù)據(jù)。但是上述方法通常是逐個處理每個樣本，導(dǎo)致數(shù)據(jù)中的技術(shù)性偏差難以消除，同時不可避免地丟失了樣本之間的相關(guān)信息。如圖1所示（數(shù)據(jù)來自3.1節(jié)中的小鼠數(shù)據(jù)集），老鼠基因APOE在其大腦組織的兩個樣本中，其讀段分布具有高度的相似性，其相關(guān)系數(shù)高達(dá)0.96。說明同一個基因的讀段分布模式在不同樣本中是具有很高相關(guān)性的。在分析轉(zhuǎn)錄組的其他任務(wù)中，如尋找差異基因和異構(gòu)體表達(dá)水平中，BDSeq方法考慮了多樣本RNA-Seq數(shù)據(jù)之間關(guān)聯(lián)獲得更為精確的結(jié)果[21]；尋找新的異構(gòu)體方法中，MITIE（mixed integertranscript identification）方法發(fā)現(xiàn)精度隨著樣本數(shù)目增多而提高[22]。這些方法證明了考慮多樣本數(shù)據(jù)之間的相關(guān)性能幫助提高精度。而在估計表達(dá)水平中，很少有方法結(jié)合多個樣本數(shù)據(jù)來計算基因和異構(gòu)體表達(dá)水平。

在多樣本數(shù)據(jù)中，同一個基因中的某個異構(gòu)體，在不同技術(shù)性重復(fù)樣本中其表達(dá)水平應(yīng)該是相近的。但由于技術(shù)性偏差以及數(shù)據(jù)噪聲，導(dǎo)致有些異構(gòu)體的表達(dá)水平出現(xiàn)偏差。其次同一個基因雖包含多個異構(gòu)體，但基因發(fā)生表達(dá)時通常由少數(shù)異構(gòu)體所體現(xiàn)出來，具有稀疏性的特點(diǎn)[23]。比如人類基因ENSG000000040597包含6個異構(gòu)體，表1是采用Cufflinks方法計算出的6個異構(gòu)體在不同大腦樣本中的表達(dá)水平。從表1中可看出，在不同樣本中基因ENSG000000040597通常有4或5個異構(gòu)體發(fā)生表達(dá)，但起主要作用的是由ENST00000000233和ENST00000489671兩個異構(gòu)體表現(xiàn)出來的，ENST00000415666和ENST00000463733兩個異構(gòu)體是不確定的低表達(dá)，ENST00000459680異構(gòu)體不表達(dá)，而ENST00000467281在不同樣本的表達(dá)水平差異很大且表達(dá)值很低，受到偏差和數(shù)據(jù)噪聲的影響。因此，逐個處理單個樣本的方法不能很好地表現(xiàn)生物體本身所蘊(yùn)含的生物特性。

基于上述問題，本文提出了一個基于多樣本RNA-Seq數(shù)據(jù)的回歸模型來估計基因和異構(gòu)體表達(dá)水平，MRSeq（multi-sample RNA-Seq model）?？紤]基因在多樣本中的讀段分布具有高度相似性，設(shè)計了一個聯(lián)合多樣本數(shù)據(jù)的偏差模型，能更為準(zhǔn)確地表現(xiàn)基因讀段分布在樣本之間的共享特征。同時增加了L2/L1稀疏約束，用來消除技術(shù)性偏差給基因和異構(gòu)體在不同樣本中帶來的差異，且保留基因和異構(gòu)體表達(dá)水平之間的稀疏性質(zhì)，更好地表現(xiàn)其生物特性，獲得更有意義的生物解釋。通過老鼠和人類的多個真實(shí)數(shù)據(jù)集來驗(yàn)證MRSeq的有效性。

2　方法

2.1MRSeq模型

因選擇性剪接在真核生物中是普遍存在的，外顯子片段的選取導(dǎo)致基因通常包含多個異構(gòu)體。假設(shè)基因g包含K個異構(gòu)體，其異構(gòu)體與基因之間的關(guān)系可使用數(shù)學(xué)模型二元矩陣A=(aik)M×K準(zhǔn)確表示，其中M表示基因g包含M個外顯子，其長度分別為l1,l2…,lm。每個異構(gòu)體對應(yīng)矩陣A中的一列，其中(aik)=1表示第k個異構(gòu)體包含第i個外顯子，反之則為0。

Fig.1　Read distribution ofAPOE gene in two mouse brain samples圖1　老鼠基因APOE在兩個大腦組織樣本中的讀段分布

假設(shè)RNA-Seq測序?qū)嶒?yàn)獲得N個樣本，對于基因g，yij表示第i個外顯子在第j個樣本中觀測到的讀段數(shù)目。根據(jù)實(shí)驗(yàn)原理，yij等同于基因g在第j個樣本中所包含異構(gòu)體中第i個外顯子上的讀段之和：

其中，xkj表示第j個樣本中第k個異構(gòu)體上期望的讀段數(shù)目。當(dāng)獲得xkj后，便可計算異構(gòu)體的FPKM（fragments per kilobase of transcript per million mapped reads）值，Wj是第j個樣本中的總讀段數(shù)目[15]。

式（1）模型是假設(shè)讀段在均勻分布的情況下，但在實(shí)際數(shù)據(jù)中，讀段具有明顯的非均勻分布特征。假設(shè)bi表示第i個外顯子的偏差權(quán)重。因?yàn)榛虻淖x段分布模式在不同樣本中是高度相似的，所以假設(shè)bi在不同樣本中是共享的。因此式（1）可改寫成：

對于基因g，其K個異構(gòu)體在N個樣本中表達(dá)水平X可通過回歸模型估計出來，其公式如下：

對于所有異構(gòu)體表達(dá)水平xkj≥0。

同一個基因中相同異構(gòu)體在不同的技術(shù)性重復(fù)樣本中，其表達(dá)水平應(yīng)該是接近的，但在實(shí)際數(shù)據(jù)中，因?yàn)槠詈蛿?shù)據(jù)噪聲的影響，造成表達(dá)水平在不同樣本中具有差異性，特別是低表達(dá)水平的異構(gòu)體和基因。其次，同一個基因雖包含多個異構(gòu)體，但在基因表達(dá)時通常由少數(shù)異構(gòu)體體現(xiàn)出來[23]，基因和異構(gòu)體表達(dá)水平之間是具有稀疏性的。為了融入數(shù)據(jù)中的先驗(yàn)信息，MRSeq方法增加了L2/L1稀疏約束，式（3）可寫成如下：

為了方便計算，式（4）可簡化成矩陣形式：

通過優(yōu)化MRSeq方法的目標(biāo)函數(shù)（5），即可獲得異構(gòu)體在不同樣本中的表達(dá)水平，其基因表達(dá)水平等于對應(yīng)異構(gòu)體表達(dá)水平之和。當(dāng)λ→+∞，X*=0表示異構(gòu)體在不同樣本中都沒有表達(dá)。隨著λ的減小，矩陣X*的某些行不再為0，表示對應(yīng)的異構(gòu)體發(fā)生表達(dá)。本文的所有實(shí)驗(yàn)都選擇λ=1。選擇一個最優(yōu)的λ是至關(guān)重要的，在3.6節(jié)中將討論λ的選擇問題。

Fig.2　Core optimization problem of MRSeq圖2　MRSeq方法的核心優(yōu)化問題

2.2多樣本數(shù)據(jù)的偏差曲線估計

在多樣本數(shù)據(jù)中，基因的讀段分布在樣本之間具有高度相似性。本文提出一種多樣本數(shù)據(jù)的偏差估計模型來表示基因讀段分布在樣本之間的共享特征。選擇只包含單個異構(gòu)體的基因，因?yàn)槎喈悩?gòu)體的基因其讀段分布受到基因結(jié)構(gòu)的影響，不能正確反映基因中讀段的分布趨勢。由于低表達(dá)水平基因的不確定性，本文排除讀段數(shù)目小于50的基因。為了避免讀段分布的局部影響，把基因的長度劃為S個區(qū)間（通常選擇20），統(tǒng)計每個區(qū)間內(nèi)的讀段數(shù)目，采用均值為一的方法來歸一化區(qū)間的讀段數(shù)目。采用多項(xiàng)式回歸模型來擬合歸一化后每個區(qū)間的讀段數(shù)目，估計出光滑的曲線來描述多樣本數(shù)據(jù)的偏差曲線。具體算法步驟如圖3中步驟1所示。

2.3模型實(shí)現(xiàn)

在獲得多樣本讀段匹配數(shù)據(jù)，MRSeq方法的實(shí)現(xiàn)分為兩個步驟：數(shù)據(jù)預(yù)處理和表達(dá)水平計算。MRSeq方法預(yù)處理過程主要包括統(tǒng)計基因外顯子的讀段數(shù)目以及計算基因的偏差權(quán)重。在表達(dá)水平計算步驟上，因?yàn)槟繕?biāo)函數(shù)是帶L2/L1稀疏約束的凸優(yōu)化問題，采用SPAMS（sparse modeling software）優(yōu)化工具箱來求解[24]。SPAMS是一個為解決各種稀疏估計問題，提供多種語言接口和跨多平臺的開源優(yōu)化工具箱。MRSeq方法的詳細(xì)流程如算法1所示。

算法1 MRSeq

為了方便用戶理解和使用MRSeq方法，本文提供了一個系統(tǒng)的RNA-Seq實(shí)驗(yàn)數(shù)據(jù)分析通道，如圖4所示。當(dāng)RNA-Seq測序?qū)嶒?yàn)獲得N個讀段數(shù)據(jù)樣本后，使用最流行的匹配軟件（Bowtie（v.2.2.3））來匹配讀段到已知的參考轉(zhuǎn)路本序列[25]。匹配成功的讀段數(shù)據(jù)作為MRSeq方法的輸入，而MRSeq方法輸出的基因和異構(gòu)體表達(dá)水平可提供給后續(xù)分析使用，比如尋找差異基因或異構(gòu)體表達(dá)，基因網(wǎng)絡(luò)分析等。

Fig.3　Workflow of bias weight calculation圖3　計算基因的偏差權(quán)重的流程圖

Fig.4　Pipeline of RNA-Seq data analysis圖4　RNA-Seq數(shù)據(jù)分析通道的流程圖

3　實(shí)驗(yàn)結(jié)果分析

本文選擇3個主流方法Cufflinks（v.2.2.1）、RSEM（v.1.2.19）和NURD（v.1.1.1）與MRSeq方法在3個數(shù)據(jù)集上進(jìn)行比較，用來驗(yàn)證基因和異構(gòu)體表達(dá)水平。

3.1實(shí)驗(yàn)數(shù)據(jù)集

MRSeq方法可同時估計基因和異構(gòu)體表達(dá)水平，本文選擇3個真實(shí)RNA-Seq數(shù)據(jù)集來評估分法性能。數(shù)據(jù)集都來自Illumina/solexa測序平臺，基因注釋信息都來自UCSCGenome Browser。

小鼠數(shù)據(jù)集包含大腦、肝臟和骨骼肌3個組織，其中每個組織分別包含兩個技術(shù)性重復(fù)實(shí)驗(yàn)樣本。使用RefSeq數(shù)據(jù)庫的基因注釋信息（GRCm38/mm10），總共包含33 608個異構(gòu)體，主要用來驗(yàn)證樣本之間的異構(gòu)體表達(dá)水平[26]。

人類大腦數(shù)據(jù)集來自美國藥品監(jiān)管局（FDA）聯(lián)合全球多所高校研究機(jī)構(gòu)進(jìn)行的“生物芯片質(zhì)量控制（MAQC）”項(xiàng)目。MAQC數(shù)據(jù)集是評估計算基因表達(dá)水平模型的標(biāo)準(zhǔn)數(shù)據(jù)集，被最為廣泛地使用。此數(shù)據(jù)集包含單末端（single-end,SE）和雙末端（pairedend,PE）兩種類型的讀段數(shù)據(jù)。單末端讀段數(shù)據(jù)集包含兩個條件，正常大腦組織（HBR）和病變大腦組織（UHR），每個條件下包括7次測序通道，相當(dāng)于7次技術(shù)性重復(fù)實(shí)驗(yàn)樣本。雙末端讀段數(shù)據(jù)集只有正常大腦組織（HBR），包括3次測序通道。使用Ensembl數(shù)據(jù)庫的基因注釋信息（GRCh37/hg19）。MAQC項(xiàng)目提供了1 000個qRT-PCR（quantitative real-time PCR）驗(yàn)證基因，根據(jù)與Ensembl注釋庫的對應(yīng)匹配，最終獲得833個qRT-PCR驗(yàn)證基因。這些基因的qRT-PCR值被認(rèn)為是基因的真實(shí)表達(dá)水平，可被用來評估模型計算基因表達(dá)水平的準(zhǔn)確程度[27]。

人類乳腺癌數(shù)據(jù)集有兩個條件，乳腺正常細(xì)胞（HME）和乳腺癌細(xì)胞（MCF-7），分別包括4個和7個測序通道。文獻(xiàn)[28]提供了5個基因中8個異構(gòu)體的qRT-PCR驗(yàn)證值，可用來評估模型計算異構(gòu)體表達(dá)水平的準(zhǔn)確程度。根據(jù)文獻(xiàn)選擇UCSC注釋庫的基因注釋信息（GRCh36/hg18）[28]。

3.2多樣本數(shù)據(jù)的偏差曲線

為了驗(yàn)證多樣本偏差曲線估計模型，選擇MAQC在HBR條件下的單末端數(shù)據(jù)集，采用圖3中步驟1的計算流程，把過濾后的基因分成20個區(qū)間進(jìn)行讀段數(shù)目統(tǒng)計和歸一化，選擇多項(xiàng)式回歸模型來擬合7個樣本獲得的數(shù)據(jù)，最終得到如圖5（a）所示的偏差曲線。圖5（a）中的偏差曲線顯示在基因的讀段分布是具有明顯非均勻分布特征的?；虻?′端和5′端是最容易受到偏差影響的[12,17]，正好對應(yīng)著偏差曲線中的兩端。

Fig.5　Bias curves of multi-samples圖5　多樣本偏差曲線

3.3多樣本數(shù)據(jù)之間的表達(dá)水平驗(yàn)證

在多樣本數(shù)據(jù)中，同一個基因中相同異構(gòu)體在不同技術(shù)性重復(fù)樣本中其表達(dá)水平應(yīng)該是相近的，且當(dāng)基因發(fā)生表達(dá)時通常是由其包含的少數(shù)異構(gòu)體所表現(xiàn)出來的，此處主要驗(yàn)證異構(gòu)體在樣本之間的表達(dá)水平。表1給出了Cufflinks方法估計的異構(gòu)體表達(dá)水平，基因ENSG000000040597通常有4或5個異構(gòu)體發(fā)生表達(dá)，且ENST00000467281在不同樣本中差異較大且表達(dá)值很低，不太符合實(shí)際情況。MRSeq方法通過L2/L1稀疏約束來考慮生物體所包含的生物特性。表2是MRSeq方法估計ENSG0000-00040597基因中6個異構(gòu)體的表達(dá)水平。發(fā)現(xiàn)ENSG-000000040597基因表達(dá)主要由ENST00000000233、ENST00000463733和ENST00000489671異構(gòu)體來體現(xiàn)，而ENST00000415666、ENST00000459680和ENST-00000467281異構(gòu)體完全不表達(dá)。這符合基因發(fā)生表達(dá)時通常是由其包含的少數(shù)異構(gòu)體所表現(xiàn)出來的生物特性。

為了進(jìn)一步驗(yàn)證基因中相同異構(gòu)體在不同技術(shù)性重復(fù)樣本中其表達(dá)水平應(yīng)該是相近的生物特性，本文選擇小鼠數(shù)據(jù)集來驗(yàn)證異構(gòu)體樣本之間的表達(dá)水平。因?yàn)镽NA-Seq測序技術(shù)的靈敏度很高，導(dǎo)致相關(guān)系數(shù)極容易受到少數(shù)高表達(dá)異構(gòu)體的影響，為了避免這個問題，對表達(dá)水平對數(shù)化后再進(jìn)行相關(guān)系數(shù)的計算。后續(xù)數(shù)據(jù)集的比較采用同樣的相關(guān)系數(shù)計算方式。在表3中，MRSeq方法在大腦、肝臟和骨骼肌3個條件下都獲得了比其他3個方法更好的結(jié)果，說明MRSeq方法能盡可能保持基因中相同異構(gòu)體在不同樣本中表達(dá)水平應(yīng)該是相近的生物特性。

3.4qRT-PCR基因的表達(dá)水平驗(yàn)證

MAQC數(shù)據(jù)集因提供了約1 000個qRT-PCR驗(yàn)證基因，被最為廣泛地用來比較各個方法的優(yōu)劣性。通過映射到Ensembl注釋庫，最終匹配到833個基因，計算不同方法估計出的基因表達(dá)水平與qRTPCR驗(yàn)證的基因表達(dá)值之間的相關(guān)系數(shù)。由表4的結(jié)果可看出，MRSeq方法在單末端數(shù)據(jù)集上稍微優(yōu)于Cufflinks和RSEM，但明顯優(yōu)于NURD，而在雙末端數(shù)據(jù)集上，MRSeq方法的優(yōu)勢就較為明顯。結(jié)果表明，MRSeq方法相比其他方法在基因表達(dá)水平估計上取得了較為準(zhǔn)確的結(jié)果。

Table 1　Six isoforms expression of ENSG000000040597 gene estimated by Cufflinks in different samples表1　Cufflinks方法估計基因ENSG000000040597中6個異構(gòu)體在不同樣本中的表達(dá)水平

Table 2　Six isoforms expression of ENSG000000040597 gene estimated by MRSeq in different samples表2　MRSeq方法估計基因ENSG000000040597中6個異構(gòu)體在不同樣本中的表達(dá)水平

3.5qRT-PCR異構(gòu)體表達(dá)水平驗(yàn)證

人類乳腺癌數(shù)據(jù)集被用來驗(yàn)證不同方法在異構(gòu)體表達(dá)水平上的準(zhǔn)確程度。文獻(xiàn)[28]對8個異構(gòu)體進(jìn)行了qRT-PCR實(shí)驗(yàn)，其qRT-PCR值被作為標(biāo)準(zhǔn)值。計算不同方法估計出的異構(gòu)體表達(dá)水平與qRT-PCR驗(yàn)證的標(biāo)準(zhǔn)值之間的相關(guān)系數(shù)。從表5中可看出，相比其他方法，MRSeq方法在兩個不同條件下獲得更高的相關(guān)系數(shù)。雖然在MCF-7條件上，4種方法的相關(guān)系數(shù)都非常低，但MRSeq方法仍能獲得比其他方法更好的結(jié)果。結(jié)果表明，MRSeq方法在計算異構(gòu)體表達(dá)水平上取得了較為準(zhǔn)確的結(jié)果。

Table 3　Correlation coefficients of estimated isoform expression within samples by various methods in mouse dataset表3　小鼠數(shù)據(jù)集上各種方法估計的異構(gòu)體表達(dá)水平在樣本間的相關(guān)系數(shù)

Table 4　Correlation coefficient of estimated gene expression by various methods with qRT-PCR results in MAQC dataset表4　MAQC數(shù)據(jù)集上各種方法與qRT-PCR驗(yàn)證基因間的相關(guān)系數(shù)

Table 5　Correlation coefficient of estimated isoform expression by various methods with qRT-PCR results in human breast dataset表5　人類乳腺癌數(shù)據(jù)集上各種方法與qRT-PCR驗(yàn)證異構(gòu)體間的相關(guān)系數(shù)

3.6λ參數(shù)選擇的分析

MRSeq方法通過L2 L1稀疏約束來考慮生物中的生物特性，參數(shù)λ的選擇對表達(dá)水平計算有著很大影響，特別是異構(gòu)體表達(dá)水平的計算。當(dāng)λ→+∞，異構(gòu)體在不同樣本中都沒有表達(dá)。隨著λ的減小，基因中將會有異構(gòu)體表達(dá)出來。選擇MAQC數(shù)據(jù)集在HBR條件下的單末端數(shù)據(jù)集來分析不同λ的選擇對結(jié)果產(chǎn)生的影響。圖6（a）顯示與qRT-PCR驗(yàn)證基因的相關(guān)系數(shù)，隨著λ增大，相關(guān)系數(shù)有著明顯下降，但是在區(qū)間[0,1]內(nèi)，相關(guān)系數(shù)變化不大，表明λ的選擇不應(yīng)過大，較大的λ約束會導(dǎo)致一些真正的低表達(dá)基因或異構(gòu)體被消除。而從圖6（b）中可以看出，基因和異構(gòu)體在樣本之間的相關(guān)系數(shù)隨著λ的增大有顯著的升高，表明λ的增大能消除樣本之間的技術(shù)性偏差而引起的噪聲。為了權(quán)衡λ對表達(dá)水平計算的影響，根據(jù)在多個數(shù)據(jù)集上的結(jié)果，λ=1是個較優(yōu)的選擇，本文所有實(shí)驗(yàn)結(jié)果都采用此設(shè)置。

Fig.6　Effect ofλselection for expression estimation圖6　不同λ選擇對表達(dá)水平估計的影響

4　結(jié)論與展望

本文提出了一個基于多樣本RNA-Seq數(shù)據(jù)的回歸模型來估計基因和異構(gòu)體表達(dá)水平，增加L2 L1稀疏約束來獲得更好的生物解釋。為了準(zhǔn)確描述基因的讀段分布在不同樣本之間的高度相似性，設(shè)計了一個多樣本數(shù)據(jù)的偏差曲線估計模型，通過多項(xiàng)回歸模型獲得偏差曲線，用來表示基因的讀段分布在樣本之間的共享特征。通過偏差曲線可以看出，基因的讀段分布具有非均勻分布性質(zhì)，在分布兩端受到技術(shù)性和生物性偏差的影響，有明顯的上升或下降的趨勢，通過計算每個基因的偏差權(quán)重把讀段分布的非均勻信息融入到表達(dá)水平的估計中。在小鼠數(shù)據(jù)集上，MRSeq方法估計的異構(gòu)體表達(dá)水平在不同樣本之間獲得更高的相關(guān)系數(shù)，表示增加L2 L1稀疏約束可消除異構(gòu)體在不同樣本中受偏差和噪聲而引起的差異，且保留了基因和異構(gòu)體表達(dá)水平之間的稀疏性質(zhì)，從而更好地保留了生物體中所包含的生物特性，獲得了更好的生物解釋性。在人類大腦和乳腺數(shù)據(jù)集上，通過與qRT-PCR基因和異構(gòu)體的驗(yàn)證，MRSeq方法相比其他3種方法能獲得更為準(zhǔn)確的基因和異構(gòu)體表達(dá)水平。在未來研究中，MRSeq方法基于已知的注釋信息，但是生物體的注釋信息遠(yuǎn)遠(yuǎn)沒有達(dá)到完備狀態(tài)，因此通過RNA-Seq測序數(shù)據(jù)來發(fā)現(xiàn)基因中未知異構(gòu)體是很有意義的工作。而MRSeq方法可簡單地擴(kuò)展到此任務(wù)中，利用讀段的跨結(jié)合區(qū)信息來預(yù)測可能存在的異構(gòu)體，增加模型中基因和異構(gòu)體的映射關(guān)系，通過L2 L1稀疏約束來尋找最有可能的未知異構(gòu)體結(jié)構(gòu)，將在后續(xù)的工作來驗(yàn)證此想法。

References:

[1]Marioni J C,Mason C E,Mane S M,et al.RNA-Seq:an assessment of technical reproducibility and comparison with gene expression arrays[J].Genome Research,2008,18(9):1509-1517.

[2]Margueat S,B?hler J.RNA-Seq:from technology to biology[J].Cellular and Molecular Life Sciences,2010,67(4): 569-579.

[3]Marguerat S,Wilhelm B T,B?hler J.Next-generation sequencing:applications beyond genomes[J].Biochemical Society Transactions,2008,36(5):1091-1096.

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

[5]Burgess D J.Gene expression:a global assessment of RNA-seq performance[J].Nature Reviews Genetics,2014,15 (10):645-645.

[6]J?nes J,Hu Fengyuan,Lewin A,et al.A comparative study of RNA-seq analysis strategies[J].Briefings in Bioinformatics, 2015,16(6):1-9.

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

[8]Hansen K D,Brenner S E,Dudoit S.Biases in illuminatranscriptome sequencing caused by random hexamer priming[J]. NucleicAcids Research,2010,28(12):e131.

[9]Risso D,Schwartz K,Sherlock G,et al.GC-content normalization for RNA-Seq data[J].BMC Bioinformatics,2011, 12(1):480.

[10]Huang Yan,Hu Yin,Jones C D,et al.A robust method for transcript quantification with RNA-Seq data[J].Journal of Computational Biology,2013,20(3):167-187.

[11]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(1):220.

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

[13]Hu Ming,Zhu Yu,Taylor J M G,et al.Using Poisson mixed-effects model to quantify transcript-level gene expression in RNA-Seq[J].Bioinformatics,2012,28(1):63-68.

[14]Srivastava S,Chen Liang.A two-parameter generalized Poisson model to improve the analysis of RNA-seq data[J]. NucleicAcids Research,2010,38(17):e170.

[15]Zhang Jing,Kuo C-C J,Chen Liang.WemIQ:an accurate and robust isoform quantification method for RNA-seq data[J]. Bioinformatics,2015,31(6):878-885.

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

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

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

[19]Li Bo,Newey C.RSEM:accurate transcript quantification from RNA-Seq data with or without a reference genome[J]. BMC Bioinformatics,2011,12(1):323.

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

[21]Zhang Li,Liu Xuejun,Chen Songcan.Detecting differential expression from RNA-seq data with expression measurement uncertainty[J].Frontiers of Computer Science,2015,9(4): 652-663.

[22]Behr J,Kahles A,Zhong Yi,et al.MITIE:simultaneous RNA-Seq-based transcript identification and quantification in multiple samples[J].Bioinformatics,2013,29(20):2529-2538.

[23]Zheng Xia,Wen Jianguo,Cheng Chung-Che,et al.NSMAP:a method for spliced isoforms identification and quantification from RNA-Seq[J].BMC Bioinformatics,2011,12(1):162.

[24]Jenatton R,Mairal J,Obozinski G,et al.Proximal methods for sparse hierarchical dictionary learning[C]//Proceedings of the 27th International Conference on Machine Learning, Haifa,Israel,Jun 21-24,2010.

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

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

[27]Canales R D,Luo Yuling,Willey J C,et al.Evaluation of DNA microarray results with quantitative gene expression platforms[J].Nature Biotechnology,2006,24(9):1115-1122.

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

張禮（1985—），男，2010年于南京航空航天大學(xué)計算機(jī)應(yīng)用專業(yè)獲得碩士學(xué)位，現(xiàn)為南京航空航天大學(xué)博士研究生，主要研究領(lǐng)域?yàn)樯镄畔W(xué)，機(jī)器學(xué)習(xí)等。

劉學(xué)軍（1976—），女，2006年于英國曼切斯特大學(xué)獲得博士學(xué)位，現(xiàn)為南京航空航天大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院教授，主要研究領(lǐng)域?yàn)樯镄畔W(xué)，機(jī)器學(xué)習(xí)等。

陳松燦（1961—），男，1997年于南京航空航天大學(xué)獲得博士學(xué)位，現(xiàn)為南京航空航天大學(xué)教授、博士生導(dǎo)師，主要研究領(lǐng)域?yàn)闄C(jī)器學(xué)習(xí)，模式識別等。發(fā)表學(xué)術(shù)論文100余篇，主持國家自然科學(xué)基金、江蘇省自然科學(xué)基金等多個項(xiàng)目。

Novel Method to Estimate Expression Level Based on Multi-Sample RNA-Seq Data*

ZHANG Li+,LIU Xuejun,CHEN Songcan
College of Computer Science&Technology,Nanjing University ofAeronautics&Astronautics,Nanjing 210016,China
+Corresponding author:E-mail:leo.zhang@nuaa.eud.cn

With the rapid development of the next-generation high-throughput sequencing technology,RNA-Seq has become the standard and important technique for transcriptome analysis.For multi-sample RNA-Seq data,the existing expression estimation methods usually deal with each single RNA-Seq sample,and ignore the read distributions with high consistency between multiple samples.This paper proposes a novel method,MRSeq,to estimate expression using multi-sample RNA-Seq data.MRSeq introduces a bias curve estimation model to capture the common features of read distributions shared among multiple samples.The common features are embedded into the model by deviation weight to correct read distributions.Meanwhile,by adding a sparse constraint,the method considers the sparsity between gene and the corresponding isoform expression.Three real datasets are used to validate the proposed method on gene and isoform expression estimation.Compared with the popular methods,MRSeq obtains more accurate gene and isoform expression estimation,and more meaningful biological explanation.

RNA-Seq;multi-sample;bias curve;sparse-specific;gene and isoform expression

2015-04,Accepted 2015-06.

ZHANG Li was born in 1985.He the M.S.degree in computer applications from Nanjing University of Aeronautics andAstronautics in 2010.Now he is a Ph.D.candidate at Nanjing University ofAeronautics andAstronautics.His research interests include bioinformatics and machine learning,etc.

LIU Xuejun was born in 1976.She the Ph.D.degree in computer science from University of Manchester, UK in 2006.Now she is a professor at College of Computer Science and Technology,Nanjing University of Aeronautics andAstronautics.Her research interests include bioinformatics and machine learning,etc.

CHEN Songcan was born in 1961.He the Ph.D.degree in communication and information systems from Nanjing University of Aeronautics and Astronautics in 1997.Now he is a Professor and Ph.D.supervisor at College of Computer Science and Technology,Nanjing University of Aeronautics and Astronautics.His research interests include machine learning and pattern recognition,etc.

10.3778/j.issn.1673-9418.1505045

*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ù)費(fèi)專項(xiàng)資金).

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

A

TP391