可變剪接的理論工作研究進展

梁棟

摘要：pre-mRNA的可變剪接是調控真核生物基因表達以及蛋白質多樣性的重要機制，使得可變剪接與組織分化或疾病發生有密切聯系。從序列水平和表達水平對可變剪接的理論工作研究進展做了介紹，主要由相關數據庫和常見的生物信息學方法兩部分組成。

關鍵詞：可變剪接；理論預測；研究進展；數據庫；生物信息學方法

中圖分類號：Q343.1 文獻標識碼：A 文章編號：0439-8114（2016）04-0820-05

DOI：10.14088/j.cnki.issn0439-8114.2016.04.002

Research Progress of Theoretical Work about Alternative Splicing

LIANG Donga，b

（a.Institute of Bioengineering and Technology，b.School of Mathematics Physics and Biological Engineering/Inner Mongolia Key Laboratory of Biomass-Energy Conversion，The Institute of Bioengineering and Technology， Baotou 014010， Inner Mongolia， China）

Abstract： The alternative splicing of pre-mRNA plays an important role in regulation of eukaryotic gene expression and diversity of protein， which made a close relation between alternative splicing and tissue differentiation or occurrence of disease. Then， this paper gives the introduction about research progress of theoretical work of alternative splicing from sequence level and expression level， which mainly is divided into relative database and common Bioinformatics Method.

Key words： alternative splicing； theoretical prediction； research progress； database； bioinformatics method

可變剪接是mRNA前體經過不同的剪接方式，形成不同的剪接異構體，最終產生不同成熟mRNA的過程。在基因可變剪接中，經過調控元件和剪接因子的相互作用，相關外顯子跳躍或者保留，形成不同的組合，進而產生不同的成熟mRNA和蛋白產物[1]。所以，可變剪接可以調控基因功能，影響組織表型而產生疾病[2]，例如神經系統的發育和程序性細胞死亡[3，4]。研究表明，大約90%的人類基因會發生可變剪接，極大地增加了蛋白質的多樣性[5，6]。

由于試驗研究周期較長，花費較大，所以通過生物信息學方法可變剪接具有獨特優勢。從pre-mRNA序列水平分析，主要涉及剪接信號和剪接因子結合位點的預測。從基因表達水平分析，主要涉及外顯子表達水平變化。表達數據主要有傳統的表達序列標簽（EST）、基因芯片和RNA-seq。EST是從一個隨機選擇的cDNA克隆進行5′端和3′端單一次測序獲得的短的cDNA部分序列，代表一個完整基因的一小部分。基因芯片是通過與一組已知序列的核酸探針雜交進行核酸序列測定的方法，即在一塊基片表面固定了序列已知的靶核苷酸的探針，與待測核酸序列互補匹配。RNA-Seq利用高通量測序技術對組織或細胞中所有RNA反轉錄而成的cDNA文庫進行測序，通過統計相關（Reads）數計算出不同RNA的表達量，發現新的轉錄本[7]；目前表達分析的趨勢是RNA-Seq技術，主要原因是RNA-Seq作為一個“開放系統”，可分別檢測有參考序列和無參考序列的轉錄組，進一步可發現和尋找新的信息[8]。

1 相關數據庫

通過文獻調研，收集了在序列水平上進行可變剪接理論工作所需的相關數據庫，見表1。

在表達水平上進行可變剪接理論工作主要用到以下數據庫：GEO、ArrayExpress和SRA（表2）。

GEO數據庫（Gene Expression Omnibus，http：//www.ncbi.nlm.nih.gov/geo）是NCBI在2000年開發的一個基因表達和雜交微陣列數據庫，包括3個相關數據庫：平臺（platform）、樣本（sample）和系列（series）[9]。ArrayExpress（http：//www.ebi.ac.uk/arrayexpress）是微陣列基因表達數據的公共知識數據庫，可以保存微陣列平臺經過注釋的數據。ArrayExpress下設兩個子數據庫：試驗數據集和基因表達圖譜[10]。SRA數據庫（Sequence Read Archive，http：//www.ncbi.nlm.nih.gov/sra）儲存來自454、IonTorrent、Illumina、SOLiD和Helicos等測序技術測得的高通量短片段測序，可以在該數據庫下載RNA-seq數據[11]。RNA-Seq數據的下載方法：登陸SRA數據庫網址，在檢索欄輸入查詢的W。

