打开文本图片集
摘要：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）表示寡核苷酸成为5’SS位点的先验概率；P（-ss）表示寡核苷酸成为类供体信号位点的先验概率；P（oligo|ss）表示寡核苷酸出现在5’SS位点的可能性；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.