保留非全長(zhǎng)讀段的ISO-seq數(shù)據(jù)轉(zhuǎn)錄組表達(dá)分析

劉學(xué)軍 瞿錫垚 張 禮

（1.南京航空航天大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，南京，211106；2.南京林業(yè)大學(xué)信息科學(xué)技術(shù)學(xué)院，南京，210037）

引 言

由于選擇性剪切（Alternative splicing，AS）的存在，一條前體mRNA可以剪切出多條mRNA并指導(dǎo)蛋白質(zhì)的生成，這種轉(zhuǎn)錄組一對(duì)多的剪切方式是造成生物多樣性的最重要原因之一。而各項(xiàng)研究都表明選擇性剪切現(xiàn)象普遍存在于高等真核生物中[1]，同時(shí)一些特異剪切方式產(chǎn)生的新型異構(gòu)體也是導(dǎo)致基因疾病的重要原因之一，因此研究選擇性剪切對(duì)于揭示人類疾病機(jī)制具有重要的意義。

基因及異構(gòu)體表達(dá)水平計(jì)算是研究選擇性剪切的一種重要途徑，具有高通量特性的第二代測(cè)序技術(shù)RNA-seq擁有量化轉(zhuǎn)錄片段的突出優(yōu)勢(shì)，在這一領(lǐng)域具有較多有效的應(yīng)用[2-3]，很多方法采用RNA-seq數(shù)據(jù)計(jì)算基因以及異構(gòu)體的表達(dá)水平。例如，基于泊松分布的PGseq[4]和NURD(Non-uniform read distribution)[5]，以及基于讀段產(chǎn)生式的Cufflinks[6]和RESM[7-8]等。但是由于讀段長(zhǎng)度短、GC誤差等缺點(diǎn)的存在使得RNA-seq技術(shù)在識(shí)別全長(zhǎng)異構(gòu)體方面顯得十分乏力；另外，在表達(dá)水平計(jì)算方面，RNA-seq技術(shù)也面臨著讀段多源映射的問(wèn)題（參考序列中大量重復(fù)和同源序列的存在，導(dǎo)致一個(gè)讀段映射至多個(gè)位置）[9]，在已有的大多數(shù)方法中，都存在嚴(yán)重依賴注釋信息的情況，而注釋信息不完善降低了表達(dá)水平計(jì)算的準(zhǔn)確度。由于統(tǒng)計(jì)每個(gè)外顯子上讀段數(shù)量是基于RNA-seq技術(shù)計(jì)算基因及異構(gòu)體表達(dá)水平的基礎(chǔ)[10]，因此即使采用相同注釋信息的模型，不同方法在估計(jì)表達(dá)值時(shí)也存在較大差異[11]。

近幾年誕生的第三代測(cè)序技術(shù)以其讀段長(zhǎng)度長(zhǎng)、無(wú)聚合酶鏈?zhǔn)椒磻?yīng)（Polymerase chain reaction,PCR）過(guò)程引入的誤差等特點(diǎn)，迅速得到研究者的關(guān)注，并應(yīng)用于RNA-seq技術(shù)不適用的場(chǎng)景。ISO-seq技術(shù)是PacBio公司開(kāi)發(fā)的用于轉(zhuǎn)錄組研究的第三代測(cè)序技術(shù)，該技術(shù)從細(xì)胞中分離出mRNA，在size-selection（長(zhǎng)度選擇）之后制備成ISO-seq文庫(kù)，用于測(cè)序儀測(cè)序。整個(gè)測(cè)序過(guò)程沒(méi)有對(duì)測(cè)序片段作任何打斷處理，這樣的測(cè)序結(jié)果可認(rèn)為是測(cè)序片段的完整讀段，因此ISO-seq技術(shù)在誕生之后，被大多數(shù)研究者應(yīng)用于轉(zhuǎn)錄組重構(gòu)和基因組組裝等領(lǐng)域[12-13]。但該技術(shù)的測(cè)序結(jié)果存在較高的錯(cuò)誤率，為了解決這一問(wèn)題，大多數(shù)方法同時(shí)使用RNA-seq和ISO-seq數(shù)據(jù)。一方面利用RNA-seq數(shù)據(jù)的準(zhǔn)確性進(jìn)行ISO-seq數(shù)據(jù)糾錯(cuò)，另一方面利用RNA-seq數(shù)據(jù)的高通量輔助預(yù)測(cè)異構(gòu)體并計(jì)算異構(gòu)體達(dá)水平。例如，IDP[14]方法是ISO-seq數(shù)據(jù)處理的代表方法，它使用混合策略，以聚類方式從RNA-seq數(shù)據(jù)和注釋庫(kù)中找出junction（外顯子剪切點(diǎn)），去除沒(méi)有junction支持的ISO-seq長(zhǎng)讀段，將得到的非冗余多外顯子長(zhǎng)讀段作為預(yù)測(cè)異構(gòu)體，再將RNA-seq數(shù)據(jù)比對(duì)至預(yù)測(cè)異構(gòu)體并計(jì)算各個(gè)預(yù)測(cè)異構(gòu)體的表達(dá)水平。

