多條件多樣本RNA-Seq 數據的剪切異構體表達水平估計

張禮，馬越，吳東洋

（1.南京林業大學 信息科學技術學院，江蘇 南京 210016;2.江蘇健康衛生職業學院 中西醫結合學院，江蘇 南京 210018）

選擇性剪切事件是導致生物體多樣性的重要原因之一。為了進一步揭示選擇性剪切的內在機制，迫切需要計算剪切異構體的表達水平。與傳統的基因芯片技術相比，高通量RNA 測序(RNA sequencing，RNA-Seq)技術具有高通量、高靈敏度、可重復性好等優勢，已成為轉錄組學分析的一個標準技術手段[1-5]。

RNA-Seq 測序實驗獲得海量讀段，將讀段與參考注釋序列進行匹配，之后便可估計剪切異構體的表達水平。但是在估計剪切異構體表達水平的過程中，面臨著兩個最大挑戰，即讀段的多源映射和數據偏差[6-7]。研究者提出了大量剪切異構體表達水平估計方法來解決上述的問題。rSeq方法把讀段映射到外顯子的過程當作一個泊松隨機過程，其泊松分布的參數對應著基因所包含剪切異構體表達水平的線性加權[8]。但是rSeq 方法假設基因上讀段分布是均勻的，這與真實數據分布特點不一致。在真實數據中，讀段分布呈現明顯的非均勻特征。讀段的非均勻分布主要是由測序數據中的各種偏差造成的，比如GC 堿基序列偏差，5 端和3 端的位置偏差以及實驗技術性偏差等。針對偏差所導致問題，NURD 方法考慮了全局和局部位置偏差所帶來的影響[9]。POME 方法考慮了序列中堿基之間的關聯性[10]。為了考慮更復雜的偏差，大量概率生成式模型被提出，其直接模擬讀段的隨機采樣過程。Cufflinks 方法設計了不同的模型來消除序列偏差和位置偏差的影響，從而更加準確地描述讀段隨機采樣過程[11]。BitSeq 和PBSeq 方法采用了與Cufflinks 同樣的偏差估計模型[12-13]。RSEM 方法考慮了讀段匹配的不確定性因素，并且使用了讀段起始位置的經驗分布來表示讀段在基因上的非均勻分布特征，但是其未考慮序列偏差這個重要因素[14]。上述方法采用不同的偏差估計模型來模擬讀段的非均勻分布特征，都能提高剪切異構體表達水平的估計準確程度。

由于數據噪聲和偏差的影響，異構體表達水平的準確性仍然有較大提高的空間[15-16]。常規的RNA-Seq 測序實驗通常會設置不同的實驗條件，比如：同一個細胞組織下參照組和對照組，不同時間點下胚胎發育狀況等。此外為了避免實驗中的技術性誤差，同一個實驗條件下會進行多次重復性技術性實驗。這使得一次測序實驗獲得的RNA-Seq 數據集是一個多條件多樣本的數據集。但是上述方法都是假設RNA-Seq 數據集中各個樣本之間是相互獨立，因此都是單獨逐個處理每個數據樣本。這導致樣本之間的相關性沒有得到充分利用。因此有少量工作開始探索聯合多樣本RNA-Seq 數據進行異構體表達水平估計[17-18]。Sequgio 方法能從多樣本數據中自動獲取位置偏差和局部序列影響，再通過對聯合統計模型添加一個光滑的正則化項，來控制讀段在多樣本的一致性[19]。MSIQ 方法考慮多樣本之間的異質性所導致的結果不穩定性，首先將同質性相近的樣本歸為同一組，然后在貝葉斯框架模型下，給同一組之內的樣本賦予較高的權重，從而獲得更加魯棒的異構體表達水平[20]。XAEM 方法采用雙線性模型同時估計異構體表達水平和數據偏差，該模型能夠自動對潛在的未知偏差進行經驗校正[21]。但是上述方法所處理的多樣本數據，僅僅是針對單條件下的多樣本，比如同一個組織細胞的對照組或者同一個時間點狀態。當處理多條件多樣本數據時，這些方法都是假設各個條件之間不相關，把多條件多樣本數據拆分為多個單條件多樣本數據集來進行異構體表達水平計算。但是基因讀段分布在不同條件下同樣具有高度相似性[22]。為了充分利用數據信息，PGSeq 方法采用泊松分布和伽瑪分布的混合模型聯合估計基因和異構體表達水平，其伽瑪分布用來模擬基因讀段分布在多條件多樣本下的偏差信息[23]。但PGSeq 方法未考慮到基因和異構體表達水平之間的稀疏特性，易受到數據噪聲的影響。