2 常見生物信息學方法

2.1 位置特異性得分矩陣

位置特異性得分矩陣（PSSM）主要對短序列內各堿基的偏好性進行量化，由于剪接因子的結合位點大多有保守的motif，所以PSSM在可變剪接中的應用是根據剪接因子結合序列內部的保守位點，對其進行判斷。PSSM是一個4×n的矩陣，行數為4代表四種常見堿基A、T、C、G；列數為n代表剪接因子結合位點長度為n。設堿基？琢∈{A，T，C，G}，則PSSM定義[12]如下：

P（？琢，i）=■ （1）

Sα=■ （2）

w（？琢，i）=log10■ （3）

score=■w（？琢，i） （4）

式中，Ni是已知相關剪接因子結合位點序列的總數。f（？琢，i）是在所有Ni條序列中，在位置i處的堿基 α出現頻率。P（？琢，i）表示在位置i處的出現頻率。P？琢是堿基α在背景序列中出現的頻率。w（？琢，i）是矩陣元素值，即堿基α在位置 i上的頻率信息，反映堿基α的位置偏好性。Sα是考慮由于觀察失誤導致矩陣元素為0的情況，而引入的數值極小的偽數目。公式（1）至公式（3）通過建立PSSM模型，利用堿基在指定位置的出現頻率，通過對數表示該堿基在指定位置的偏好性，最后利用公式（4）完成打分過程，即得到每條位點結合序列的各堿基偏好性分布。最后設定閾值T，來判斷各條序列是否成為結合位點，若序列打分S≥T，則認為該序列可能是剪接因子的結合位點。

2.2 貝葉斯網絡

貝葉斯網絡（Bayesian network，BN） 是一系列變量的聯合概率分布的圖形表示，是一種系統描述變量之間關系的工具。BN由兩部分組成：①貝葉斯網絡結構圖，是一個有向無環圖（DAG）。圖中節點代表相應變量，節點之間聯系代表條件獨立關系。②節點之間的條件概率表（CPT）。BN理論基礎由公式（5）至公式（7）組成：

P（X1，X2，…，Xn）=P（X1）P（X2|X1）…P（Xn|X1X2，…，Xn）

（5）

表示貝葉斯網絡之間的鏈規則（Chain rule）。貝葉斯定理見公式（6）：

P（A|B）=■（6）

公式（7）表示貝葉斯網路中各變量間的條件獨立性：

P（X1，X2，…，Xn）=■P（X1|（π（Xi）） （7）

李驁等[13]通過網絡建模，預測真核生物DNA序列的剪接位點。首先其利用受體和供體上下游位點構建受體和供體貝葉斯網絡模型，再利用貝葉斯網絡學習的最大似然法求解各節點在父節點條件下的條件概率分布，即公式（8）中得參數？茲j：

L（？茲）=■P（Yj|YPA（i），？茲j） （8）