IDP在除去冗余全長(zhǎng)讀段時(shí)，不僅去掉了信息一致的全長(zhǎng)讀段，還刪除了全長(zhǎng)讀段之間的包含情況，即當(dāng)一條較短的全長(zhǎng)讀段包含于一條較長(zhǎng)全長(zhǎng)讀段時(shí)，只保留較長(zhǎng)的讀段。但從ISO-seq測(cè)序技術(shù)原理來(lái)看，在制備cDNA文庫(kù)時(shí)并沒(méi)有進(jìn)行類似RNA-seq測(cè)序技術(shù)的隨機(jī)打斷，因此本文認(rèn)為一條全長(zhǎng)讀段等價(jià)一個(gè)異構(gòu)體，去掉包含關(guān)系的讀段會(huì)遺漏異構(gòu)體，最終影響到異構(gòu)體的預(yù)測(cè)結(jié)果。另外，大多數(shù)研究工作認(rèn)為ISO-seq數(shù)據(jù)的通量低，不適合計(jì)算表達(dá)水平，其主要原因在于這些研究工作大都丟棄了占到所有數(shù)據(jù)50%～60%的非全長(zhǎng)讀段[15]。而非全長(zhǎng)讀段的產(chǎn)生是由于測(cè)序時(shí)酶失活，導(dǎo)致測(cè)序過(guò)程無(wú)法繼續(xù)進(jìn)行，但非全長(zhǎng)讀段仍然具有ISO-seq數(shù)據(jù)超長(zhǎng)讀長(zhǎng)的特性，也包含外顯子信息，能夠反映樣本中轉(zhuǎn)錄本的濃度，丟棄非全長(zhǎng)讀段將會(huì)直接影響基因和異構(gòu)體表達(dá)水平的計(jì)算，故目前的方法大多采用ISO-seq數(shù)據(jù)和RNA-seq數(shù)據(jù)相結(jié)合的方式進(jìn)行異構(gòu)體表達(dá)水平的計(jì)算。

本文首次提出僅利用ISO-seq數(shù)據(jù)，且保留非全長(zhǎng)讀段進(jìn)行基于狄利克雷采樣的探測(cè)與預(yù)測(cè)（Dirichlet sampling for isoform detection and prediction，DSIDP）方法。同時(shí)，第三代測(cè)序技術(shù)雖然擁有超長(zhǎng)讀長(zhǎng)測(cè)序，但也無(wú)法保證全長(zhǎng)讀段數(shù)據(jù)涵蓋所有表達(dá)異構(gòu)體，針對(duì)一些沒(méi)有全長(zhǎng)讀段數(shù)據(jù)的異構(gòu)體預(yù)測(cè)問(wèn)題，本文在沿用DSIDP預(yù)測(cè)異構(gòu)體思想的基礎(chǔ)之上，還提出了一種基于馬爾科夫鏈的異構(gòu)體探測(cè)與預(yù)測(cè)（Markov chain for isoform detection and prediction,MCIDP）方法。兩種模型均在模擬數(shù)據(jù)集和真實(shí)數(shù)據(jù)上得到了有效驗(yàn)證。

1 實(shí)驗(yàn)方法

1.1 數(shù)據(jù)特性

圖1顯示了ISO-seq數(shù)據(jù)中全長(zhǎng)讀段和非全長(zhǎng)讀段長(zhǎng)度分布直方圖，圖中數(shù)據(jù)來(lái)自PacBio公司公開(kāi)數(shù)據(jù)集 MCF-7（http：//www.pacb.com/blog/data-release-human-mcf-7-transcriptome/）。本文統(tǒng)計(jì)了 6個(gè)cell的原始數(shù)據(jù)（如表1所示），其中按照ISO-seq技術(shù)的size-selection原則，對(duì)樣本長(zhǎng)度1～2 Kb，2～3 Kb和＞3 Kb三個(gè)范圍的Cell各選取兩個(gè)。從統(tǒng)計(jì)結(jié)果可以看出全長(zhǎng)讀段和非全長(zhǎng)讀段的長(zhǎng)度分布具有相似的模態(tài)，數(shù)據(jù)多集中在長(zhǎng)度為1～3 Kb的區(qū)間內(nèi)，這說(shuō)明非全長(zhǎng)讀段數(shù)據(jù)也具有遠(yuǎn)超過(guò)RNA-seq數(shù)據(jù)的長(zhǎng)度，從第三代測(cè)序數(shù)據(jù)超長(zhǎng)讀長(zhǎng)這一本質(zhì)特征來(lái)說(shuō)，非全長(zhǎng)讀段與全長(zhǎng)讀段一樣，也包含關(guān)于異構(gòu)體的有效信息。并且隨著樣本序列長(zhǎng)度的增加，非全長(zhǎng)讀段也隨之增加，并達(dá)到接近60%。如果在異構(gòu)體的構(gòu)建中不考慮這部分?jǐn)?shù)據(jù)，相當(dāng)于丟棄了大部分實(shí)驗(yàn)數(shù)據(jù)。因此將ISO-seq數(shù)據(jù)應(yīng)用于轉(zhuǎn)錄組學(xué)研究領(lǐng)域時(shí)，保留非全長(zhǎng)讀段具有重要意義。

