摘要：可變剪接是遺傳信息傳遞與表達(dá)的重要環(huán)節(jié)，真核生物多外顯子基因pre-mRNA的可變剪接現(xiàn)象非常普遍。可變剪接增加了蛋白質(zhì)的多樣性，影響諸如細(xì)胞組織特異性分化、個(gè)體發(fā)育、疾病發(fā)生等重要的生物學(xué)過程?；趦?nèi)含子的k-mer（k=1…5）信息，利用多樣性增量結(jié)合二次判別的方法，對(duì)保留型內(nèi)含子和組成型內(nèi)含子進(jìn)行了區(qū)分。五折交叉檢驗(yàn)顯示的預(yù)測(cè)總精度、敏感性和特異性均大于70%。表明內(nèi)含子的序列信息可能對(duì)保留型內(nèi)含子的剪接有重要調(diào)控作用。

關(guān)鍵詞：可變剪接；內(nèi)含子保留；多樣性增量；二次判別函數(shù)

中圖分類號(hào)：Q61 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：0439-8114（2016）01-0208-05

DOI：10.14088/j.cnki.issn0439-8114.2016.01.055

20世紀(jì)80年代Walter Gilbert研究Adenovirus hexon基因時(shí)提出了可變剪接概念，解釋了一條RNA序列為何能編碼多個(gè)蛋白質(zhì)的問題。早期的研究認(rèn)為只有5%的高等真核生物基因存在可變剪接現(xiàn)象[1]，但隨著深度測(cè)序方法的應(yīng)用，研究發(fā)現(xiàn)超過95%的人類基因存在可變剪接現(xiàn)象[2]，其中70%～80%的可變剪接與蛋白質(zhì)的功能相關(guān)?？勺兗艚邮沁z傳信息傳遞與表達(dá)的重要環(huán)節(jié)，通過選擇不同外顯子產(chǎn)生不同的mRNA剪接變體，影響諸如細(xì)胞組織特異性分化、個(gè)體發(fā)育、疾病發(fā)生等重要的生物學(xué)過程或現(xiàn)象[3，4]。由于剪接過程的復(fù)雜性，可變剪接的調(diào)控機(jī)制仍未闡明。目前試驗(yàn)識(shí)別剪接位點(diǎn)的方法有表達(dá)序列標(biāo)簽比對(duì)[5-9]、基因芯片[10]和RNA-seq[11]等技術(shù)。采用RNA-seq等高通量技術(shù)識(shí)別和研究可變剪接顯然是很好的途徑，但是通過試驗(yàn)手段確定所有的可變剪接位點(diǎn)費(fèi)時(shí)、費(fèi)力，而發(fā)展優(yōu)良的理論模型對(duì)可變剪接位點(diǎn)進(jìn)行精確預(yù)測(cè)并探討其剪接機(jī)制既省時(shí)省力又可以指導(dǎo)試驗(yàn)工作的進(jìn)行。

目前，已經(jīng)發(fā)現(xiàn)的可變剪接形式大致可以分為5種[5，12]：盒式外顯子、可變5′剪接、可變3′剪接、內(nèi)含子保留和互斥外顯子（圖1）。研究者已經(jīng)利用基因序列的保守性等特征，綜合數(shù)學(xué)、計(jì)算機(jī)手段設(shè)計(jì)了一些可變剪接的理論預(yù)測(cè)方法，當(dāng)前已經(jīng)用于剪接位點(diǎn)預(yù)測(cè)的包括支持向量機(jī)（SVM）[13，14]、隱馬爾科夫模型（HMM）[15]、神經(jīng)網(wǎng)絡(luò)模型（NNM）[16，17]、權(quán)重矩陣（WAM）[18]及概率模型等。這些方法能夠不依賴于轉(zhuǎn)錄豐度，直接在全基因組范圍內(nèi)快速識(shí)別潛在的剪接位點(diǎn)。目前，預(yù)測(cè)剪接位點(diǎn)的工作主要針對(duì)可變盒式外顯子、可變5′和可變3′進(jìn)行，而對(duì)試驗(yàn)數(shù)據(jù)相對(duì)缺乏的內(nèi)含子保留型可變剪接的預(yù)測(cè)工作較少[19]。研究者基于結(jié)合多樣性指標(biāo)的二次判別分析法（Diversity measure combinied with quadratic discriminant analysis，IDQD）區(qū)分了保留型內(nèi)含子剪接位點(diǎn)和臨近區(qū)假剪接位點(diǎn)[20]，本工作統(tǒng)計(jì)分析了保留型內(nèi)含子序列的長(zhǎng)度，計(jì)算其k-mer（k=1～5）的頻次，并使用結(jié)合多樣性指標(biāo)的二次判別分析法（IDQD）[21，22]對(duì)保留型內(nèi)含子和組成型內(nèi)含子進(jìn)行了分類預(yù)測(cè)，為今后基于序列信息預(yù)測(cè)內(nèi)含子保留型可變剪接提供了參考。

