豐繼華， 牟 錦， 郭亞茹

(云南民族大學(xué) 電氣信息工程學(xué)院， 昆明 650500)

Dekker等人于2002年就利用3c技術(shù)[1]用于重構(gòu)酵母染色體的空間結(jié)構(gòu)，由于受限于當(dāng)時(shí)的技術(shù)條件，他們只進(jìn)行了染色體區(qū)域內(nèi)一對(duì)一的相互作用研究。近年來(lái)基于此技術(shù)發(fā)展出了4C[2]，5C[3]，Hi-c以及TCC技術(shù)[4]，為染色體三維結(jié)構(gòu)研究奠定了基礎(chǔ)。4C是將DNA利用連接酶連接成環(huán)狀后用特異PCR引物進(jìn)行反向PCR，再對(duì)其產(chǎn)物進(jìn)行分析得到染色體相互作用數(shù)據(jù)；5C技術(shù)是在3c基礎(chǔ)上增加了測(cè)序通量，主要用于研究染色體高頻率的空間結(jié)構(gòu)。而Hi-c技術(shù)[5]就是一種在3C基礎(chǔ)上發(fā)展的高通量測(cè)序染色體全基因組的交互數(shù)據(jù)的生物信息學(xué)技術(shù)。隨著基因組學(xué)的深入研究，人們發(fā)現(xiàn)染色體的相互作用數(shù)據(jù)在一定程度上可以反應(yīng)基因組在三維空間的表達(dá)狀況[6]。另一方面，相較于其他的3C衍生技術(shù)，利用Hi-c技術(shù)捕獲到的全基因組交互數(shù)據(jù)為利用二維接觸矩陣構(gòu)建三維空間結(jié)構(gòu)提供了可能。通過(guò)基因空間結(jié)構(gòu)解析與傳統(tǒng)轉(zhuǎn)錄數(shù)據(jù)相結(jié)合，研究人員可以更深入的闡釋生物在基因調(diào)控以及表觀遺傳中的真實(shí)性狀[7]。因此預(yù)測(cè)和重構(gòu)染色體的三維結(jié)構(gòu)對(duì)于后基因組時(shí)代研究具有指導(dǎo)意義。盡管目前部分染色體的組織結(jié)構(gòu)可以通過(guò)電子顯微鏡進(jìn)行研究。顯微鏡提供了單個(gè)細(xì)胞的信息，但分辨率相對(duì)較低；而染色體構(gòu)象捕獲能獲得高分辨率的染色體三維信息，極大拓展了我們對(duì)基因組的全面認(rèn)識(shí)。

文中使用Duan等人研究使用的酵母數(shù)據(jù)樣本[8]，根據(jù)酵母16條染色體Hi-c數(shù)據(jù)構(gòu)建了數(shù)據(jù)統(tǒng)計(jì)分布模型，在此基礎(chǔ)上利用梯度優(yōu)化算法預(yù)測(cè)并繪制了酵母染色體的三維結(jié)構(gòu)圖。

1 理論與方法

利用Hi-c技術(shù)獲得的基因組高通量染色體交互頻率數(shù)據(jù)，通過(guò)特定數(shù)學(xué)模型預(yù)測(cè)基因組空間結(jié)構(gòu)是重構(gòu)染色體三維結(jié)構(gòu)普遍采取的方法。其預(yù)測(cè)流程主要分為：Hi-c數(shù)據(jù)歸一化處理；Hi-c數(shù)據(jù)轉(zhuǎn)化為距離接觸矩陣；染色體三維重構(gòu)模型；和模型結(jié)果分析等。其中，Lieberman-Aiden等人曾對(duì)染色體上兩個(gè)片段的接觸頻率和基因線性以及空間距離做了開(kāi)創(chuàng)性研究，發(fā)現(xiàn)染色體片段之間的接觸頻率值與兩個(gè)片段之間空間的距離成反比關(guān)系，即空間越接近則接觸頻率值越大，空間距離越遠(yuǎn)接觸頻率越小[9]，在此原理上提出了以下距離轉(zhuǎn)換關(guān)系式：

(1)

式(1)中，Dij是表示酵母染色體上兩個(gè)片段之間的通過(guò)轉(zhuǎn)換的空間距離值，F(xiàn)ij表示酵母染色體片段間的接觸頻率值。