圖1 ISO-seq數(shù)據(jù)中讀段長(zhǎng)度分布直方圖Fig.1 Histograms of ISO-seq reads

表1 MCF-7數(shù)據(jù)集讀段統(tǒng)計(jì)Tab.1 Read statistics of MCF-7 data set

為了進(jìn)一步說(shuō)明ISO-seq數(shù)據(jù)用于計(jì)算表達(dá)水平的可行性，隨機(jī)選擇了4個(gè)基因并統(tǒng)計(jì)它們?cè)贗SO-seq數(shù)據(jù)和RNA-seq數(shù)據(jù)中外顯子上讀段數(shù)分布情況，結(jié)果如圖2所示，第一行顯示了RNA-seq數(shù)據(jù)結(jié)果，第二行顯示了ISO-seq數(shù)據(jù)結(jié)果。從圖2中可以看出ISO-seq數(shù)據(jù)與RNA-seq數(shù)據(jù)具有極為相似的分布模態(tài)，這說(shuō)明ISO-seq數(shù)據(jù)與RNA-seq數(shù)據(jù)一樣，能通過(guò)讀段計(jì)數(shù)反映出樣本中相應(yīng)轉(zhuǎn)錄本的濃度，進(jìn)而可以用來(lái)計(jì)算轉(zhuǎn)錄組的表達(dá)水平。

圖2 4個(gè)基因在ISO-seq數(shù)據(jù)和RNA-seq數(shù)據(jù)中外顯子上的讀段分布Fig.2 Distribution of reads on exons of four genes in ISO-seq data and RNA-seq data

1.2 原始數(shù)據(jù)處理

ISO-seq的下機(jī)數(shù)據(jù)中全長(zhǎng)和非全長(zhǎng)讀段是混合在一起的，需要按照一定的準(zhǔn)則將其區(qū)分開(kāi)。PacBio公司提供的SMRT Analysis軟件根據(jù)讀段數(shù)據(jù)兩端是否均存在接頭序列，將其分為全長(zhǎng)讀段和非全長(zhǎng)讀段，但這樣的方式需要原始讀段數(shù)據(jù)之間的比對(duì)，對(duì)計(jì)算效率影響較大。圖3展示了PacBio測(cè)序原理，ISO-seq數(shù)據(jù)的cDNA文庫(kù)是兩端接上接頭的啞鈴狀結(jié)構(gòu)，測(cè)序時(shí)會(huì)在整個(gè)結(jié)構(gòu)上循環(huán)進(jìn)行[16]。根據(jù)這樣的測(cè)序原理，本文使用一種簡(jiǎn)單高效的方法區(qū)分全長(zhǎng)和非全長(zhǎng)讀段。當(dāng)一個(gè)零膜導(dǎo)波管（Zero mode waveguide，ZMW）中出現(xiàn)多條Subread時(shí)，從中選擇最長(zhǎng)的讀段劃分至全長(zhǎng)讀段集合，否則，將唯一的一條Subread劃分至非全長(zhǎng)讀段集合。整個(gè)過(guò)程不僅區(qū)分出了全長(zhǎng)和非全長(zhǎng)讀段，還去掉了多條Subread中的冗余讀段。

圖3 PacBio測(cè)序原理Fig.3 PacBio sequencing principle

ISO-seq數(shù)據(jù)的另一特征是較高的測(cè)序錯(cuò)誤率，目前的大多數(shù)研究工作均使用RNA-seq數(shù)據(jù)對(duì)其進(jìn)行糾正，本文使用LoRDEC[17]對(duì)非冗余讀段（包括全長(zhǎng)讀段和非全長(zhǎng)讀段）進(jìn)行糾錯(cuò)處理，參數(shù)值k設(shè)置為21，s設(shè)置為3。糾錯(cuò)后的數(shù)據(jù)使用BWA-MEM[18]比對(duì)至參考基因組，借助基因注釋庫(kù)信息從比對(duì)結(jié)果中獲取到讀段中的外顯子序列。對(duì)于沒(méi)有RNA-seq數(shù)據(jù)的情況，同樣可以采用僅使用ISO-seq數(shù)據(jù)進(jìn)行自糾錯(cuò)的方法，例如Chen等[19]使用最長(zhǎng)的讀段作為種子來(lái)收集其他所有讀段，構(gòu)建高度準(zhǔn)確的讀段數(shù)據(jù)。

