基于多樣本RNA-Seq數據的表達水平估計方法*

張　禮，劉學軍，陳松燦

南京航空航天大學 計算機科學與技術學院,南京 210016

基于多樣本RNA-Seq數據的表達水平估計方法*

張禮+，劉學軍，陳松燦

南京航空航天大學 計算機科學與技術學院,南京 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測序的快速發展，RNA-Seq測序已成為轉錄組學分析的標準技術。在處理多樣本RNA-Seq數據時，現有表達水平估計方法通常基于單個樣本逐個處理，忽略了基因讀段分布在樣本間高度相似的特點。因此，提出了一個基于多樣本RNA-Seq數據的表達水平估計方法，稱為MRSeq。其關鍵是通過建立偏差曲線估計模型獲得基因讀段分布在樣本之間的共享特征，通過偏差權重將共享特征嵌入到模型中，用來修正讀段數據，同時通過增加稀疏約束來表現基因和異構體表達水平之間的稀疏性。進而將該模型應用到多個真實數據集進行評測，與目前主流方法的比較結果表明：MRSeq不僅能得到準確的基因和異構體表達水平，同時也獲得了更有意義的生物解釋。

RNA-Seq；多樣本；偏差曲線；稀疏；基因和異構體表達水平

1　引言

近幾年來，下一代高通量測序技術得到快速發展，RNA-Seq（RNA sequencing）測序技術被廣泛應用到轉錄組學的研究上[1-2]。與傳統的基因芯片技術相比，RNA-Seq可在全基因組范圍內進行測序，具有高通量，高靈敏度，可重復性好，樣本需求低等特點，具有廣泛的應用場景，比如估計基因或異構體表達水平，尋找差異表達基因或異構體，發現未知的異構體等，正快速成為研究轉錄組學的基本實驗手段[3-4]。

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

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

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

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

2　方法

2.1MRSeq模型

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

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

假設RNA-Seq測序實驗獲得N個樣本，對于基因g，yij表示第i個外顯子在第j個樣本中觀測到的讀段數目。根據實驗原理，yij等同于基因g在第j個樣本中所包含異構體中第i個外顯子上的讀段之和：

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

式（1）模型是假設讀段在均勻分布的情況下，但在實際數據中，讀段具有明顯的非均勻分布特征。假設bi表示第i個外顯子的偏差權重。因為基因的讀段分布模式在不同樣本中是高度相似的，所以假設bi在不同樣本中是共享的。因此式（1）可改寫成：

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

對于所有異構體表達水平xkj≥0。

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

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

通過優化MRSeq方法的目標函數（5），即可獲得異構體在不同樣本中的表達水平，其基因表達水平等于對應異構體表達水平之和。當λ→+∞，X*=0表示異構體在不同樣本中都沒有表達。隨著λ的減小，矩陣X*的某些行不再為0，表示對應的異構體發生表達。本文的所有實驗都選擇λ=1。選擇一個最優的λ是至關重要的，在3.6節中將討論λ的選擇問題。

Fig.2　Core optimization problem of MRSeq圖2　MRSeq方法的核心優化問題

2.2多樣本數據的偏差曲線估計

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

2.3模型實現

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

算法1 MRSeq

為了方便用戶理解和使用MRSeq方法，本文提供了一個系統的RNA-Seq實驗數據分析通道，如圖4所示。當RNA-Seq測序實驗獲得N個讀段數據樣本后，使用最流行的匹配軟件（Bowtie（v.2.2.3））來匹配讀段到已知的參考轉路本序列[25]。匹配成功的讀段數據作為MRSeq方法的輸入，而MRSeq方法輸出的基因和異構體表達水平可提供給后續分析使用，比如尋找差異基因或異構體表達，基因網絡分析等。

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

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

3　實驗結果分析

本文選擇3個主流方法Cufflinks（v.2.2.1）、RSEM（v.1.2.19）和NURD（v.1.1.1）與MRSeq方法在3個數據集上進行比較，用來驗證基因和異構體表達水平。

3.1實驗數據集

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

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

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

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

3.2多樣本數據的偏差曲線

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

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

3.3多樣本數據之間的表達水平驗證

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

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

3.4qRT-PCR基因的表達水平驗證

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

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

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

3.5qRT-PCR異構體表達水平驗證

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

Table 3　Correlation coefficients of estimated isoform expression within samples by various methods in mouse dataset表3　小鼠數據集上各種方法估計的異構體表達水平在樣本間的相關系數

Table 4　Correlation coefficient of estimated gene expression by various methods with qRT-PCR results in MAQC dataset表4　MAQC數據集上各種方法與qRT-PCR驗證基因間的相關系數

Table 5　Correlation coefficient of estimated isoform expression by various methods with qRT-PCR results in human breast dataset表5　人類乳腺癌數據集上各種方法與qRT-PCR驗證異構體間的相關系數

3.6λ參數選擇的分析

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

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

4　結論與展望

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

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年于南京航空航天大學計算機應用專業獲得碩士學位，現為南京航空航天大學博士研究生，主要研究領域為生物信息學，機器學習等。

劉學軍（1976—），女，2006年于英國曼切斯特大學獲得博士學位，現為南京航空航天大學計算機科學與技術學院教授，主要研究領域為生物信息學，機器學習等。

陳松燦（1961—），男，1997年于南京航空航天大學獲得博士學位，現為南京航空航天大學教授、博士生導師，主要研究領域為機器學習，模式識別等。發表學術論文100余篇，主持國家自然科學基金、江蘇省自然科學基金等多個項目。

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(國家自然科學基金);the Qinglan Project of Jiangsu Province(江蘇省青藍工程);the Fundamental Research Funds for the Central Universities of China under Grant No.CXZZ11_0217 (中央高校基本科研業務費專項資金).

CNKI網絡優先出版：2015-06-18,http://www.cnki.net/kcms/detail/11.5602.TP.20150618.1646.001.html

A

TP391