基于上述問題，本文提出了一個多條件多樣本RNA-Seq 測序數據異構體表達水平估計方法，MCMS-Seq(multi-condition multi-sample RNASeq)。該模型考慮了基因讀段分布在不同條件下的樣本具有高度相似性，設計一個聯合多條件多樣本數據的偏差估計模型，同時考慮了基因讀段分布受全局偏差和局部偏差的影響。此外，MCMS-Seq方法增加了L2/L1組稀疏約束和L1稀疏約束兩個正則化項，用來體現基因和剪切異構體之間存在稀疏特性，以及消除技術性誤差和數據噪聲的影響。最后，通過3 個多條件多樣本RNA-Seq 數據集來評估MCMS-Seq 方法的性能。

1 MCMS-Seq 方法

1.1 MCMS-Seq 模型表示

由于選擇性剪切事件在真核生物中普遍存在，這給計算剪切異構體表達水平帶來了一個最大問題，即如何定量確定匹配到共享外顯子上的讀段來自哪個剪接異構體。圖1 中顯示的基因包含4 個外顯子(Exon)和3 個剪切異構體。其中一個外顯子可以同時被多個剪切異構體共享，比如外顯子1 被3 個剪切異構體共享，但是剪切異構體2 僅共享了外顯子1 的部分序列。針對這類部分共享情況，可將外顯子1 分割為2 個不重疊的外顯子片段。因此該基因的4 個外顯子被分割成7 個完全不重疊的外顯子片段。映射矩陣A表示圖1 中剪切異構體與外顯子片段的關系，其中矩陣元素a12=1 表示異構體1 包含外顯子片段2。當測序讀段匹配到基因上，外顯子片段上的讀段數目即可被統計出來。假設某個數據集有2 個條件每個條件包括2 個樣本，總計4 個樣本，那么圖1 中基因在不同樣本中讀段數據可用數據矩陣D表示。每一行表示該基因在一個樣本中外顯子片段的讀段數目。

圖1 剪切異構體中外顯子片段劃分示例Fig.1 Example of exon segmentation in an isoform

假設測序實驗獲得RNA-Seq 數據包含C個條件，每個條件包含N個樣本。對于基因g，該基因包含K個剪切異構體和M個外顯子片段，其與外顯子片段的映射關系由映射矩陣AM×K表示。yci j表示基因第j個外顯子片段在第c個條件中第i個樣本上的讀段數目。根據實驗原理，基因外顯子片段的讀段數目等于共享該外顯子片段的剪切異構體上讀段之和，其數學模型為

式中：xcik表示基因第k個剪切異構體在第c個條件中第i個樣本；ajk表示剪切異構體與外顯子片段之間的映射關系；wci表示第c個條件中第i個樣本的讀段總數；lj是第j個外顯子片段的長度。

式(1)模型是基于基因讀段是均勻分布假設的前提，但是實際數據中，基因讀段分布呈現明顯的非均勻特征。由于基因讀段分布模式在不同條件不同樣本下具有高度相似性，因此假設bj表示第j個外顯子片段的偏差權重，其值在樣本之間是共享的。現將偏差bj融入到式(1)中，得到如下模型：

對于多條件多樣本的RNA-Seq 數據集，基因g所包含的K個剪切異構體的表達水平X可以通過回歸模型計算，其公式如下：