1.3 真實(shí)數(shù)據(jù)有效性驗(yàn)證

MCF-7數(shù)據(jù)集中共有119個(gè)Cell，且測(cè)序時(shí)間并不都是一致，因此本文在建模之前驗(yàn)證了數(shù)據(jù)的有效性。選取6個(gè)Cell的真實(shí)數(shù)據(jù)，分為兩組，每組均包含Size-selection的3個(gè)長(zhǎng)度范圍，即Cell 1,Cell 3,Cell 5為一組記為Group 1，Cell 2,Cell 4,Cell 6為另一組記為Group 2。對(duì)兩組數(shù)據(jù)中的5 665個(gè)公共基因通過(guò)計(jì)算RPKM[20]值得到表達(dá)水平，在對(duì)數(shù)刻度上驗(yàn)證兩組數(shù)據(jù)獲得的基因表達(dá)水平的吻合性。結(jié)果如圖4所示，相關(guān)系數(shù)為0.900 6。可以看出這兩組重復(fù)實(shí)驗(yàn)在基因表達(dá)值上具有很高的一致性，表明多次測(cè)量得到的讀段具有較好的可重復(fù)性，因此數(shù)據(jù)集的有效性得以驗(yàn)證。

圖4 MCF-7數(shù)據(jù)集重復(fù)實(shí)驗(yàn)結(jié)果對(duì)比Fig.4 Comparison of repeated experiment results of MCF-7 data set

1.4 DSIDP模型

PacBio測(cè)序技術(shù)從細(xì)胞中提取到mRNA后沒(méi)有進(jìn)行分子隨機(jī)打斷，因此本文認(rèn)為經(jīng)過(guò)對(duì)下機(jī)原始數(shù)據(jù)的處理后，得到的所有非冗余全長(zhǎng)讀段數(shù)據(jù)即為細(xì)胞中表達(dá)的異構(gòu)體集合，并將這個(gè)集合作為模型的異構(gòu)體預(yù)測(cè)結(jié)果。

借鑒RNA-seq數(shù)據(jù)計(jì)算表達(dá)水平時(shí)統(tǒng)計(jì)外顯子上讀段數(shù)量的方式，本文將所有ISO-seq讀段數(shù)據(jù)映射至異構(gòu)體預(yù)測(cè)集，并統(tǒng)計(jì)各預(yù)測(cè)異構(gòu)體上讀段數(shù)量，在總讀段數(shù)上做歸一化得到其表達(dá)水平。映射過(guò)程中，ISO-seq數(shù)據(jù)也將面臨相同基因下多個(gè)異構(gòu)體之間的多源映射問(wèn)題。MCF-7數(shù)據(jù)集的6個(gè)Cell中，35%的讀段存在多源映射的情況，這遠(yuǎn)低于二代數(shù)據(jù)70%讀段的多源映射情況[9]。這里的ISO-seq數(shù)據(jù)多源映射是指一條非全長(zhǎng)讀段映射至多條預(yù)測(cè)異構(gòu)體，如何分配這樣的非全長(zhǎng)讀段是計(jì)算異構(gòu)體表達(dá)水平中要解決的核心問(wèn)題。為了解決這個(gè)問(wèn)題，本文提出了DSIDP模型。DSIDP是一個(gè)基于Dirichlet分布，對(duì)該問(wèn)題進(jìn)行建模求解的算法，使用隨機(jī)采樣方法將發(fā)生多源映射的非全長(zhǎng)讀段映射到概率最大的異構(gòu)體,通過(guò)這樣的方式利用非全長(zhǎng)讀段進(jìn)行異構(gòu)體表達(dá)比例的計(jì)算。具體算法過(guò)程如下：

算法1DSIDP

輸入：全長(zhǎng)讀段數(shù)據(jù)XFL，非全長(zhǎng)讀段數(shù)據(jù)XnFL以及異構(gòu)體預(yù)測(cè)集合T，和均代表一條讀段數(shù)據(jù)，Ti代表一個(gè)異構(gòu)體，|T|=k；

輸出：預(yù)測(cè)異構(gòu)體的表達(dá)水平向量E，每一維代表相應(yīng)預(yù)測(cè)異構(gòu)體的表達(dá)水平。

(1)將讀段數(shù)據(jù)矩陣XFL和XnFL映射至異構(gòu)體預(yù)測(cè)集合矩陣T，統(tǒng)計(jì)每個(gè)異構(gòu)體上唯一映射讀段數(shù)量，得到一個(gè)k維向量，并對(duì)每一維在所有維度上做歸一化記為τ。