1.1 酵母染色體Hi-c數(shù)據(jù)分布擬合函數(shù)模型

首先，需要對(duì)根據(jù)酵母染色體交互數(shù)據(jù)建立統(tǒng)計(jì)分布模型，為此，分別對(duì)酵母16條染色體的Hi-c數(shù)據(jù)分布情況進(jìn)行高斯擬合，對(duì)每條染色體的數(shù)據(jù)我們都分別與高斯8個(gè)線性組合核函數(shù)進(jìn)行擬合，再最終選取出擬合指標(biāo)SSE,RMSE,R-square最優(yōu)的高斯核函數(shù)，最終選取核函數(shù)的擬合指標(biāo)結(jié)果如表1所示。

表1 16條染色體擬合情況表Table 1 Fitting of 16 chromosomes

在最終確定了每條染色體擬合出對(duì)應(yīng)的高斯核函數(shù)后，繪制了16條染色體Hi-c數(shù)據(jù)分布的擬合曲線(見(jiàn)圖1)。

通過(guò)與酵母16條染色體交互數(shù)據(jù)分布擬合后，獲得目標(biāo)函數(shù)如下：

(2)

1.2 染色體三維模型建立

在似然估計(jì)中，用S表示酵母染色體結(jié)構(gòu)，D表示從染色體交互數(shù)據(jù)導(dǎo)出的接觸矩陣,似然函數(shù)P(Di|S) 表示在結(jié)構(gòu)S條件下D中數(shù)據(jù)點(diǎn)概率[11]，在此，真實(shí)的Hi-c數(shù)據(jù)分布由擬合得到的組合高斯模型代替，因此P(Di|S)可以表示為：

(3)

我們的目的是找到一個(gè)最大化似然函數(shù)的結(jié)構(gòu)S*。式(3)中的目標(biāo)函數(shù)僅依賴染色體結(jié)構(gòu)中的(x,y,z)坐標(biāo)。

1.3梯度上升優(yōu)化算法

利用梯度上升算法對(duì)式(3)進(jìn)行迭代優(yōu)化，直到算法收斂為止。具體過(guò)程：如果使用新的(x,y,z)坐標(biāo)計(jì)算的似然函數(shù)值和前一步的差值小于一個(gè)閾值，就認(rèn)為算法收斂[11]。

梯度上升迭代優(yōu)化中。利用等式 (3)計(jì)算偏導(dǎo)數(shù)，再根據(jù)偏導(dǎo)采用梯度上升優(yōu)化算法對(duì)各坐標(biāo)進(jìn)行調(diào)整，并按下式更新似然概率：

S(t+1)=S(t)+λ(t)L(S(t))

(4)

式(4)中t是迭代索引指標(biāo)，S(t)是迭代索引指標(biāo)t的結(jié)構(gòu)坐標(biāo)，λ(t)是在t處的學(xué)習(xí)速率，隨著迭代的進(jìn)行可能發(fā)生變化，L(S(t))是結(jié)構(gòu)中坐標(biāo)的似然偏導(dǎo)數(shù)。式(5)表示的是在時(shí)間步長(zhǎng)t處參數(shù)Si的似然梯度。因此，式(4)中的隨機(jī)梯度上升可以用式(6)表示。Si是S中的參數(shù)。

gt,i=(Si(t))

(5)

(6)

在Ada Grad的迭代規(guī)則中，根據(jù)式(7)，在參數(shù)Si的之前計(jì)算的梯度基礎(chǔ)上，修正了每個(gè)參數(shù)Si在每一時(shí)間步長(zhǎng)上的學(xué)習(xí)速率 。

(7)

式中，Gt是一個(gè)對(duì)角元素為i的對(duì)角矩陣，i是Si的梯度平方和，如式(8)所示。其中Gt,ii是Gt中染色體片段i對(duì)應(yīng)的值，而ε是一個(gè)平滑項(xiàng)，它是避免函數(shù)除以零(通常是1×10-6)。

(8)

Ada Grad的主要優(yōu)點(diǎn)之一是它不需要在每次迭代時(shí)手動(dòng)調(diào)整學(xué)習(xí)速率。

表2 酵母染色體Hi-c數(shù)據(jù)分布模型參數(shù)Table 2 Parameters of yeast chromosome Hi-c data distribution model