1 材料和方法

1.1 可變剪接數(shù)據(jù)集的構(gòu)建

從ASTD數(shù)據(jù)庫(kù)中下載文件AltSplice-rel3.events.txt和AltSplice-rel3.intron.txt（ftp：//ftp.ebi.ac.uk/pub/databases/astd/）[23]。從AltSplice-rel3.events.txt中挑選出長(zhǎng)度范圍為130～519 nt的保留型內(nèi)含子（Simple Intron retention，SIR）序列2 199條，其中1 759條保留內(nèi)含子序列作為訓(xùn)練正集，440條保留內(nèi)含子作為測(cè)試正集。從AltSplice-rel3.intron.txt文件內(nèi)識(shí)別組成型內(nèi)含子（Constitutive intron，CI），剔除非GT-AG剪接類型后，共獲得3 299條與內(nèi)含子保留長(zhǎng)度范圍一致的組成型內(nèi)含子，其中2 639條作為訓(xùn)練負(fù)集，660條作為測(cè)試負(fù)集（表1）。

1.2 選取特征參數(shù)

特征參數(shù)的選取對(duì)分類器的性能是非常關(guān)鍵的，反映堿基保守性的單堿基組成和反映序列緊鄰關(guān)聯(lián)性的二聯(lián)體頻數(shù)以及序列近程關(guān)聯(lián)性的三聯(lián)體、四聯(lián)體及五聯(lián)體頻數(shù)是較為合理的序列特征參數(shù)。

選取供體端除GT外上游30個(gè)位點(diǎn)，內(nèi)含子及受體端除AG外下游3個(gè)位點(diǎn)，定義以下參數(shù)：①統(tǒng)計(jì)分析A、T、G、C的頻數(shù)，定義多樣性源參數(shù)；②統(tǒng)計(jì)分析16種二聯(lián)體的數(shù)目，定義多樣性源參數(shù)；③統(tǒng)計(jì)分析64種三聯(lián)體的數(shù)目，定義多樣性源參數(shù)；④統(tǒng)計(jì)分析256種四聯(lián)體數(shù)目，定義多樣性源參數(shù)；⑤統(tǒng)計(jì)分析1 024種五聯(lián)體數(shù)目，定義多樣性源參數(shù)。

2 結(jié)果與分析

2.1 內(nèi)含子保留的長(zhǎng)度統(tǒng)計(jì)分析結(jié)果

去除AltSplice-rel3.genes.txt文件中的重復(fù)序列，剔除非GT-AG剪切類型的內(nèi)含子，得到4 525條內(nèi)含子保留序列，統(tǒng)計(jì)其序列長(zhǎng)度并繪制內(nèi)含子保留的長(zhǎng)度分布圖（圖2）。

從圖2可以看出，長(zhǎng)度為100～200 bp的保留型內(nèi)含子在所有保留型內(nèi)含子中所占的比例最大，約為26.5%；長(zhǎng)度小于400 bp的大約占70%；長(zhǎng)度超過1 000 bp的部分所占比例不到總數(shù)的10%。這說明剪接效率可能與相應(yīng)的內(nèi)含子長(zhǎng)度有關(guān)；保留型內(nèi)含子最終會(huì)和外顯子一同被保留下來，如果內(nèi)含子過長(zhǎng)就會(huì)形成一個(gè)很長(zhǎng)的外顯子從而使剪接后的外顯子變得不穩(wěn)定。保留型內(nèi)含子的長(zhǎng)度分布特征也與內(nèi)含子漫長(zhǎng)的進(jìn)化過程有關(guān)。

2.2 序列特征分析

2.2.1 組成型內(nèi)含子與保留型內(nèi)含子受體端序列分析 截取全部組成型內(nèi)含子和保留型內(nèi)含子的受體端除AG外上游30 nt，下游10 nt的序列片段，利用在線軟件WebLogo 3[24]進(jìn)行序列保守性分析，結(jié)果如圖3所示。