(2)將發(fā)生多源映射的讀段數(shù)據(jù)合并記為Xm，則每一個(gè)均對(duì)應(yīng)一個(gè)tj(tj?T)和一個(gè)τj(τj?τ)，τj是歸一化的結(jié)果，Isoform～Dirichlet(τj)。

(3)從Isoform～Dirichlet(τj)中采樣得到變量isoform，其中各維度上的概率值表示屬于對(duì)應(yīng)異構(gòu)體的可能性，選擇概率最大的異構(gòu)體，在其讀段計(jì)數(shù)上加一。

(4)遍歷完所有Xmj之后得到新的異構(gòu)體讀段計(jì)數(shù)向量，歸一化處理結(jié)果記為表達(dá)水平E。

1.5 MCIDP模型

在RNA-seq數(shù)據(jù)的處理中，LeGault等[21]使用概率連接圖（Probabilistic splice graphs,PSGs）方法對(duì)異構(gòu)體結(jié)構(gòu)進(jìn)行預(yù)測(cè)。在固定基因結(jié)構(gòu)的情況下，量化基因的選擇性剪接事件，從測(cè)序數(shù)據(jù)中找到異構(gòu)體的junction，通過(guò)junction之間的跳轉(zhuǎn)做出異構(gòu)體結(jié)構(gòu)的預(yù)測(cè)。本文將這樣的思想運(yùn)用到第3代測(cè)序數(shù)據(jù)中，提出了MCIDP模型。由于ISO-seq數(shù)據(jù)長(zhǎng)讀段的特點(diǎn)，在junction跨越很長(zhǎng)的區(qū)域時(shí)，也能有讀段的支持，因此這樣找到的junction較之LeGault等使用RNA-seq數(shù)據(jù)要更為精確和全面。

MCIDP使用馬爾科夫鏈對(duì)異構(gòu)體junction之間的跳轉(zhuǎn)進(jìn)行建模。一個(gè)基本的馬爾科夫鏈包含3元素：狀態(tài)節(jié)點(diǎn)(V)、初始狀態(tài)概率向量(π)、狀態(tài)轉(zhuǎn)移概率矩陣(A)，因此，模型可以表示為G=(V,A,π)。其中狀態(tài)節(jié)點(diǎn)由基因結(jié)構(gòu)決定，將基因外顯子由編號(hào)從小到大進(jìn)行排列，并在該排列的兩端加上起始點(diǎn)V0=0和終止點(diǎn)VM=M(M=|V|-1)，即為狀態(tài)節(jié)點(diǎn)集合；Aij表示狀態(tài)節(jié)點(diǎn)i轉(zhuǎn)移至狀態(tài)節(jié)點(diǎn)j的概率值，且有Aij∈[0,1],?i,；πi表示由狀態(tài)節(jié)點(diǎn)i作為路徑起始點(diǎn)的概率值，且有從模型建立的整個(gè)過(guò)程可以看出，MCIDP方法只需要知道基因外顯子組成，不依賴注釋庫(kù)中異構(gòu)體的注釋信息。圖5為MCIDP建模示意圖。

圖5 MCIDP建模示意圖Fig.5 Modeling diagram of MCIDP

模型的一條路徑σt即代表一個(gè)可能存在的異構(gòu)體，例如圖5中的路徑V0V1V3V5V6。在任意一條沒(méi)有全長(zhǎng)讀段的路徑上，如果存在其他的讀段能夠拼接出該路徑，那么該路徑也能夠被模型模擬并預(yù)測(cè)。表示路徑t中的第k個(gè)狀態(tài)節(jié)點(diǎn)，根據(jù)之前的設(shè)定，=V0,=VM，其中e=|σ|，表示路徑所包含的狀態(tài)節(jié)點(diǎn)數(shù),下標(biāo)指示路徑節(jié)點(diǎn)，上標(biāo)指代路徑。路徑中轉(zhuǎn)移概率的累積乘積Expr(σt)可表示為

使用極大似然估計(jì)出馬爾科夫鏈模型的參數(shù)π和A，令Nij表示狀態(tài)節(jié)點(diǎn)i與狀態(tài)節(jié)點(diǎn)j之間的junction總個(gè)數(shù)，則對(duì)于Aij和πi的極大似然估計(jì)有

MCIDP沿用了DSIDP從全長(zhǎng)讀段中建立異構(gòu)體預(yù)測(cè)集的思想，將所有非冗余全長(zhǎng)讀段數(shù)據(jù)作為模型預(yù)測(cè)異構(gòu)體的初始集合。由于構(gòu)造的圖模型中有些路徑的junction結(jié)構(gòu)較為相似，可以進(jìn)行合并計(jì)算，所以需對(duì)所有其他可能存在的異構(gòu)體根據(jù)定義的距離公式，將其路徑概率累加到結(jié)構(gòu)最近的預(yù)測(cè)異構(gòu)體中。這里距離公式的定義同時(shí)考慮到了兩個(gè)junction之間局部跨區(qū)域的差異和所有junction之間累積起來(lái)的全局差異，具體描述如下：

