劉國慶

(內(nèi)蒙古科技大學 生命科學與技術(shù)學院，內(nèi)蒙古 包頭 014010)

核小體是真核生物染色質(zhì)的基本結(jié)構(gòu)單元，是DNA雙螺旋纏繞在組蛋白八聚體上形成的復合物。標準的組蛋白八聚體由進化上高度保守的H2A、H2B、H3和H4各兩個拷貝組成[1]。并非所有的核小體都由標準的組蛋白組裝，全基因組范圍內(nèi)還富含一些組蛋白變體，如H2A.Z，H3.3等。組蛋白變體與標準組蛋白間存在一定的序列差異，對染色質(zhì)結(jié)構(gòu)和基因轉(zhuǎn)錄有不同的調(diào)控作用[2-4]。標準組蛋白和非標準組蛋白均由基因組上的若干個基因表達[3-4]。核小體核心DNA的長度約為147 bp[5-9]，而相鄰的核小體之間的鏈接DNA的長度并不恒定(約20-80 bp)。核小體核心顆粒在組蛋白H1的作用下形成30-nm結(jié)構(gòu)，進一步組裝成更高級結(jié)構(gòu)[10]，使基因組DNA包裝到狹小的細胞核中。核小體占據(jù)真核基因組的絕大部分(約75%～90%)[11-12]。

核小體區(qū)域的DNA纏繞于組蛋白上，相比鏈接DNA不易于相關(guān)蛋白因子與之相接觸并結(jié)合，從而導致核小體參與基因轉(zhuǎn)錄、DNA復制、修復、重組以及RNA剪接等眾多生物學過程[11-15]。核小體在DNA序列以外的因素(如重塑蛋白、組蛋白修飾酶、細胞微環(huán)境變化等)或內(nèi)在信號(DNA序列突變)的擾動下其位置時有發(fā)生變化并參與上述生物學過程[11-22]。體外的核小體定位只決定于DNA序列和相鄰核小體之間的空間位阻效應(yīng)[12]。而在體內(nèi)，核小體則與一些DNA結(jié)合蛋白競爭結(jié)合基因組DNA，可能會導致DNA序列信號在核小體定位中的作用受到不同程度的影響；而且染色質(zhì)重塑酶的作用也不可小覷，有時發(fā)揮單純的催化功能影響核小體的組裝效率，而有時ATP依賴的染色質(zhì)重塑酶能使核小體發(fā)生位移[21-22]。盡管DNA序列的內(nèi)在性質(zhì)和序列以外的因素(如染色質(zhì)重塑酶、蛋白因子與DNA序列的競爭結(jié)合等)在核小體定位中的重要性存在一定的爭議[23-24]，但體內(nèi)和體外核小體定位圖譜的高度相似性足以說明DNA序列是影響核小體定位的重要因素[12]。

過去幾十年間發(fā)展了不少核小體定位的理論預測模型，這些模型與核小體實驗圖譜相結(jié)合促進了核小體定位機制與功能的研究[11-12,18,25-27]。理論預測模型中最常見的是基于DNA序列的預測模型，而這類預測模型又可大致分為基于DNA序列的生物信息學模型[12,18,28-36]和基于DNA變形能的生物物理學模型[37-45]。生物信息學模型中利用機器學習算法構(gòu)建的模型不占少數(shù)。機器學習模型的建立與預測效果依賴于訓練集數(shù)據(jù)，而生物物理學模型的建立主要借助于DNA物理特性和核小體晶體結(jié)構(gòu)數(shù)據(jù)，并不需要訓練集。生物物理學方法能夠計算出DNA雙螺旋纏繞組蛋白八聚體的變形能，從而預測DNA序列形成核小體的能力、全基因組水平上的核小體占據(jù)率和核小體形成自由能[37-45]，也能夠預測出核小體在DNA序列上的較為準確的位置(或中心位置)[40]。核小體的準確位置涉及核小體的兩種定位方式，即平移定位和旋轉(zhuǎn)定位[46]，前者描述DNA序列與核小體核心區(qū)相對線性位置，而后者描述DNA雙螺旋與組蛋白八聚體相對方向。一般來講，旋轉(zhuǎn)定位信號強的DNA區(qū)域具有較高的彎曲各向異性，即容易朝某一個特定方向彎曲纏繞組蛋白八聚體形成核小體。顯然，旋轉(zhuǎn)定位信號強的區(qū)域也是平移定位較穩(wěn)定的區(qū)域。兩種定位方式緊密關(guān)聯(lián)，因此核小體在單堿基水平上的精確位置可借助旋轉(zhuǎn)定位信息來預測。本文介紹預測核小體定位的生物物理學方法及其應(yīng)用，旨在幫助人們更好地理解核小體定位，并建議選擇性地使用這些模型。

