多條件多樣本RNA-Seq 數(shù)據(jù)的剪切異構(gòu)體表達(dá)水平估計(jì)

張禮，馬越，吳東洋

（1.南京林業(yè)大學(xué) 信息科學(xué)技術(shù)學(xué)院，江蘇 南京 210016;2.江蘇健康衛(wèi)生職業(yè)學(xué)院 中西醫(yī)結(jié)合學(xué)院，江蘇 南京 210018）

選擇性剪切事件是導(dǎo)致生物體多樣性的重要原因之一。為了進(jìn)一步揭示選擇性剪切的內(nèi)在機(jī)制，迫切需要計(jì)算剪切異構(gòu)體的表達(dá)水平。與傳統(tǒng)的基因芯片技術(shù)相比，高通量RNA 測序(RNA sequencing，RNA-Seq)技術(shù)具有高通量、高靈敏度、可重復(fù)性好等優(yōu)勢，已成為轉(zhuǎn)錄組學(xué)分析的一個標(biāo)準(zhǔn)技術(shù)手段[1-5]。

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

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

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

1 MCMS-Seq 方法

1.1 MCMS-Seq 模型表示

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

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

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

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

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

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

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

式中D表示歸一化后的數(shù)據(jù)矩陣。

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

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

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

圖2 MCMS-Seq 方法的優(yōu)化問題Fig.2 Optimization problem of MCMS-Seq method

1.2 多條件多樣本偏差估計(jì)模型

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

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

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

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

1.3 MCMS-Seq 模型實(shí)現(xiàn)

MCMS-Seq 方法的實(shí)現(xiàn)可以分為3 個部分：讀段數(shù)據(jù)預(yù)處理、基因偏差估計(jì)和表達(dá)水平估計(jì)。

1) 讀段數(shù)據(jù)預(yù)處理，是從匹配成功的讀段數(shù)據(jù)中統(tǒng)計(jì)基因每個外顯子片段的讀段計(jì)數(shù)，以及從注釋文件中獲得外顯子片段和剪切異構(gòu)體之間的映射關(guān)系矩陣。

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

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

MCMS-Seq 方法的詳細(xì)流程如算法1 所示，采用Python 和MATLAB 混合編程實(shí)現(xiàn)。

算法1MCMS-Seq 方法

輸入多條件多樣本數(shù)據(jù)，注釋文件；

輸出每個基因的剪切異構(gòu)體表達(dá)水平。

1)數(shù)據(jù)預(yù)處理：統(tǒng)計(jì)外顯子片段讀段數(shù)目矩陣D，構(gòu)建映射關(guān)系矩陣A。

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

3)表達(dá)水平估計(jì)：計(jì)算所有基因的X?。

1.4 多條件多樣本數(shù)據(jù)分析通道

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

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

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

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

2.1 數(shù)據(jù)集

3 個多條件多樣本的RNA-Seq 數(shù)據(jù)集被用來驗(yàn)證MCMS-Seq 方法估計(jì)剪切異構(gòu)體表達(dá)水平的準(zhǔn)確性。3 個數(shù)據(jù)集分別是小鼠數(shù)據(jù)集、人類大腦的SEQC 和MAQC-II 數(shù)據(jù)集，它們都來自Illumina/solexa 測序平臺。

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

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

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

2.2 多條件多樣本偏差估計(jì)模型驗(yàn)證

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

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

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

2.3 多條件多樣本下剪切異構(gòu)體表達(dá)水平的驗(yàn)證

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

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

表1 Cufflinks 估計(jì)基因Nhp2 中3 個剪切異構(gòu)體表達(dá)水平Table 1 Expression level of three isoforms in Nph2 gene estimated using cufflinks

表2 XAEM 估計(jì)基因Nhp2 中3 個剪切異構(gòu)體表達(dá)水平Table 2 Expression level of three isoforms in Nph2 gene estimated using XAEM

表3 MCMS-Seq 估計(jì)基因Nhp2 中3 個剪切異構(gòu)體表達(dá)水平Table 3 Expression level of three isoforms in Nph2 gene estimated using MCMS-Seq

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

表4 PGSeq 估計(jì)基因Nhp2 中3 個剪切異構(gòu)體表達(dá)水平Table 4 Expression level of three isoforms in Nph2 gene estimated using PGSeq

2.4 剪切異構(gòu)體的可重復(fù)性驗(yàn)證

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

表5 在小鼠數(shù)據(jù)集上不同方法估計(jì)的剪切異構(gòu)體表達(dá)水平在樣本之間的相關(guān)系數(shù)Table 5 Correlation coefficients between isoform expression levels estimated using various methods in the mouse dataset

2.5 PCR 剪切異構(gòu)體的表達(dá)水平驗(yàn)證

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

表6 在 SEQC 數(shù)據(jù)集上不同方法與 qRT-PCR 驗(yàn)證剪切異構(gòu)體之間的相關(guān)系數(shù)Table 6 Correlation coefficients between qRT-PCR values and isoform expression levels estimated using various methods in SEQC dataset

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

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

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

2.7 稀疏參數(shù)的選擇

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

圖7 參數(shù) λ1和 λ2 對剪切異構(gòu)體表達(dá)水平的影響Fig.7 Effect of isoform expression levels by parameters λ1and λ2

3 結(jié)束語

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

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