所有剪切異構體在不同樣本中的表達水平都要求xcjk≥0。為了便于理解和計算，式(3)可以簡化為矩陣形式：

式中D表示歸一化后的數據矩陣。

一個基因雖然包含多個剪切異構體，但是在不同條件下，少數剪切異構體的表達水平決定了該基因的表達。因此基因和剪切異構體表達水平之間具有稀疏特性。通過對剪切異構體表達水平X增加L1范數來保留稀疏特性，式(4)可改寫為

雖然模型增加了L1范數的稀疏約束，但仍然會發現出現大量低表達的剪切異構體，而這部分剪切異構體不全是真實的低表達。當一個剪切異構體在同一個條件下的所有重復樣本都是低表達水平，那么可認為此剪切異構體是真實的低表達。而對于零散出現的低表達剪切異構體，則受到數據噪聲和偏差的影響，不是真實的低表達。為了消除虛假的低表達剪切異構體的影響，在式(5)的基礎上增加了L2/L1組稀疏約束得到了MCMSSeq 方法的最終形式：

式中 λ1和 λ2分別是L2/L1和L1約束的系數。通過兩個稀疏約束項，MCMS-Seq 方法不僅考慮了基因和剪切異構體表達水平之間的稀疏性質，同時也可以消除數據噪聲和偏差對低表達剪切異構體的影響。圖2 顯示了MCMS-Seq 方法的優化問題。

圖2 MCMS-Seq 方法的優化問題Fig.2 Optimization problem of MCMS-Seq method

1.2 多條件多樣本偏差估計模型

在多條件多樣本數據中，圖3 顯示了基因的讀段分布無論在不同條件下，還是在同一個條件的重復樣本中，其分布模式具有高度相似性。MCMS-Seq方法提出了一個基于多條件多樣本的讀段非均勻偏差估計模型。該偏差估計模型由兩部分構成：全局偏差 βglobal和局部偏差 βlocal。全局偏差 βglobal的讀段非均勻分布模式是從數據集中所有表達基因中獲得。由于讀段多源映射會影響基因讀段分布，全局偏差估計僅僅選擇只包含單個剪切異構體的基因。此外，由于低表達水平基因的不確定性，讀段計數小于50 的基因被排除。將篩選后的基因均分為20 個等長度的區間，統計并歸一化每個區間內讀段數目。最后采用多項式回歸來擬合基因每個區間上的讀段數目，得到的擬合曲線表示基因讀段分布的全局偏差特征。而局部偏差βlocal僅僅統計基因每個外顯子片段在多條件多樣本數據上的讀段數目，再進行均一化處理，其反映了單個基因自身的讀段分布特征。

圖3 小鼠數據集中基因Utrn 讀段分布Fig.3 Read distributions of gene Utrn in the mouse dataset

一旦獲得數據集的全局偏差曲線和單個基因的局部偏差特性，便可以計算出基因上每個外顯子片段的偏差值：

式中：α 是權重參數，用來權衡全局偏差和局部偏差的影響。本文選擇α=0.5，表示全局偏差和局部偏差對基因讀段分布具有相同的影響[9]，不僅僅能反映讀段非均勻分布在多條件多樣本之間具有高度相似的特征，同時還可以體現出每個基因獨有讀段分布特點。

1.3 MCMS-Seq 模型實現

MCMS-Seq 方法的實現可以分為3 個部分：讀段數據預處理、基因偏差估計和表達水平估計。

1) 讀段數據預處理，是從匹配成功的讀段數據中統計基因每個外顯子片段的讀段計數，以及從注釋文件中獲得外顯子片段和剪切異構體之間的映射關系矩陣。

2) 基因偏差估計，是計算數據集的全局偏差和基因的局部偏差，從而獲得基因每個外顯子片段的基因偏差值。

3) 剪切異構體表達水平估計，由于模型是針對多條件多樣本數據集，同時模型包含L2/L1和L1約束，MCMS-Seq 方法采用SPAMS 優化工具箱來求解[24-25]。