1 基于DNA變形能的核小體定位預測方法

結(jié)合核小體的晶體衍射結(jié)構(gòu)數(shù)據(jù)，可以計算任意一條147 bp DNA片段在形成核小體的假定下的DNA變形能，并以此判斷該DNA片段形成核小體的能力：變形能越小，DNA越容易彎曲，形成核小體的可能性越大。

為了計算變形能(或彈性能)，首先得科學描述DNA的幾何結(jié)構(gòu)。描述DNA構(gòu)象的方法主要包括(見表1)：自由連接鏈模型[47-48]、蠕蟲鏈模型[48-52]、圓形截面彈性桿模型[53-56]和堿基對梯階模型[57-58]。自由連接鏈模型視高分子為由自由鉸鏈(即連接處的彎曲方向與角度不受任何約束)連接多個獨立的剛性片段而形成的分子。雙鏈DNA分子中存在堿基對的氫鍵作用和堿基堆積作用，因此自由連接鏈模型不太適用于雙鏈DNA。該模型本身也有一些缺陷，導致其應(yīng)用受限[48]。蠕蟲鏈模型被認為是更加接近真實高分子的粗?；叻肿渔溎Ｐ停鋺?yīng)用較廣，但對于短于持久長度的DNA分子適用與否存在爭議[51,58]。國際上公認的另一種描述DNA結(jié)構(gòu)的方法是堿基對階梯模型(即劍橋協(xié)議方法[57])，即每一個堿基對上建立直角坐標系，并以相鄰的坐標系之間的三個線位移和三個角位移描述DNA雙螺旋結(jié)構(gòu)(見圖1)。研究DNA結(jié)構(gòu)的傳統(tǒng)彈性桿模型(如蟲鏈模型)中，將DNA看作是沿著序列連續(xù)彎曲的柔性桿，以DNA螺旋的骨架位形大致表示DNA的結(jié)構(gòu)。而劍橋協(xié)議方法則極大地豐富了DNA結(jié)構(gòu)的幾何學，促進了DNA柔性估計[58]、DNA結(jié)構(gòu)預測[59]、核小體定位[31,37,39-45]、啟動子識別[60-62]、剪接位點識別[63]、重組熱點識別[64-65]等多個與DNA結(jié)構(gòu)相關(guān)的生物學問題的研究(見表1)。

圖1 DNA雙螺旋結(jié)構(gòu)的堿基對梯階模型Fig.1 Base-pair step model for the DNA double helix structure

注：兩個平板表示兩個相鄰的堿基對，其中一個板相對于另一板繞三個軸旋轉(zhuǎn)的角位移(Tilt, Roll, Twist)和沿著三個軸方向平移的線位移(Shift, Slide, Rise)是模型中的六個結(jié)構(gòu)參數(shù).

表1 常用的DNA構(gòu)象描述方法與應(yīng)用

計算DNA變形能需預先知道描述二核苷酸彎曲能力的力常數(shù)。不同的二核苷酸有不同的力常數(shù)，同一個二核苷酸的不同的結(jié)構(gòu)變化(如Slide, Shift, Rise, Roll, Tilt, Twist的變化)對應(yīng)不同的力常數(shù)。力常數(shù)的計算主要有三種方法：基于DNA解鏈溫度的測定[37-38]、基于DNA結(jié)構(gòu)數(shù)據(jù)的計算[40,43,66-67]和分子動力學模擬方法[45,68]，下面重點介紹其中最常用的一種，即基于DNA結(jié)構(gòu)數(shù)據(jù)的力常數(shù)計算[67]。

令θn代表劍橋協(xié)議下參數(shù)符號，n=1,2,3,4,5,6分別對應(yīng)Ω、ρ、τ、DX、Dy和DZ。利用晶體結(jié)構(gòu)數(shù)據(jù)庫(如DNA-蛋白復合物的晶體結(jié)構(gòu)數(shù)據(jù))[69]計算每一種二核苷酸的結(jié)構(gòu)參數(shù)相對于所有二核苷酸的結(jié)構(gòu)參數(shù)平均值的漲落：

Δθn=θn-θ0

力常數(shù)的逆與結(jié)構(gòu)參數(shù)漲落的協(xié)方差之間的關(guān)系為[52]：

[F-1]nm=<ΔθnΔθm>/kT

(1)

其中k為玻爾茲曼常數(shù)，T為室溫(單位為開爾文)，<ΔθnΔθm>=<θnθm>-<θn><θm>。