2 模型評(píng)價(jià)

為了評(píng)估染色體三維結(jié)構(gòu)模型的準(zhǔn)確性，我們使用Pearson相關(guān)系數(shù)(PCC)、Spearman相關(guān)系數(shù)(SCC)這兩個(gè)參數(shù)作為評(píng)價(jià)指標(biāo)。假設(shè)兩個(gè)模型的成對(duì)距離數(shù)據(jù)集，其中{di,...,dn}有n個(gè)值，另一數(shù)據(jù)集{Di,...,Dn}也含n個(gè)值，那么DPCC、DSCC可以使用以下公式來(lái)計(jì)算。

(1)距離Pearson相關(guān)系數(shù)(DPCC)

定義為：

(9)

(2)距離Spearman相關(guān)系數(shù)(DSCC)

定義為：

(10)

DSCC測(cè)量了兩個(gè)三維結(jié)構(gòu)距離剖面的相似性。DSCC值在-1.0和1.0之間變化，DSCC值越高，這兩個(gè)結(jié)構(gòu)就越相似。

3 實(shí)驗(yàn)結(jié)果

根據(jù)酵母16條染色體的Hi-c數(shù)據(jù)建立特定的分布函數(shù)，在此基礎(chǔ)上構(gòu)建目標(biāo)函數(shù)，然后利用梯度上升算法對(duì)每條染色體Hi-c數(shù)據(jù)目標(biāo)函數(shù)進(jìn)行迭代，迭代的最大次數(shù)為2 000次，而收斂閾值設(shè)置為0.000 01，酵母16條染色體目標(biāo)函數(shù)收斂曲線如圖2所示。

圖2 酵母16條染色體目標(biāo)函數(shù)收斂曲線Fig.2 Convergence curve of 16 chromosome objective functions in yeast

從圖2中可以看出酵母16條染色體目標(biāo)函數(shù)最終都達(dá)到收斂，說(shuō)明該目標(biāo)函數(shù)模型是有效可行的。文中使用梯度上升優(yōu)化算法對(duì)酵母16條染色體數(shù)據(jù)進(jìn)行三維空間結(jié)構(gòu)重構(gòu)，其模型的評(píng)價(jià)指標(biāo)如表3所示。

表3 16條染色體結(jié)構(gòu)Spearman和Pearson系數(shù)Table 3 Spearman and Pearson coefficients of 16 chromosomal structures

從表中可以看出，經(jīng)過(guò)對(duì)每條染色體的Hi-c數(shù)據(jù)分布特征進(jìn)行擬合出具體函數(shù)作為目標(biāo)函數(shù)的模型的Spearman系數(shù)都達(dá)到了0.95以上，Pearson系數(shù)平均值也能達(dá)到0.71以上，說(shuō)明對(duì)每條染色體進(jìn)行Hi-c數(shù)據(jù)分布特征進(jìn)行擬合出不同目標(biāo)函數(shù)的方法來(lái)預(yù)測(cè)其結(jié)構(gòu)是有效可行的。通過(guò)不同目標(biāo)函數(shù)模型預(yù)測(cè)出的染色體結(jié)構(gòu)如圖3所示。

4 結(jié)論與展望

目前，利用染色體Hi-c交互數(shù)據(jù)預(yù)測(cè)三維空間結(jié)構(gòu)大多根據(jù)單一分布模型，并沒(méi)有考慮每條染色體數(shù)據(jù)的具體分布情況，而本文通過(guò)分析酵母16條染色體Hi-c數(shù)據(jù)的實(shí)際分布，從而擬合出更真實(shí)的分布模型，在此基礎(chǔ)上利用梯度優(yōu)化算法預(yù)測(cè)出較準(zhǔn)確的染色體三維結(jié)構(gòu)。但是為了分析方便，在對(duì)酵母16條染色體Hi-c數(shù)據(jù)進(jìn)行距離轉(zhuǎn)換時(shí)使用了統(tǒng)一的參數(shù)，后續(xù)我們將會(huì)針對(duì)具體染色體不同數(shù)據(jù)對(duì)轉(zhuǎn)換函數(shù)的參數(shù)進(jìn)行優(yōu)化，增強(qiáng)模型的自適應(yīng)性，從而進(jìn)一步提高模型預(yù)測(cè)的準(zhǔn)確性。