基于序列信息預測內含子保留型可變剪接

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

關鍵詞：可變剪接；內含子保留；多樣性增量；二次判別函數

中圖分類號：Q61 文獻標識碼：A 文章編號：0439-8114（2016）01-0208-05

DOI：10.14088/j.cnki.issn0439-8114.2016.01.055

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

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

1 材料和方法

1.1 可變剪接數據集的構建

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

1.2 選取特征參數

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

選取供體端除GT外上游30個位點，內含子及受體端除AG外下游3個位點，定義以下參數：①統計分析A、T、G、C的頻數，定義多樣性源參數；②統計分析16種二聯體的數目，定義多樣性源參數；③統計分析64種三聯體的數目，定義多樣性源參數；④統計分析256種四聯體數目，定義多樣性源參數；⑤統計分析1 024種五聯體數目，定義多樣性源參數。

2 結果與分析

2.1 內含子保留的長度統計分析結果

去除AltSplice-rel3.genes.txt文件中的重復序列，剔除非GT-AG剪切類型的內含子，得到4 525條內含子保留序列，統計其序列長度并繪制內含子保留的長度分布圖（圖2）。

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

2.2 序列特征分析

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

從圖3可以看出，在內含子的3′端上游大約27 nt范圍內T和C含量非常豐富，這個區域為多聚嘧啶序列（PPT）。比較組成型內含子3′端和保留型內含子3′端序列保守性后發現，保留型內含子3′端信號強度明顯弱于組成型內含子的信號強度。組成型內含子在其3′端AG上游第3位、第7～12位等七個位點上胸腺嘧啶出現的頻率高于胞嘧啶，而其他位點則更加偏向于胞嘧啶，在兩者含量上的差別并不明顯。保留型內含子的多聚嘧啶序列相對于胞嘧啶更加偏好使用胸腺嘧啶，且胸腺嘧啶含量明顯高于胞嘧啶。

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

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

2.3 預測結果

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

為檢驗正負集樣本序列長度對預測結果是否有影響。本研究分析在內含子上、下游選取不同序列長度對預測結果的影響，結果見表3。由表3可以看出，在供體端上游堿基數目選取30 nt，受體端下游堿基數目選取3 nt時，預測結果最接近最佳值。

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

3 討論

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

應用貝葉斯二次判別函數判斷真假時，其分界值并不一定為零。經過多次調整閾值ξ最終獲得較為滿意的預測結果，即當閾值ξ取-0.45時，敏感性、特異性、總精度3個指標均處于合理范圍，可以認為此時預測值是最接近最佳值的。

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

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

參考文獻：

[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] 王科俊，呂俊杰，馮偉興，等.可變剪接與疾病的生物信息學研究概況[J].生命科學研究，2011，15（1）：86-94.

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

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

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

[8] 蔡鈺深.可變剪接研究的主流方法[D].臺北：臺灣大學，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] 呂俊杰.采用智能方法的可變剪接調控機制與相關疾病研究[D].哈爾濱：哈爾濱工程大學，2012.

[11] 何 濤，王端青，胡亞歐，等.基于RNA-Seq數據識別果蠅剪接位點和可變剪接事件[J].中國科學：生命科學，2011，41（10）： 1016-1023.

[12] 李燕青.基于支持向量機方法的剪接位點預測[D].福州：福建農林大學，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] 朱紅梅，王家廒，趙燕南，等.延時HMM在基因剪接供體位點識別中的應用[J].計算機工程，2007，33（5）：1-3.

[15] 程國建，趙 斐，吳曉怡.神經網絡在基因序列預測中的應用研究[J].微計算機信息，2008，24（11）：264-2645.

[16] 閆曉強.RNA剪接識別研究[J].長春：吉林大學，2009.

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

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

[19] 金 鷹.選擇性剪接的理論預測[J].激光生物學報，2008，17（2）：283-285.

[20] 張利絨，羅遼復.多樣性指標用于基因中剪切位點的識別[J].生物化學與生物物理進展，2004，31（1）：77-82.

[21] 邢永強，張利絨，羅遼復.人類基因組盒式外顯子和內含子保留的可變剪接位點預測[J].生物物理學報，2008，24（5）：393-400.

[22] 張利絨，羅遼復，邢永強，等.人類基因組中可變和組成性剪接位點的預測[J].生物化學與生物物理進展，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].計算機工程與應用，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.