從圖3可以看出，在內(nèi)含子的3′端上游大約27 nt范圍內(nèi)T和C含量非常豐富，這個(gè)區(qū)域?yàn)槎嗑坂奏ば蛄校≒PT）。比較組成型內(nèi)含子3′端和保留型內(nèi)含子3′端序列保守性后發(fā)現(xiàn)，保留型內(nèi)含子3′端信號(hào)強(qiáng)度明顯弱于組成型內(nèi)含子的信號(hào)強(qiáng)度。組成型內(nèi)含子在其3′端AG上游第3位、第7～12位等七個(gè)位點(diǎn)上胸腺嘧啶出現(xiàn)的頻率高于胞嘧啶，而其他位點(diǎn)則更加偏向于胞嘧啶，在兩者含量上的差別并不明顯。保留型內(nèi)含子的多聚嘧啶序列相對(duì)于胞嘧啶更加偏好使用胸腺嘧啶，且胸腺嘧啶含量明顯高于胞嘧啶。

2.2.2 組成型內(nèi)含子與保留型內(nèi)含子供體端序列分析 截取組成型內(nèi)含子和保留型內(nèi)含子5′端GT上游19 nt（不包括GT）、下游19 nt的序列片段，WebLogo分析結(jié)果如圖4所示。

比較圖4A和圖4B，保留型內(nèi)含子與組成型內(nèi)含子的供體端前兩個(gè)位點(diǎn)和后四個(gè)位點(diǎn)的保守型很強(qiáng)，在5′端上游第二個(gè)位點(diǎn)上胸腺嘧啶和鳥嘌呤雖然含量大致相當(dāng)，但使用偏好明顯不同，組成型內(nèi)含子偏向于使用胸腺嘧啶，而保留的內(nèi)含子傾向選擇鳥嘌呤。5′端下游第一個(gè)位點(diǎn)上兩者的腺嘌呤和鳥嘌呤含量大致相同，而組成型內(nèi)含子偏向于腺嘌呤，保留型內(nèi)含子偏向于鳥嘌呤。5′端下游第三個(gè)位點(diǎn)情況類似，腺嘌呤和胞嘧啶含量大致一樣，對(duì)于組成型內(nèi)含子偏好腺嘌呤，而在保留型內(nèi)含子中偏好使用胞嘧啶。

2.3 預(yù)測(cè)結(jié)果

利用5-fold交叉檢驗(yàn)來評(píng)判預(yù)測(cè)方法的優(yōu)劣，即將數(shù)據(jù)集分為五組，一組作為測(cè)試集，其他四組作為訓(xùn)練集；用訓(xùn)練集中提取出的序列特征來預(yù)測(cè)測(cè)試集中的序列。經(jīng)過5次交叉檢驗(yàn)后取平均值。采用總精度（TA）、敏性指標(biāo)（Sn）、特異性指標(biāo)（Sp）及相關(guān)系數(shù)（MCC）共4個(gè)指標(biāo)評(píng)價(jià)預(yù)測(cè)結(jié)果。在閾值ξ為-0.45時(shí)，預(yù)測(cè)結(jié)果如表2所示。從表2結(jié)果可以看出，利用IDQD算法分類預(yù)測(cè)組成型內(nèi)含子和內(nèi)含子保留的總精度超過了70%，同時(shí)，敏感性指標(biāo)和特異性指標(biāo)也分別達(dá)到了70.23%和70.30%，相關(guān)系數(shù)為0.398 5。

為檢驗(yàn)正負(fù)集樣本序列長(zhǎng)度對(duì)預(yù)測(cè)結(jié)果是否有影響。本研究分析在內(nèi)含子上、下游選取不同序列長(zhǎng)度對(duì)預(yù)測(cè)結(jié)果的影響，結(jié)果見表3。由表3可以看出，在供體端上游堿基數(shù)目選取30 nt，受體端下游堿基數(shù)目選取3 nt時(shí)，預(yù)測(cè)結(jié)果最接近最佳值。