令St表示一個(gè)異構(gòu)體的外顯子序列，基因的第i個(gè)外顯子包含于其中，則=1，否則=0,|St|=M。對(duì)于兩個(gè)異構(gòu)體外顯子序列St1,St2，若＜i＜M，則認(rèn)為節(jié)點(diǎn)i處存在這兩個(gè)異構(gòu)體的相似junction，且為開(kāi)始位置，若=,i＜j＜M，則認(rèn)為節(jié)點(diǎn)j處為該相似junction的結(jié)束位置。由此，兩個(gè)異構(gòu)體中相似junction的初步距離定義為J(St1,St2),可表示為

式中λJ為度量因子，作為指數(shù)距離公式中的底數(shù)?？紤]到差異區(qū)域長(zhǎng)度對(duì)距離的影響，令li表示基因第i個(gè)外顯子的長(zhǎng)度，L(St)表示異構(gòu)體St的長(zhǎng)度，則兩個(gè)junction的差異區(qū)域長(zhǎng)度對(duì)距離的影響可以定義為Iij(St1,St2)，可表示為

式中λI可視為懲罰因子，所以兩個(gè)異構(gòu)體中相似junction的最終距離定義為Dij(St1,St2)，可表示為

則兩個(gè)可能存在的異構(gòu)體的距離D(St1,St2)可表示為

對(duì)于超出距離閾值的可能存在的異構(gòu)體，對(duì)其作Kmeans聚類處理，距離公式使用式（6），并將聚類中心作為新的預(yù)測(cè)異構(gòu)體添加至異構(gòu)體預(yù)測(cè)集合，該預(yù)測(cè)異構(gòu)體的表達(dá)比率等價(jià)于以其為聚類中心的所有可能存在的異構(gòu)體路徑概率之和。最終，模型輸出異構(gòu)體預(yù)測(cè)集合及集合元素各自的概率值，該概率值即為該基因每個(gè)可能存在的異構(gòu)體的表達(dá)比例。

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

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

本文使用了一個(gè)模擬數(shù)據(jù)集和一個(gè)真實(shí)數(shù)據(jù)集來(lái)驗(yàn)證兩個(gè)模型的有效性。模擬數(shù)據(jù)集中，假設(shè)了一個(gè)擁有10個(gè)外顯子和4個(gè)異構(gòu)體的基因，并設(shè)置異構(gòu)體的表達(dá)比例分別為t1=0.3，t2=0.3，t3=0.2和t4=0.2，如圖6所示。按照設(shè)定的比例，采樣生成了100條全長(zhǎng)讀段數(shù)據(jù)，根據(jù)后續(xù)實(shí)驗(yàn)的需求從中隨機(jī)選取n條全長(zhǎng)讀段，作隨機(jī)打斷處理，生成非全長(zhǎng)讀段。真實(shí)數(shù)據(jù)集來(lái)自PacBio公開(kāi)數(shù)據(jù)MCF-7，本文選取了其中6個(gè)Cell的數(shù)據(jù)，各Cell數(shù)據(jù)讀段的統(tǒng)計(jì)情況見(jiàn)表1。

圖6 模擬數(shù)據(jù)結(jié)構(gòu)Fig.6 Structure of simulation data

2.2 非全長(zhǎng)讀段有效性驗(yàn)證