MCMS-Seq 方法的詳細流程如算法1 所示，采用Python 和MATLAB 混合編程實現。

算法1MCMS-Seq 方法

輸入多條件多樣本數據，注釋文件；

輸出每個基因的剪切異構體表達水平。

1)數據預處理：統計外顯子片段讀段數目矩陣D，構建映射關系矩陣A。

2)基因偏差估計：計算外顯子片段偏差值。

3)表達水平估計：計算所有基因的X?。

1.4 多條件多樣本數據分析通道

為了方便用戶使用MCMS-Seq 方法，本文提供了一個多條件多樣本RNA-Seq 測序數據分析通道，如圖4 所示。當獲得RNA-Seq 測序數據樣本后，使用經典讀段匹配軟件Bowtie[26]，將每個數據樣本的讀段匹配到參考轉錄組參考序列上。每個樣本匹配結果作為輸入數據一并輸入到MCMSSeq 分析通道中，從而可獲得剪切異構體在不同樣本中的表達水平。一旦獲得剪切異構體的表達水平，可提供給高層次的后續分析使用。

圖4 多條件多樣本RNA-Seq 數據分析通道Fig.4 Analysis pipeline of multi-condition multi-sample RNA-Seq data

2 實驗結果與分析

本文選擇了經典方法Cufflinks(v.2.2.1)和PGSeq(v.1.0)，以及最新方法XAEM(v.0.1.1),分別在 3個數據集上與MCMS-Seq 方法進行比較，用來驗證剪切異構體表達水平的性能。針對多條件多樣本數據集，Cufflinks 是每個樣本單獨處理，而PGSeq、XAEM 和MCMS-Seq 都是多個樣本聯合處理。

2.1 數據集

3 個多條件多樣本的RNA-Seq 數據集被用來驗證MCMS-Seq 方法估計剪切異構體表達水平的準確性。3 個數據集分別是小鼠數據集、人類大腦的SEQC 和MAQC-II 數據集，它們都來自Illumina/solexa 測序平臺。

小鼠數據集包含3 個條件，分別是肝臟、大腦和骨骼肌3 個組織，其中每個組織分別包含了 2個重復實驗樣本。使用 RefSeq 數據庫的基因注釋信息(GRCm38/mm10)，總共包含 33608 個剪切異構體，主要用來驗證同條件下重復樣本之間剪切異構體表達水平的可重復性[27]。

MAQC(micorarray quality control)來自美國藥品監管局的生物芯片質量控制項目。該項目分為三期實施，即MAQC-I、MAQC-II 和MAQC-III，其產生的數據集被廣泛應用于評估不同測序平臺下不同方法的性能。本文主要利用了MAQC-II 和MAQC-III 兩期項目提供的數據。MAQC-III 也被稱為SEQC(sequencing quality control)。SEQC 包括兩個實驗條件UHRR(universal human reference rna)和HBRR(human brain reference RNA)，每個條件分別有8 個重復實驗樣本。SEQC 數據集提供了兩萬多個經qRT-PCR 實驗驗證的剪切異構體。與Ensembl 注釋信息(GRCh37/hg19)相匹配后，最終得到16603 個剪切異構體。這些剪切異構體的qRT-PCR 值被當作真實表達水平值，可用來評估模型計算剪切異構體表達水平的準確性[28]。

基因表達水平是由其包含的剪切異構體所構成，因此基因表達水平可用來進一步驗證剪切異構體表達水平的準確性。MAQC-II 數據集同樣包含UHRR 和HBRR 兩個實驗條件，每個條件下包含7 個重復性實驗。該數據提供了1000 個經qRT-PCR 實驗驗證的基因。根據與Ensembl 注釋信息(GRCh37/hg19)相匹配，最終獲得838 個基因。這些基因的qRT-PCR 值被當作真實基因表達水平值，用來間接評估模型計算剪切異構體表達水平的準確程度[29]。