式中，Yj表示受體和供體模型中各節點；YPA（i）表示Yj的父節點；N為貝葉斯網絡節點數目，？茲={？茲1，…，？茲j）；對公式（8）進行log變換求解，見公式（9）：

logL（？茲）=■logP（Yj|YPA（i），？茲j） （9）

最后通過似然估計函數求得參數？茲′=argmax（logL（？茲））。

Alexander Churbanov等[14]設計了7-mer寡核苷酸傳感器（47共16 384中寡核苷酸組合）來計算splice信號和splice-like信號處的寡核苷酸概率。傳感器拓撲分布見圖1。

人們使用GIGOgene[15]選取179 079個供體信號和34 258 282個類供體信號 ，以及179 076個受體信號44 353 884個類受體信號。在位點附近給定寡核苷酸條件下，可以通過公式（10）計算該位點成為剪接信號的概率：

P（ss|oligo）=■ （10）

式中，P（ss）表示寡核苷酸成為5SS位點的先驗概率；P（-ss）表示寡核苷酸成為類供體信號位點的先驗概率；P（oligo|ss）表示寡核苷酸出現在5SS位點的可能性；P（oligo|-ss）表示寡核苷酸出現在類供體信號位點的可能性。最后進行計算結果對數化，圖1中的（a）對數化block0的計算結果；針對（b）則求取所有block的計算結果對數化之和。

2.3 支持向量機

支持向量機（Support vector machine，SVM）是一種從少量樣本中提取分類信息而實現對大量樣本進行分類的機器學習方法[16]，適用于解決小樣本、非線性及高維問題，其基本思想是將輸入空間的樣本通過某種非線性函數關系映射到一個特征空間中，使兩類樣本在此特征空間中線性可分，最后在此特征空間中求取最優先線性分類面。最優線性分類面求取的原理是分類線不僅將兩類樣本無差分開，且兩類之間的距離最大[17]。

設線性分類樣本集（xi，yi），i=1，2，…，x∈Rd，y {-1，+1}，d代表空間維度，則線性分類判別函數為公式（11）：

g（x）=w·x+b （11）

式中，X為輸入向量，w為權值向量，b為分類閾值；則：

w·x+b>0 yi=+1；

w·x+b<0 yi=-1 （12）

要求分類線對所有樣本正確分類，則有公式（13）[由公式（12）轉化而來]：

yi（w·xi+b）-1≥0 （13）

分類面方程為公式（14）：

w·x+b=0 （14）

由解析幾何知識可知（r為各分類樣本與分類面的距離）：

r=■ （15）

將判別函數歸一化，使兩類所有樣本均滿足|g（x）≥1|，指離分類面最近樣本|g（x）|=1。則兩類樣本分類間隔為■，因此，是分類間隔最大，只需使||w||最小。其中離分類面最近點且平行于最優分類面的超平面分類樣本稱作支持向量。利用拉格朗日函數將最優分類面求解問題轉化為對偶問題，見公式（16）：

f（x）=sgn{∑？琢iyi（xi·x）+b} （16）

式中，sgn（）為符號函數，b是分類閾值，給定待判別樣本x，利用公式（16）計算，可確定該樣本種類。

對于非線性可分樣本，須引入松弛變量和懲罰因子來允許一定錯誤。用內積K（xi·x）代替最優分類面中的點積yi（xi·x），最優分類函數變為公式（17）：

f（x）=sgn{∑？琢iKi（xi·x）+b} （17）

常見的內積函數如下：

線性內積函數K（xi·xj）=xixj （18）

多項式內積函數K（xi·xj）=exp（？姿xixj+r），？姿>0（19）

核函數型內積函數K（xi·xj）=exp（-？姿||xi-xj||2）（20）

Dror等[18]選取228個序列特征（包括外顯子長度，內含子序列上下游保守性等）作為輸入向量xi=（x■■，…，x■■），yi=-1或yi=+1表示是組成型或者可變外顯子。最后使用軟件SVM light （http：//svmlight.joachims.org）來完成運算過程。

2.4 剪接指數

剪接指數算法最早由Clark等[19]提出，目前也可用于分析Affymetrix公司開發的外顯子芯片。由于外顯子的表達會受到所在基因的差異表達影響，而剪接指數算法通過將外顯子表達信號與基因表達信號做比值[20][見公式（21）]，實現對外顯子表達信號的標準化，若比值較大，則外顯子剪接趨勢越大。

標準化相對信號=■ （21）

標準化之后，將外顯子標準表達信號樣本分為處理組（Treatment）和對照組（Control），剪接指數的定義就是二者比值的對數，見公式（22）：

剪接指數=log2■ （22）

剪接指數的意義是反映外顯子差異表達的定量，最后利用t檢驗對剪接指數進行假設檢驗，若統計顯著性P值越小，外顯子差異表達越大，可變剪接發生的可能性越大。

2.5 ARH-seq

ARH（AS robust prediction method based on entropy）是依賴于熵值預測可變剪接的一種方法。 Rasche等[21]利用ARH的方法和RNA-Seq數據來進行預測。獲取從不同人類組織中得到的RNA-Seq數據，用bowtie將讀段（Reads）與基因組比對（Map），根據基因組結構得到外顯子連接處（Exon junctions）。然后計算外顯子及其連接處（Exon combi count values）RPKM值，見公式（23）：

RPKM=■（23）

式中，exon length即比對到基因組的長度；mapped reads即所有長度比對到基因組的reads總數；total exon reads即每個長度比對到基因組的reads數目。RPKM值的意義是對外顯子表達進行衡量。之后ARH的計算流程如下：

令基因g有m個外顯子，在兩種生物條件c和t下，則？準g，e，t，e=1，…，m，代表基因g的第e個外顯子在條件t下得RPKM值。ζg，e是剪接偏差，表示各外顯子表達倍數差異與所有外顯子表達倍數差異中位數的偏差，計算公式見公式（24）（加1是避免分母為0）：

ζg，e=log2（■）-■（log2（■））（24）

Pg，e的意義是各外顯子的剪接概率，其中∑ePg，e=1（表示可變剪接一定發生），見公式（25）：

Pg，e=■ （25）

為確定各外顯子剪接概率是否均勻分布，為各外顯子定義熵值（參考申農熵定義），見公式（26）：

Hg（Pg，1，……，Pg，m）=-■Pg，e*log2（Pg，e） （26）

熵值Hg（Pg，1，……，Pg，m）受到外顯子數量影響，利用Hg（Pg，1，……，Pg，m）理論最大值與其做差來解決這個問題，見公式（27）：

max（Hg）-Hg=log2（m）-Hg（Pg，1，……，Pg，m） （27）

除熵值之外，還需利用分位數之比來描述可變剪接發生概率。將基因g各外顯子表達差異倍數按遞增順序排列，取第一個和第三個四分位數，即Q.25，g，e=1，……，m（■）和Q.75，g，e=1，……，m（■），將二者做比值，得■，若可變剪接發生概率較低，則比值趨近于1。最后，結合分位數比值和公式（11），得ARH剪接預測結果，見公式（28）：

ARHg=■×（max（Hg）-Hg） （28）

若ARH越大，則可變剪接發生的可能性越大。

3 小結與展望

本文從兩個層次對可變剪接的一些理論工作或者對可變剪接理論研究有參考價值的工作進行了歸納，即：①序列水平，其一指預測序列上的剪接位點和剪接調控元件，例如貝葉斯網絡和PSSM；其二指利用序列特征對外顯子的類型進行預測，例如支持向量機（SVM）；②表達水平，主要指通過芯片數據和RNA-seq數據（尤其涉及外顯子的數據）的表達變化，對剪接事件發生的可能性進行預測，例如剪接指數和ARH-seq。

本文認為今后可變剪接的研究主要集中在以下幾個方面：

1）構建可變剪接調控網絡。可變剪接網絡涉及多個因子和調控元件之間的相互作用，目前的工作主要集中于預測蛋白質與核酸的結合，例如馬昕等[22]通過提出特征PSSM-PP（包含影響蛋白質與RNA結合的氨基酸理化特征），選用隨機森林樹算法來預測蛋白質中的RNA結合殘基。Ahmad等[23]利用從蛋白質骨干原子的坐標計算得到的電荷、偶極距、四極距等特征，采用神經網絡（NN）方法構建RNA結合蛋白的預測模型。Zhao等[24]通過將體積分數校正引入DFIRE統計能量函數中，來預測蛋白質與RNA的結合親和性。但是，關于蛋白質與核酸結合后的綜合調控作用的研究還較為缺乏。

2）深入探究可變剪接的致病機制。目前很多工作通過研究基因芯片和RNA-seq數據中蘊含的基因表達信息，來闡述基因與疾病的相關性。但可變剪接對基因表達的影響還不是很清楚，這其中可能涉及可變剪接產物的形成，以及信號通路調節基因表達等過程。

綜上所述，可變剪接對真核生物表達調控和蛋白質多樣性起到重要作用。可變剪接的理論預測工作仍有很大上升空間，通過優化現有的理論模型，獲取更優質精確的數據，可變剪接的本質最終將被人們揭曉，轉化為實際應用而造福人類。

參考文獻：

[1] MITTENDORF K F， DEATHERAGE C L， OHI R， et al. Tailoring of membrane proteins by alternative splicing of pre-mRNA[J]. Biochemistry，2012，51（28）：5541-5556.

[2] LU Z X，JIANG P，XING Y. Genetic variation of pre-mRNA alternative splicing in human populations[J]. Wiley Interdisciplinary Reviews： RNA，2012，3（4）：581-592.

[3] LI Q， LEE J， BLACK D L. Neuronal regulation of alternative pre-mRNA splicing[J]. Nature Reviews Neuroscience，2007，8（11）：819-831.

[4] JIANG Z H， WU J Y. Alternative splicing and programmed cell death[J]. Proceedings of the Society for Experimental Biology and Medicine，1999，220（2）：64-72.

[5] PAN Q， SHAI O， LEE L J， et al. Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing[J]. Nat Genet，2008，40（12）：1413-1415.

[6] WANG E T， SANDBERG R， LUO S， et al. Alternative isoform regulation in human tissue transcriptomes[J]. Nature， 2008，456（7221）：470-476.

[7] 祁云霞，劉永斌，榮威恒.轉錄組研究新技術：RNA-Seq及其應用[J].遺傳，2011，33（11）：1191-1202.

[8] 高 山，歐劍虹，肖 凱.R語言與Bioconductor生物信息學應用[M].天津：天津科技翻譯出版有限公司，2014.168-169.

[9] 羅烈偉.基因表達與雜交陣列數據庫：GEO[EB/OL].中國科技論文在線，2008.1-9.

[10] 周大瓊，曹繼華，任力鋒，等.基因芯片數據庫GEO與ArrayExpress的使用及比較分析[J].中國現代醫學雜志，2014，24（12）：38-42.

[11] 熊筱晶.NCBI高通量測序數據庫SRA介紹[J].生命的化學， 2010，30（6）：959-963.

[12] 楊科利.基于序列信息的轉錄因子結合位點和啟動子理論預測[D].呼和浩特：內蒙古大學，2007.

[13] 李 驁，王 濤，馮煥清，等.基于貝葉斯網絡的DNA序列剪接位點預測[J].生物物理學報，2003，19（4）：431-436.

[14] CHURBANOV A，ROGOZIN I B，DEOGUN J S，et al.Method of predicting splice sites based on signal interactions[J].Biol Direct，2006，DOI：10.1186/1745-6150-1-10.

[15] CHURBANOV A，PAULEY M，QUEST D，et al.A method of precise mRNA/DNA homology-based gene structure prediction[J].BMC Bioinformatics，2005，DOI：10.1186/1471-2105-6-261.

[16] 孫 嘯，陸祖宏，謝建明.生物信息學基礎[M].北京：清華大學出版社，2005.

[17] 張 鴿，陳書開.基于SVM的手寫體阿拉伯數字識別[J].軍民兩用技術與產品，2005（9）：41-43.

[18] DROR G，SOREK R，SHAMIR R.Accurate identification of alternatively spliced exons using support vector machine[J]. Bioinformatics，2005，21（7）：897-901.

[19] CLARK T A，SUGNET C W，ARES M J R. Genomewide analysis of mRNA processing in yeast using splicing-specific microarrays[J]. Science，2002，296（5569）：907-910.

[20] 杭興宜，鄧明華，孫志賢，等.外顯子芯片的數據分析和可變剪接預測的新算法[J].軍事醫學，2009，33（5）：465-468.

[21] RASCHE A，LIENHARD M，YASPO M L，et al.ARH-seq： Identification of differential splicing in RNA-seq data[J].Nucleic Acids Res，2014，42（14）：e110.

[22] 馬 昕，郭 靜，孫 嘯.蛋白質中RNA-結合殘基預測的隨機森林模型[J].東南大學學報（自然科學版），2012，42（1）：50-54.

[23] AHMAD S，SARAI A. Analysis of electric moments of RNA-binding proteins：implications for mechanism and prediction[J]. BMC Structural Biology，2011，DOI：10.1186/1472-6807-11-8.

[24] ZHAO H，YANG Y，ZHOU Y.Structure-based prediction of DNA-binding proteins by structural alignment and a volume-fraction corrected DFIRE-based energy function[J].Bioinformatics，2010，26（15）：1857-1863.