可變剪接的理論工作研究進(jìn)展

梁棟

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

關(guān)鍵詞：可變剪接；理論預(yù)測；研究進(jìn)展；數(shù)據(jù)庫；生物信息學(xué)方法

中圖分類號：Q343.1 文獻(xiàn)標(biāo)識碼：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前體經(jīng)過不同的剪接方式，形成不同的剪接異構(gòu)體，最終產(chǎn)生不同成熟mRNA的過程。在基因可變剪接中，經(jīng)過調(diào)控元件和剪接因子的相互作用，相關(guān)外顯子跳躍或者保留，形成不同的組合，進(jìn)而產(chǎn)生不同的成熟mRNA和蛋白產(chǎn)物[1]。所以，可變剪接可以調(diào)控基因功能，影響組織表型而產(chǎn)生疾病[2]，例如神經(jīng)系統(tǒng)的發(fā)育和程序性細(xì)胞死亡[3，4]。研究表明，大約90%的人類基因會發(fā)生可變剪接，極大地增加了蛋白質(zhì)的多樣性[5，6]。

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

1 相關(guān)數(shù)據(jù)庫

通過文獻(xiàn)調(diào)研，收集了在序列水平上進(jìn)行可變剪接理論工作所需的相關(guān)數(shù)據(jù)庫，見表1。

在表達(dá)水平上進(jìn)行可變剪接理論工作主要用到以下數(shù)據(jù)庫：GEO、ArrayExpress和SRA（表2）。

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

2 常見生物信息學(xué)方法

2.1 位置特異性得分矩陣

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

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

Sα=■ （2）

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

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

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

2.2 貝葉斯網(wǎng)絡(luò)

貝葉斯網(wǎng)絡(luò)（Bayesian network，BN） 是一系列變量的聯(lián)合概率分布的圖形表示，是一種系統(tǒng)描述變量之間關(guān)系的工具。BN由兩部分組成：①貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)圖，是一個有向無環(huán)圖（DAG）。圖中節(jié)點代表相應(yīng)變量，節(jié)點之間聯(lián)系代表條件獨立關(guān)系。②節(jié)點之間的條件概率表（CPT）。BN理論基礎(chǔ)由公式（5）至公式（7）組成：

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

（5）

表示貝葉斯網(wǎng)絡(luò)之間的鏈規(guī)則（Chain rule）。貝葉斯定理見公式（6）：

P（A|B）=■（6）

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

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

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

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

式中，Yj表示受體和供體模型中各節(jié)點；YPA（i）表示Yj的父節(jié)點；N為貝葉斯網(wǎng)絡(luò)節(jié)點數(shù)目，？茲={？茲1，…，？茲j）；對公式（8）進(jìn)行l(wèi)og變換求解，見公式（9）：

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

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

Alexander Churbanov等[14]設(shè)計了7-mer寡核苷酸傳感器（47共16 384中寡核苷酸組合）來計算splice信號和splice-like信號處的寡核苷酸概率。傳感器拓?fù)浞植家妶D1。

人們使用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）表示寡核苷酸出現(xiàn)在5SS位點的可能性；P（oligo|-ss）表示寡核苷酸出現(xiàn)在類供體信號位點的可能性。最后進(jìn)行計算結(jié)果對數(shù)化，圖1中的（a）對數(shù)化block0的計算結(jié)果；針對（b）則求取所有block的計算結(jié)果對數(shù)化之和。

2.3 支持向量機(jī)

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

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

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

式中，X為輸入向量，w為權(quán)值向量，b為分類閾值；則：

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

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

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

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

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

w·x+b=0 （14）

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

r=■ （15）

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

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

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

對于非線性可分樣本，須引入松弛變量和懲罰因子來允許一定錯誤。用內(nèi)積K（xi·x）代替最優(yōu)分類面中的點積yi（xi·x），最優(yōu)分類函數(shù)變?yōu)楣剑?7）：

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

常見的內(nèi)積函數(shù)如下：

線性內(nèi)積函數(shù)K（xi·xj）=xixj （18）

多項式內(nèi)積函數(shù)K（xi·xj）=exp（？姿xixj+r），？姿>0（19）

核函數(shù)型內(nèi)積函數(shù)K（xi·xj）=exp（-？姿||xi-xj||2）（20）

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

2.4 剪接指數(shù)

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

標(biāo)準(zhǔn)化相對信號=■ （21）

標(biāo)準(zhǔn)化之后，將外顯子標(biāo)準(zhǔn)表達(dá)信號樣本分為處理組（Treatment）和對照組（Control），剪接指數(shù)的定義就是二者比值的對數(shù)，見公式（22）：

剪接指數(shù)=log2■ （22）

剪接指數(shù)的意義是反映外顯子差異表達(dá)的定量，最后利用t檢驗對剪接指數(shù)進(jìn)行假設(shè)檢驗，若統(tǒng)計顯著性P值越小，外顯子差異表達(dá)越大，可變剪接發(fā)生的可能性越大。

2.5 ARH-seq

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

RPKM=■（23）

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

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

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

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

Pg，e=■ （25）

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

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

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

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

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

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

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

3 小結(jié)與展望

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

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

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

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

綜上所述，可變剪接對真核生物表達(dá)調(diào)控和蛋白質(zhì)多樣性起到重要作用?？勺兗艚拥睦碚擃A(yù)測工作仍有很大上升空間，通過優(yōu)化現(xiàn)有的理論模型，獲取更優(yōu)質(zhì)精確的數(shù)據(jù)，可變剪接的本質(zhì)最終將被人們揭曉，轉(zhuǎn)化為實際應(yīng)用而造福人類。

參考文獻(xiàn)：

[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] 祁云霞，劉永斌，榮威恒.轉(zhuǎn)錄組研究新技術(shù)：RNA-Seq及其應(yīng)用[J].遺傳，2011，33（11）：1191-1202.

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

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

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

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

[12] 楊科利.基于序列信息的轉(zhuǎn)錄因子結(jié)合位點和啟動子理論預(yù)測[D].呼和浩特：內(nèi)蒙古大學(xué)，2007.

[13] 李 驁，王 濤，馮煥清，等.基于貝葉斯網(wǎng)絡(luò)的DNA序列剪接位點預(yù)測[J].生物物理學(xué)報，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] 孫 嘯，陸祖宏，謝建明.生物信息學(xué)基礎(chǔ)[M].北京：清華大學(xué)出版社，2005.

[17] 張 鴿，陳書開.基于SVM的手寫體阿拉伯?dāng)?shù)字識別[J].軍民兩用技術(shù)與產(chǎn)品，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] 杭興宜，鄧明華，孫志賢，等.外顯子芯片的數(shù)據(jù)分析和可變剪接預(yù)測的新算法[J].軍事醫(yī)學(xué)，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] 馬 昕，郭 靜，孫 嘯.蛋白質(zhì)中RNA-結(jié)合殘基預(yù)測的隨機(jī)森林模型[J].東南大學(xué)學(xué)報（自然科學(xué)版），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.