2.2 多條件多樣本偏差估計模型驗證

MCMS-Seq 方法提出了一個基于多條件多樣本偏差估計模型，同時考慮了讀段分布受到全局偏差和局部偏差的影響，用來獲取讀段分布在樣本之間的高度相似性特征。SEQC 數據集被用來驗證偏差估計模型的有效性。圖5 顯示使用該模型對SEQC 數據集的偏差估計流程。從圖5(a)中可以看出，在SEQC 數據集中，基因的讀段分布呈現明顯的非均勻分布特征，特別是在基因的兩端。這個現象符合基因的 3′端和 5′端最容易受到RNA-Seq 測序技術影響的事實。選擇基因Cdca4來展示估計全局偏差和局部偏差的過程。基因Cdca4 包含3 個剪切異構體和5 個外顯子片段，其結構如圖5(b) 所示。圖5(c)是通過多項式回歸擬合圖5(a)讀段分布所得到SEQC 數據集全局偏差曲線。曲線上黑點表示基因Cdca4 外顯子片段長度的比率。通過長度比率在曲線上的取值，可得到Cdca4 基因中每個外顯子片段的全局偏差值。統計并歸一化基因Cdca4 的外顯子片段在所有樣本中的讀段數目，即可獲得該基因的局部偏差，如圖5(d) 所示。從圖5 中可以看出，該基因在3′端和 5′端受到的局部偏差影響要略小于全局偏差。為了進一步驗證基因的局部偏差，從SEQC 數據集中隨機選擇4 個基因：Plagl1、Eif4a、Sv2b 和Whrn，其分別包含5、6、7、8 個剪切異構體。從圖6 中可以看出，不同基因的局部偏差整體上都呈現明顯非均勻分布特征，但是單個基因之間存在一定差異，比如基因Whrn 中間外顯子的偏差值表現出由高到低的趨勢。因此，MCMSSeq 方法提出的多條件多樣本偏差估計模型，不僅能反映在多條件多樣本數據中讀段非均勻分布具有高度相似性的特征，同時還可以體現出單個基因獨有讀段分布特點。

圖5 MCMS-Seq 方法的偏差估計流程Fig.5 Bias estimation process of the MCMS-Seq method

圖6 基因的局部偏差分布Fig.6 Local bias distribution of genes

2.3 多條件多樣本下剪切異構體表達水平的驗證

MCMS-Seq 方法處理多條件多樣本數據集時是聯合所有樣本同時處理，通過增加稀疏約束，不僅可以消除數據噪聲的影響，同時也能體現基因和剪切異構體之間存在的稀疏特性。選擇小鼠數據集的基因Nph2 來驗證，該基因包含3 個剪切異構體。

在小鼠數據集中，同一個剪切異構體在同一個條件下的多個重復樣本中，其表達水平應該是相近的。若一個剪切異構體在重復樣本中零散地出現低表達，則此剪切異構體的表達水平受到數據噪聲的影響。傳統方法Cufflinks 都是每個樣本依次單獨處理，其表達水平值如表1 所示。NM_001 364736 表達水平在Muscle 條件兩個重復樣本中就可能受到數據噪聲的影響，NM_157294 在Liver 條件下也存在同樣的情況。表2 中XAEM 方法獲得的NM_001364736 和NM_157294 表達水平都是極低值，極大可能是受到數據噪聲的干擾。MCMS-Seq 方法聯合處理多條件多樣本數據集。從表3 中可以看出，NM_001364736 在3 個組織條件下都未表達，NM_157294 在大腦和骨骼肌組織條件下具有真實的低表達，而在肝臟組織條件下未表達，能有效消除數據噪聲的影響。

表1 Cufflinks 估計基因Nhp2 中3 個剪切異構體表達水平Table 1 Expression level of three isoforms in Nph2 gene estimated using cufflinks

表2 XAEM 估計基因Nhp2 中3 個剪切異構體表達水平Table 2 Expression level of three isoforms in Nph2 gene estimated using XAEM