正負(fù)集樣本容量的不同也會(huì)直接影響最終的預(yù)測(cè)結(jié)果。在以上分析的基礎(chǔ)上，分析樣本容量變化對(duì)預(yù)測(cè)結(jié)果的影響，列舉5組比較有代表性的數(shù)據(jù)進(jìn)行說明，具體見表4。從表4中可以看到，隨著負(fù)集樣本容量逐漸增大，敏感性指標(biāo)（Sn）逐漸降低，而特異性指標(biāo)（Sp）基本呈增長(zhǎng)趨勢(shì)，且在負(fù)集樣本容量為正集的1.5倍時(shí)特異性指標(biāo)突增為83.16%，綜合考慮評(píng)價(jià)指標(biāo)，在正負(fù)集樣本數(shù)目比值為1∶1.5時(shí)，預(yù)測(cè)效果更接近最佳值。

3 討論

利用Weblogo 3在線軟件初步分析內(nèi)含子保留型與組成型的供體端和受體端的堿基保守性特征，保留型內(nèi)含子與組成型內(nèi)含子的序列保守性存在明顯差異。在此基礎(chǔ)上，選取單堿基、二聯(lián)體、三聯(lián)體、四聯(lián)體、五聯(lián)體5個(gè)多樣性指標(biāo)對(duì)兩者進(jìn)行區(qū)分，從預(yù)測(cè)結(jié)果來看，敏感性、特異性、總精度3個(gè)指標(biāo)均超過了70%。因此，從序列信息角度來區(qū)分組成型內(nèi)含子和保留型內(nèi)含子是可行的，說明序列信息在內(nèi)含子保留可變剪接過程中有著重要作用。由于內(nèi)含子上下游序列的長(zhǎng)度以及正負(fù)集樣本的容量對(duì)內(nèi)含子保留的預(yù)測(cè)有一定影響，本研究通過調(diào)整內(nèi)含子上下游序列的長(zhǎng)度以及正負(fù)集樣本比例得到在供體端上游堿基數(shù)目選取30 nt，受體端下游堿基數(shù)目選取3 nt時(shí)，負(fù)集樣本容量為正集的1.5倍時(shí)預(yù)測(cè)結(jié)果最接近最佳值。

應(yīng)用貝葉斯二次判別函數(shù)判斷真假時(shí)，其分界值并不一定為零。經(jīng)過多次調(diào)整閾值ξ最終獲得較為滿意的預(yù)測(cè)結(jié)果，即當(dāng)閾值ξ取-0.45時(shí)，敏感性、特異性、總精度3個(gè)指標(biāo)均處于合理范圍，可以認(rèn)為此時(shí)預(yù)測(cè)值是最接近最佳值的。

本研究只選擇了符合GT-AG規(guī)則的剪接位點(diǎn)，沒有考慮其他非標(biāo)準(zhǔn)剪接位點(diǎn)，而且只選取了單堿基、二聯(lián)體、三聯(lián)體、四聯(lián)體、五聯(lián)體頻數(shù)5個(gè)特征參數(shù)，選取的特征參數(shù)相對(duì)較少。研究顯示剪接效率依賴于內(nèi)含子的長(zhǎng)度[25，26]，也可以嘗試將內(nèi)含子長(zhǎng)度作為IDQD算法的一個(gè)參數(shù)；組成型內(nèi)含子受體端上游的多聚嘧啶序列與保留型內(nèi)含子的多聚嘧啶序列特征明顯不同，該段序列的GC含量[27]也可以作為IDQD算法的一個(gè)參數(shù)，也可整合翻譯水平保留型內(nèi)含子的特征，如將兩類內(nèi)含子中終止密碼子出現(xiàn)的頻率作為IDQD算法的輸入?yún)?shù)。如果將這些特征整合到IDQD算法中，預(yù)測(cè)精度可能會(huì)有進(jìn)一步的提高。

致謝：工作中梁棟和魏官云給予很多幫助和很好的建議，在此深表謝意。

參考文獻(xiàn)：

[1] SHARP P A.Split genes and RNA splicing[J].Cell，1994，77（6）：805-815.

[2] WAHL M C，WILL C L，L？譈HRMANN R. The spliceosome： Design principles of a dynamic RNP machine[J].Cell，2009， 136（4）：701-718.

[3] LUCO R F，ALLO M，SCHOR I E，et al. Epigenetics in alternative pre-mRNA splicing[J].Cell，2011，144（1）：16-26.

[4] 王科俊，呂俊杰，馮偉興，等.可變剪接與疾病的生物信息學(xué)研究概況[J].生命科學(xué)研究，2011，15（1）：86-94.