根據(jù)式(1)計算力常數(shù)矩陣。力常數(shù)矩陣中對角元表示各種結(jié)構(gòu)參數(shù)對應(yīng)的力常數(shù)，而非對角元反映同一種二核苷酸的不同結(jié)構(gòu)參數(shù)間的耦合關(guān)聯(lián)，是屬于二階微量。

彈性桿模型中，任一147 bp長度的DNA序列的彈性能可近似用平衡態(tài)附近的二次勢能函數(shù)表示為[52]：

(2)

Tolstorokov等利用該模型預測了體外組裝的4個核小體在DNA序列上的位置，發(fā)現(xiàn)能量最低值(能量越低，越容易形成核小體)與核小體位置非常吻合[42]。他們的研究又表明roll和slide是決定核小體DNA超螺旋結(jié)構(gòu)的最主要參數(shù)。

Morozov等提出的預測核小體定位的彈性能模型中[43]，核小體DNA可從最初的超螺旋結(jié)構(gòu)進一步調(diào)整(微調(diào))。故核小體DNA的自由能由兩部分組成，

E=Eel+wEsh

(3)

其中Eel表示序列特異的DNA彈性能，Esh表示組蛋白與DNA的相互作用能，即由組蛋白與DNA的相互作用導致的核小體DNA相對于理想超螺旋結(jié)構(gòu)的偏離所對應(yīng)的能量。在此基礎(chǔ)上定義核小體起始概率和占據(jù)率[43]，成功預測了6個體外核小體的中心位置(能量的最小值較好地吻合核小體的中心位置)。

Miele等的基于能量的核小體定位預測模型[41]是依據(jù)與核小體晶體衍射結(jié)構(gòu)(1kx5)最吻合的理想超螺旋結(jié)構(gòu)建立的。假定DNA序列是不可伸長且不發(fā)生切變的彈性棒，而核小體形成時DNA的變形只由twist, roll和tilt造成，變形能表示為：

(4)

(5)

用該模型預測酵母和果蠅核小體占據(jù)率時，能夠準確預測轉(zhuǎn)錄起始等調(diào)控區(qū)域的核小體缺乏區(qū)，但對整個基因組范圍的核小體占據(jù)率的預測精度不是很高(如酵母第3號染色體的核小體占據(jù)率的預測結(jié)果與實驗數(shù)據(jù)[70]的相關(guān)系數(shù)R=0.45,P<10-15)。

De Santis等提出的預測核小體形成自由能的統(tǒng)計熱力學模型[37-38]與Miele的模型[41]相似，同樣是計算twist, roll和tilt對應(yīng)的自由能，但不包含核小體形成對應(yīng)的熵變這一項。用該模型預測的一些酵母基因組區(qū)域的形成核小體的自由能與實驗測定的核小體占據(jù)率相吻合，而且預測的100個核小體DNA的自由能與實驗測定的自由能之間有很強的正相關(guān)性(R=0.92，P<0.001)[37]。

Deniz等建立的DNA變形能模型[45]，形式上與式(2)相同，但其物理意義與前幾種模型不同。該模型計算的變形能是用來表示使自由DNA片段的結(jié)構(gòu)變?yōu)楹诵◇wDNA結(jié)構(gòu)的變形能。這里，所謂的核小體DNA的結(jié)構(gòu)指的是描述核小體DNA結(jié)構(gòu)的每個堿基對梯階(Base-pair step)對應(yīng)的六個坐標自由度(見圖1)，這些坐標通過對多個X射線晶體衍射核小體DNA結(jié)構(gòu)進行平均后獲得。而初始的平衡態(tài)自由DNA結(jié)構(gòu)則通過對少量的雙鏈DNA(但包含所有約化的10種二核苷酸)在水環(huán)境中的分子動力學模擬獲得。相關(guān)的力常數(shù)通過計算平衡態(tài)自由DNA結(jié)構(gòu)參數(shù)的協(xié)方差獲得(注：力常數(shù)矩陣是協(xié)方差矩陣的逆)。利用該模型，作者發(fā)現(xiàn)核小體缺乏區(qū)域的DNA變形能明顯高于其側(cè)翼序列，說明該區(qū)域從DNA序列的物理特性對核小體的形成有重要影響。

我們的模型[40,71-74]中主要考慮DNA彎曲和切變對應(yīng)的變形能，稱之為彎曲能和切變能。

纏繞組蛋白八聚體時DNA的彎曲主要取決于roll和tilt。假定導致DNA彎曲的扭力Fb均勻分布在DNA鏈上，則roll和tilt角的偏離平衡態(tài)的程度用下式表示：

(6)