從表1可以看出當(dāng)讀段長(zhǎng)度越長(zhǎng)時(shí)，非全長(zhǎng)讀段的數(shù)量就越多。因此本文在模擬數(shù)據(jù)集上做了非全長(zhǎng)讀段不同占比的對(duì)照實(shí)驗(yàn)，將100條全長(zhǎng)讀段按25%，50%和75%的比例隨機(jī)抽取，并作隨機(jī)打斷，產(chǎn)生相應(yīng)比例的非全長(zhǎng)讀段，剩下的全長(zhǎng)讀段作為對(duì)照組，全長(zhǎng)讀段加上非全長(zhǎng)讀段作為實(shí)驗(yàn)組。在各比例的對(duì)照實(shí)驗(yàn)中，實(shí)驗(yàn)組與對(duì)照組均使用DSIDP方法計(jì)算結(jié)果，并采用計(jì)算值與真實(shí)值之間的歐式距離作為誤差度量。如表2和表3所示，在加入非全長(zhǎng)讀段數(shù)據(jù)后，各比例實(shí)驗(yàn)的表達(dá)水平計(jì)算值均比只使用全長(zhǎng)讀段數(shù)據(jù)更為精確，表中FL(Full length)表示全長(zhǎng)讀段，nFL（non-full length）表示非全長(zhǎng)讀段。值得指出的是，在非全長(zhǎng)讀段數(shù)據(jù)占75%的比例時(shí)，誤差有大幅度的下降，但誤差本身仍然要比其他比例只用全長(zhǎng)讀段數(shù)據(jù)結(jié)果值大，這說(shuō)明在計(jì)算異構(gòu)體表達(dá)水平時(shí)，保留非全長(zhǎng)讀段數(shù)據(jù)能夠降低只使用全長(zhǎng)讀段數(shù)據(jù)的計(jì)算誤差。另外，模擬數(shù)據(jù)集構(gòu)建的假設(shè)前提是該基因的所有異構(gòu)體均來(lái)自細(xì)胞內(nèi)當(dāng)前表達(dá)且被測(cè)序到的mRNA分子，與注釋庫(kù)中的信息無(wú)關(guān)，因此可以認(rèn)為當(dāng)細(xì)胞內(nèi)出現(xiàn)新型異構(gòu)體時(shí)，也能被DSIDP預(yù)測(cè)出。例如，假設(shè)t4為新型異構(gòu)體，且100條讀段數(shù)據(jù)中包含有t4，則會(huì)被DSIDP預(yù)測(cè)出其結(jié)構(gòu)和表達(dá)值。

表2 模擬數(shù)據(jù)各比例非全長(zhǎng)讀段計(jì)算結(jié)果Tab.2 Calculation results on simulation data with different nFL read proportions

表3 模擬數(shù)據(jù)各比例非全長(zhǎng)讀段計(jì)算誤差Tab.3 Calculation error on simulation data with different nFL read proportions

2.3 MCIDP預(yù)測(cè)異構(gòu)體驗(yàn)證

MCIDP的提出是為了預(yù)測(cè)出數(shù)據(jù)中沒(méi)有全長(zhǎng)讀段的超長(zhǎng)異構(gòu)體，本文將模擬數(shù)據(jù)中t1異構(gòu)體的所有全長(zhǎng)讀段隨機(jī)打斷，這時(shí)t1異構(gòu)體即可作為沒(méi)有全長(zhǎng)讀段的超長(zhǎng)異構(gòu)體，檢驗(yàn)?zāi)Ｐ偷念A(yù)測(cè)能力。實(shí)驗(yàn)結(jié)果如表4所示，可以看出模型能預(yù)測(cè)出t1這樣的超長(zhǎng)異構(gòu)體，但在表達(dá)水平計(jì)算上，DSIDP要比MCIDP更精確，原因在于基于馬爾科夫鏈的MCIDP會(huì)產(chǎn)生較多低概率的可能路徑。如何把這些低概率路徑合并至真實(shí)異構(gòu)體中是該類模型后續(xù)研究的一個(gè)重點(diǎn)。

表4 MCIDP在模擬數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果Tab.4 MCIDP results on simulation data

2.4 真實(shí)數(shù)據(jù)集實(shí)驗(yàn)結(jié)果

在異構(gòu)體表達(dá)水平上，雖然ISO-seq數(shù)據(jù)和RNA-seq數(shù)據(jù)均反映出樣本中原始轉(zhuǎn)錄本的濃度，但是由于測(cè)序技術(shù)本身和數(shù)據(jù)特性的較大差異，尤其讀段長(zhǎng)度的差異導(dǎo)致異構(gòu)體構(gòu)建上的明顯差別，造成兩種數(shù)據(jù)在異構(gòu)體表達(dá)比例計(jì)算上的不一致，故無(wú)法采用RNA-seq數(shù)據(jù)的計(jì)算結(jié)果對(duì)ISO-seq分析結(jié)果進(jìn)行驗(yàn)證。因此，對(duì)本文中6個(gè)cell數(shù)據(jù)進(jìn)行分組，分為兩次技術(shù)性重復(fù)實(shí)驗(yàn)，具體分組方式和1.3節(jié)中的相同。其中Group 1包含139 116個(gè)全長(zhǎng)讀段，147 190個(gè)非全長(zhǎng)讀段；Group 2包含109 969個(gè)全長(zhǎng)讀段，95 431個(gè)非全長(zhǎng)讀段。將本文提出的兩個(gè)模型應(yīng)用到這兩組重復(fù)實(shí)驗(yàn)數(shù)據(jù)中，檢驗(yàn)在公共異構(gòu)體上獲得的表達(dá)比例的吻合程度，驗(yàn)證本文方法的有效性。