[5] 李稚鋒，王正志，張成崗.真核基因可變剪接研究現(xiàn)狀與展望[J]. 生物信息學(xué)，2004，2（5）：35-39.

[6] 章天驕.可變剪接的生物信息數(shù)據(jù)分析綜述[J].生物信息學(xué)，2012，10（1）：61-64.

[7] 林魯萍，馬 飛，王義權(quán).基因選擇性剪接的生物信息學(xué)研究概況[J].遺傳，2005，27（6）：1001-1006.

[8] 蔡鈺深.可變剪接研究的主流方法[D].臺(tái)北：臺(tái)灣大學(xué)，2012.

[9] WANG L，XI Y，YU J，et al. A statistical method for the detection of alternative splicing using RNA-seq[J]. PloS One， 2010，5（1）：e8529.

[10] 呂俊杰.采用智能方法的可變剪接調(diào)控機(jī)制與相關(guān)疾病研究[D].哈爾濱：哈爾濱工程大學(xué)，2012.

[11] 何 濤，王端青，胡亞歐，等.基于RNA-Seq數(shù)據(jù)識(shí)別果蠅剪接位點(diǎn)和可變剪接事件[J].中國(guó)科學(xué)：生命科學(xué)，2011，41（10）： 1016-1023.

[12] 李燕青.基于支持向量機(jī)方法的剪接位點(diǎn)預(yù)測(cè)[D].福州：福建農(nóng)林大學(xué)，2012.

[13] CUI Y， HAN J Q， ZHONG D X，et al.A novel computational method for the identification of plant alternative splice sites[J].Biochem Biophys Res Commun，2013， 431（2）：221-224.

[14] 朱紅梅，王家廒，趙燕南，等.延時(shí)HMM在基因剪接供體位點(diǎn)識(shí)別中的應(yīng)用[J].計(jì)算機(jī)工程，2007，33（5）：1-3.

[15] 程國(guó)建，趙 斐，吳曉怡.神經(jīng)網(wǎng)絡(luò)在基因序列預(yù)測(cè)中的應(yīng)用研究[J].微計(jì)算機(jī)信息，2008，24（11）：264-2645.

[16] 閆曉強(qiáng).RNA剪接識(shí)別研究[J].長(zhǎng)春：吉林大學(xué)，2009.

[17] 周艷紅，王 卉，楊 雷.基于特征挖掘與融合的剪接位點(diǎn)識(shí)別[J].華中科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2006，34（12）：117-120.

[18] 蔡 祿.表觀遺傳學(xué)前沿[M].北京：清華大學(xué)出版社，2012.

[19] 金 鷹.選擇性剪接的理論預(yù)測(cè)[J].激光生物學(xué)報(bào)，2008，17（2）：283-285.

[20] 張利絨，羅遼復(fù).多樣性指標(biāo)用于基因中剪切位點(diǎn)的識(shí)別[J].生物化學(xué)與生物物理進(jìn)展，2004，31（1）：77-82.

[21] 邢永強(qiáng)，張利絨，羅遼復(fù).人類基因組盒式外顯子和內(nèi)含子保留的可變剪接位點(diǎn)預(yù)測(cè)[J].生物物理學(xué)報(bào)，2008，24（5）：393-400.

[22] 張利絨，羅遼復(fù)，邢永強(qiáng)，等.人類基因組中可變和組成性剪接位點(diǎn)的預(yù)測(cè)[J].生物化學(xué)與生物物理進(jìn)展，2008，35（10）：1188-1194.

[23] LERIVRAY H， MEREAU A， OSBORNE H B.Our favourite alternative splice site[J].Biol Cell， 2006， 5（98）： 317-321.

[24] CROOKS G E， HON G， CHANDONIA J M， et al. WebLogo： A sequence logo generator[J]. Genome research， 2004， 14（6）：1188-1190.

[25] 馬 猛，汪 洋.利用計(jì)算方法識(shí)別定義內(nèi)含子保留的基因組特征[J].計(jì)算機(jī)工程與應(yīng)用，2012，48（28）：36-41.

[26] WANG W H， MAUCUER A， GUPTA A， et al.Structure of phosphorylated SF1 bound to U2AF65 in an essential splicing factor complex[J]. Structure， 2013， 2（21）：197-208.

[27] PASTUSZAK A W， JOACHIMIAK M P， BLANCHETTE M， et al.An SF1 affinity model to identify branch point sequences in human introns[J]. Nucleic Acids Res，2011，39（6）： 2344-2356.