其中i表示堿基對出現(xiàn)的位置，ρ0表示DNA的平均轉(zhuǎn)角，τ0表示DNA的平均傾角，kρ和kτ分別是轉(zhuǎn)角和傾角對應(yīng)的力常數(shù)，均由DNA-蛋白復合物的晶體結(jié)構(gòu)數(shù)據(jù)獲得[40]。Ωi表示從核小體中心位置(二分軸)向兩側(cè)計算的累加的扭角，其中每種堿基對梯階對應(yīng)的扭角是來自大量DNA-蛋白復合物的平均扭角。其實，所有堿基對梯階的扭角均取為w(i) =360°/10.4(核小體結(jié)構(gòu)NCP147的平均扭角)對結(jié)果的影響微乎其微。

在第i堿基對梯階中DNA的彎曲能為:

(7)

長度為L堿基對的DNA片段彎曲時，總的彎曲能為:

(8)

其中扭力Fb通過以下約束條件獲得[9,40]。在核小體上，扭力Fb的作用下129 bp長度的DNA(注：核小體DNA兩端各9 bp的區(qū)域是相對直的[9]，計算彎曲能時不考慮)在組蛋白八聚體上纏繞579度。而這579度(α)的彎曲(或纏繞)是由描述DNA雙螺旋結(jié)構(gòu)的ρ和τ角共同造成的。因此形成核小體的約束條件為:

(9)

結(jié)合式(6)可得

(10)

同理，結(jié)合核小體DNA的另一結(jié)構(gòu)約束條件(螺距)可得DNA切變能。

基于變形能可用玻爾茲曼分布近似估計序列形成核小體的潛能[40]。更普遍的做法是，將每一條染色體上的核小體分布看作是多個全同粒子(組蛋白八聚體)在一長鏈DNA分子上的分布，以變形能為基礎(chǔ)，用巨正則系綜理論計算出核小體在每一個147 bp序列片段上形成的概率，再通過概率加和計算出DNA位點上的核小體占據(jù)率[40,43]。具體計算方法如下：

(11)

其中conf表示與DNA結(jié)合的非重疊粒子的任一構(gòu)象，μ表示化學勢，E(conf)和n(conf)分別表示任一構(gòu)象對應(yīng)的總能量和總粒子數(shù)。

巨配分函數(shù)由一系列正向配分函數(shù)(從DNA序列的某一端開始計算)的遞進求解計算：

(12)

用同樣的方式計算反向的配分函數(shù)(從DNA序列的另一端開始計算)：

(13)

粒子(核小體)從第j個位點起始的概率(考慮了空間位阻效應(yīng))為:

(14)

第j個位點被核小體占據(jù)的概率為:

(15)

核小體占據(jù)率的計算以統(tǒng)計物理學巨正則系綜理論為基礎(chǔ)，考慮相鄰核小體空間位阻效應(yīng)。核小體占據(jù)率的另一種估算方法以Percus方程為基礎(chǔ)[75]。若s位點的核小體DNA的變形能為E(s)，分析核小體的組裝(即組蛋白八聚體沿DNA鏈的組裝)過程時，可將染色質(zhì)模擬為在外勢場E(s)中的有限長度l(核小體核心DNA長度，約等于147)的一維桿系統(tǒng)(流體)。熱力學巨正則系統(tǒng)描述中，系統(tǒng)被視為處于溫度項為β的熱浴環(huán)境和化學勢為μ的組蛋白八聚體熱庫中，系統(tǒng)在外勢場E(s)中達到熱力學平衡時其密度服從非線性積分方程(即Percus方程)：

(16)

獲得核小體密度的基礎(chǔ)上利用窗口大小為147 bp的矩形函數(shù)計算核小體占據(jù)率：

OCC(s)=ρ(s)·Π147(s)

(17)

我們的變形能模型能較好地預測體外核小體占據(jù)率、體外組裝的核小體在DNA序列上的準確位置、以及體外組裝核小體的自由能(即核小體的穩(wěn)定性)[40,72]。尤其是DNA彎曲能在核小體準確位置[40]和核小體滑動模式[72]的研究中有很好的應(yīng)用前景。模型中用到的描述DNA結(jié)構(gòu)的堿基對梯階參數(shù)也能較好地區(qū)分核小體富含區(qū)和核小體缺乏區(qū)[76]。

還有一些模型中除了DNA序列依賴的簡諧能以外還考慮了組蛋白和DNA之間的物理接觸位點上的相互作用[77]。另外，有一些模型中雖然用到DNA螺旋結(jié)構(gòu)參數(shù)[30]或其他DNA物理特征[31]，抑或是涉及DNA彎曲度[33]，但這些模型本質(zhì)上是屬于生物信息學方法，因為模型中主要是利用從序列提取的特征信息或DNA物理特征預測核小體的位置，而不涉及DNA變形能的計算。