表5給出了兩種方法所預(yù)測(cè)的異構(gòu)體數(shù)量以及注釋庫(kù)異構(gòu)體數(shù)量（注釋庫(kù)為GENCODE數(shù)據(jù)庫(kù)中GRCh37-mapped Releases.26），圖7展示了表5數(shù)據(jù)的韋恩圖，可以看出MCIDP預(yù)測(cè)出了更多的異構(gòu)體，與注釋庫(kù)中已有異構(gòu)體的交集也更多。因此，MCIDP較適用于注重預(yù)測(cè)異構(gòu)體數(shù)量的問(wèn)題中。

表5 模型異構(gòu)體數(shù)量預(yù)測(cè)結(jié)果Tab.5 Number of predicted isoforms

圖7 模型預(yù)測(cè)異構(gòu)體數(shù)量韋恩圖Fig.7 Venn diagram of isoform numbers predicted from various methods

另外，注釋庫(kù)是對(duì)人類基因的所有已知異構(gòu)體進(jìn)行注釋，在一些分化后的人體細(xì)胞中并不是所有基因都表達(dá)，所以圖7中注釋庫(kù)中有較多的異構(gòu)體未被兩種模型預(yù)測(cè)出，而兩個(gè)模型都預(yù)測(cè)出了大量不在注釋庫(kù)中的異構(gòu)體，這在一定程度上也說(shuō)明了當(dāng)前注釋庫(kù)還很不完善。圖8展示了所提出的兩種方法在兩次重復(fù)實(shí)驗(yàn)中計(jì)算得到的公共異構(gòu)體（共4 914個(gè)）表達(dá)比例的散點(diǎn)圖。為了更好地呈現(xiàn)大部分低比例異構(gòu)體數(shù)據(jù)的分布情況，本文采用log(105x+1)函數(shù)對(duì)異構(gòu)體表達(dá)比例進(jìn)行了變換。經(jīng)過(guò)函數(shù)變換，DSIDP結(jié)果的相關(guān)系數(shù)為0.681 7，MCIDP的為0.665 0?？梢钥闯鲈趦山M實(shí)驗(yàn)數(shù)據(jù)量不完全一致的情況下，兩個(gè)模型計(jì)算的異構(gòu)體表達(dá)比例也能有較好的一致性，這在一定程度上驗(yàn)證了本文方法的有效性。其中，DSIDP計(jì)算的異構(gòu)體比例在重復(fù)實(shí)驗(yàn)中的吻合性要高于MCIDP，顯示了其更為準(zhǔn)確的表達(dá)水平和計(jì)算能力。

圖8 MCF-7數(shù)據(jù)集異構(gòu)體層面重復(fù)實(shí)驗(yàn)結(jié)果對(duì)比Fig.8 Comparison of repeated experiment results of isoforms in MCF-7 data set

3 結(jié)束語(yǔ)

本文在保留ISO-seq數(shù)據(jù)非全長(zhǎng)讀段的基礎(chǔ)上提出了兩個(gè)適用于不同場(chǎng)景的異構(gòu)體預(yù)測(cè)和表達(dá)比例計(jì)算模型，DSIDP和MCIDP。兩個(gè)模型首次僅采用PacBio第三代測(cè)序數(shù)據(jù)用于異構(gòu)體預(yù)測(cè)以及表達(dá)水平的計(jì)算。DSIDP從全長(zhǎng)讀段中建立異構(gòu)體預(yù)測(cè)集合，將所有讀段映射至這個(gè)集合之中，統(tǒng)計(jì)集合元素各自的讀段數(shù)量，進(jìn)而計(jì)算表達(dá)水平，采用Dirichlet采樣的方法解決了讀段多源映射的問(wèn)題。實(shí)驗(yàn)結(jié)果表明DSIDP在異構(gòu)體表達(dá)水平計(jì)算上具有較好的準(zhǔn)確性。MCIDP是基于馬爾科夫鏈的一個(gè)概率模型，通過(guò)構(gòu)造概率圖模型，考慮了轉(zhuǎn)錄本中所有可能的轉(zhuǎn)錄路徑，以獲得所有可能的異構(gòu)體。在一些超長(zhǎng)異構(gòu)體無(wú)法獲得全長(zhǎng)讀段的情況下，使用MCIDP可以有效地預(yù)測(cè)出超長(zhǎng)異構(gòu)體。與IDP相比，該模型不依賴異構(gòu)體的注釋信息，只需獲取基因的外顯子組成即可預(yù)測(cè)出數(shù)據(jù)中的異構(gòu)體，但該模型在計(jì)算異構(gòu)體表達(dá)水平上具有一定不足，這與模型存在的低概率相似路徑合并這一難點(diǎn)有關(guān)。模型中使用的最近距離劃分和聚類處理，實(shí)際上都是對(duì)相似路徑的合并，在后續(xù)的工作中，擬嘗試采用二階馬爾科夫鏈模型提高相似路徑聚類的準(zhǔn)確性，以進(jìn)一步提高異構(gòu)體比例計(jì)算的準(zhǔn)確性。