表3 MCMS-Seq 估計基因Nhp2 中3 個剪切異構體表達水平Table 3 Expression level of three isoforms in Nph2 gene estimated using MCMS-Seq

此外，基因外在表現通常是由其包含的少數剪切異構體決定的，因此基因和剪切異構體之間存在稀疏特性。在表4 中，PGSeq 方法得到的3 個剪切異構體表達水平都存在較高的表達值，無法體現稀疏特性。而Cufflinks 和XAEM 受數據噪聲影響，同樣很難體現出該數據特性。MCMSSeq 方法增加了L2/L1組稀疏和L1稀疏約束來考慮上述生物特性。NM_026631 在所有組織條件下中都有較高的表達水平，說明基因Nph2 的表達主要由NM_02663 所決定。NM_001364736 在3 個組織條件下都未表達，特別在肝臟組織條件下NM_157294 和NM_001364736 同時未表達，這表明基因Nph2 在肝臟中只有剪切異構體NM_026631 參與基因表達。因此，MCMS-Seq 方法能體現基因表達是由少數剪切異構體所決定的生物特性，提供了更好的生物可解釋性。

表4 PGSeq 估計基因Nhp2 中3 個剪切異構體表達水平Table 4 Expression level of three isoforms in Nph2 gene estimated using PGSeq

2.4 剪切異構體的可重復性驗證

在多條件多樣本測序實驗中，同一個條件下設計多重復性樣本是為了避免技術性誤差所帶來的影響。這使得同一個剪切異構體在同一個條件下的重復樣本之間的表達水平是相近的。小鼠數據集被用來驗證剪切異構體表達水平在樣本之間的可重復性。采用Person 相關系數來評估可重復性，其值越高說明能更加有效地消除技術性誤差所造成的偏差。由于RNA-Seq 測序技術得到表達水平其幅度跨度很大，Person 相關系數易受到少數高表達的剪切異構體影響。因此在計算相關系數之前，對所有剪切異構體表達水平進行對數轉換，從而避免上述問題。表5 中顯示不同方法在小鼠數據集上不同條件下的相關系數值。從表中可以看出，MCMS-Seq 方法在肝臟、大腦和骨骼肌3 個條件下都獲得了比其他3 個方法更好的結果。盡管MCMS-Seq 方法是面向處理多條件多樣本數據集，但仍然可以保證剪切異構體在同一個條件中下樣本之間具有高度的可重復性。這也符合RNA-Seq 測序實驗中設計重復實驗的目的。

表5 在小鼠數據集上不同方法估計的剪切異構體表達水平在樣本之間的相關系數Table 5 Correlation coefficients between isoform expression levels estimated using various methods in the mouse dataset

2.5 PCR 剪切異構體的表達水平驗證

SEQC 數據集被用來驗證不同方法估計剪切異構體表達水平的準確性。該數據集提供了16 603個經過qRT-PCR 驗證的剪切異構體，這些剪切異構體被當作基準數據。計算不同方法得到剪切異構體表達水平與qRT-PCR 值之間的相關系數。從表6 中結果可以看出，MCMS-Seq 方法在UHRR條件上稍微優于PGSeq 方法，而在HBRR 條件上獲得較為明顯的提升。盡管XAEM 方法是多樣本處理，但獲得最差的性能，其可能是該方法對數據偏差考慮得不夠。整體上說，MCMS-Seq 方法估計的剪切異構體表達水平能取得較為準確的結果。

表6 在 SEQC 數據集上不同方法與 qRT-PCR 驗證剪切異構體之間的相關系數Table 6 Correlation coefficients between qRT-PCR values and isoform expression levels estimated using various methods in SEQC dataset

2.6 PCR 基因的表達水平驗證