2 討 論

預測核小體定位的生物信息學方法和生物物理學方法的主要區(qū)別在于：(1)生物信息學方法通常是使用大量的可靠數(shù)據(jù)來訓練模型[28-30]，但生物物理學方法是基于DNA的物理化學性質(zhì)(如二核苷酸的彎曲特性等)[37-45]；(2)由于體內(nèi)核小體定位還與其它非DNA因素有關(guān)，而且這種非DNA因素也可能是物種特異的(如物種特異的核小體定位模體)，基于不同物種核小體數(shù)據(jù)訓練的生物信息學方法的預測結(jié)果可能會優(yōu)于單純基于DNA物理化學性質(zhì)的生物物理學方法；(3)生物物理學方法能夠很好地預測核小體的中心位置及其可能的旋轉(zhuǎn)定位[40]，但只有少量的、設(shè)計巧妙的生物信息學方法才能做到這一點[34]；(4)與生物信息學方法相比，生物物理學方法的物理意義更加明了，有助于理解問題本質(zhì)。總的來說，生物信息學方法和生物物理學方法各有利弊。

各種生物物理學方法的主要差異包括：(1)使用的力常數(shù)不同，如有的用DNA-蛋白復合物的結(jié)構(gòu)數(shù)據(jù)基礎(chǔ)上計算的力常數(shù)[40,43]，而有的用二核苷酸的解鏈溫度表征其力常數(shù)[37]；(2)使用的DNA結(jié)構(gòu)參數(shù)不同，如有的用twist, roll, tilt, 有的用twist, roll, slide, 有的用roll, tilt, slide, 而有的用所有6種結(jié)構(gòu)參數(shù)；(3)預測核小體形成能力的最終指標中所包含的成分不同，如有的包含變形能和熵變，有的包含變形能和DNA-組蛋白相互作用能，有的只有變形能成分；(4)使用的核小體模型不同，如有的用核小體核心顆粒的真實DNA結(jié)構(gòu)模型[43]，而有的用與核小體核心顆粒DNA擬合最好的理想的超螺旋結(jié)構(gòu)模型。

我們認為，預測核小體定位時有以下問題值得注意：(1)力常數(shù)、平衡結(jié)構(gòu)參數(shù)的估算準確與否直接影響模型的預測結(jié)果，而這些參數(shù)的估計中需要注意DNA螺旋不同結(jié)構(gòu)參數(shù)之間的耦合相互作用[78-79]，計算原理的可靠性(如基于核小體DNA結(jié)構(gòu)的實驗數(shù)據(jù)的力常數(shù)計算和基于分子動力學模擬的力常數(shù)計算)和計算用的實驗數(shù)據(jù)量(如DNA-蛋白復合物或核小體DNA的晶體結(jié)構(gòu)數(shù)據(jù))；(2)預測能力從以下三個角度評價：核小體占據(jù)率的預測、核小體中心位置的預測、核小體裝配自由能的預測和核小體移動位置的預測，而不能只看其中一方面；(3)基于DNA序列的核小體定位預測模型應(yīng)該用體外核小體定位序列訓練模型，預測結(jié)果應(yīng)與體外核小體定位數(shù)據(jù)相比較，這樣能夠挖掘核小體定位對DNA序列的依賴本質(zhì)；(4)核小體定位預測模型的最終目標應(yīng)該是準確預測體內(nèi)核小體的位置，因為只有這樣才能使我們的預測本身更具有生物學意義。

總體而言，基于DNA序列預測酵母核小體定位的生物物理學方法取得了很大的進展，而且該類方法相對機器學習類生物信息學方法而言其物理意義清晰，但仍存在一些問題有待解決，例如：模型中的參數(shù)強烈依賴于DNA序列片段的物理性質(zhì)，不同長度、不同序列模式其物理意義可能存在很大的差異，這些都會直接地影響模型中參數(shù)的準確估計以及模型的應(yīng)用效果；模型中的有些近似假設(shè)需要更堅實的依據(jù)；不同的模型其側(cè)重點不同，適用問題(如核小體占據(jù)率的預測、核小體在序列上準確位置的預測、核小體穩(wěn)定性的預測以及核小體形成能力的預測)也略有不同；特定環(huán)境或過程中(如RNA聚合酶競爭性結(jié)合、組蛋白的化學修飾、染色質(zhì)重塑等)核小體定位會發(fā)生變化，這需要更專一的生物物理模型才能回答。