現實中包含qRT-PCR 驗證的剪切異構體數據集很少，而基因的表達水平是由其所包含的剪切異構體所決定的，因此可以通過驗證qRTPCR 驗證基因的表達水平來間接驗證剪切異構體表達水平的準確性。MAQC-II 數據集被廣泛地應用于評估不同方法估計基因表達水平的性能。MAQC-II 數據集提供了838 個qRT-PCR 驗證的基因，這些基因總共包含了6927 個剪切異構體。Cufflinks 和PGSeq 方法提供了基因的表達水平，XAEM 和MCMS-Seq 方法的基因表達水平由所對應的剪切異構體表達水平求和得到。表7 顯示了不同方法得到的基因表達水平與qRT-PCR 值之間的相關系數。從表7 中可以看出，相比其他方法，MCMS-Seq 方法得到了更好的準確性。

表7 在 MAQC-II 數據集上不同方法與 qRT-PCR 驗證基因之間的相關系數Table 7 Correlation coefficients between qRT-PCR values and isoform expression levels estimated using various methods in the MAQC-II dataset

2.7 稀疏參數的選擇

MCMS-Seq 方法包含了L2/L1組稀疏約束和L1稀疏約束兩個正則化項，不僅用來考慮基因和剪切異構體之間的稀疏特性，同時用來消除虛假低表達剪切異構體帶來的影響。在式(6)中，參數 λ1和 λ2分別對應著L2/L1組稀疏約束和L1稀疏約束，其值的選擇能影響到剪切異構體表達水平的準確性。當 λ1或 λ2→+∞ 時，都會導致剪切異構體出現不表達情況，區別在于，λ1→+∞ 會導致同一個剪切異構體在不同條件下都沒有表達。當λ1→0 時，剪切異構體的表達水平容易受到數據噪聲的影響，產生虛假的低表達。而 λ2減小時，基因與剪切異構體之間的稀疏特性將減弱。選擇SEQC 數據集中HBRR 條件來分析參數選擇對剪切異構體表達水平準確性的影響。假設參數λ1和 λ2分別選擇0.1、1、10 和100 這4 個值，圖7顯示了在取不同參數值時，MCMS-Seq 方法估計的剪切異構體表達水平與qRT-PCR 驗證的剪切異構體之間的相關系數。從圖7 可以看出，當 λ1和 λ2同時增大時，其相關系數都顯著下降，因為有大量真正表達的剪切異構體被估計成未表達。而 λ1和 λ2在取值1 附近能獲得較為穩定的結果，因此本文中所有實驗都是設定 λ1和 λ2為1。

圖7 參數 λ1和 λ2 對剪切異構體表達水平的影響Fig.7 Effect of isoform expression levels by parameters λ1and λ2

3 結束語

本文提出了一個基于多條件多樣本RNASeq 測序數據的剪切異構體表達水平估計方法。為了考慮基因讀段分布在不同條件下的高度相似性，MCMS-Seq 方法設計一個聯合多條件多樣本的偏差估計模型，同時考慮了基因讀段分布的全局偏差和局部偏差所帶來的影響。從數據分析可以看出，該偏差估計模型能較為準確地描述出基因讀段非均勻分布特性。此外，MCMS-Seq 方法增加了L2/L1組稀疏約束和L1稀疏約束兩個正則化項，體現了基因和剪切異構體之間存在稀疏的生物特性，同時消除了技術性誤差和數據噪聲的影響。在小鼠數據集上，MCMS-Seq 方法估計的剪切異構體表達水平能獲得更好的可重復性。通過與SEQC 數據集中qRT-PCR 剪切異構體和MAQC-II 數據中qRT-PCR 基因的驗證，MCMSSeq 方法比其他3 個對比方法更佳的性能。

由于大量多條件多樣本數據集是時序數據集，蘊含了時間信息，但是MCMS-Seq 模型未考慮到數據中的時間信息。在未來的研究中，可以考慮在模型中融入時間信息，從而進一步提高剪切異構體的表達水平的準確性。此外，可將MCMS-Seq 模型推廣到單細胞測序數據分析，可提供更好的生物解釋性。