基于Hi-c數(shù)據(jù)的酵母染色體三維結(jié)構(gòu)重構(gòu)

豐繼華， 牟 錦， 郭亞茹

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

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

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

1 理論與方法

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

(1)

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

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

首先，需要對根據(jù)酵母染色體交互數(shù)據(jù)建立統(tǒng)計分布模型，為此，分別對酵母16條染色體的Hi-c數(shù)據(jù)分布情況進(jìn)行高斯擬合，對每條染色體的數(shù)據(jù)我們都分別與高斯8個線性組合核函數(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

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

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

(2)

1.2 染色體三維模型建立

在似然估計中，用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)

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

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

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

梯度上升迭代優(yōu)化中。利用等式 (3)計算偏導(dǎo)數(shù)，再根據(jù)偏導(dǎo)采用梯度上升優(yōu)化算法對各坐標(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)表示的是在時間步長t處參數(shù)Si的似然梯度。因此，式(4)中的隨機(jī)梯度上升可以用式(6)表示。Si是S中的參數(shù)。

gt,i=(Si(t))

(5)

(6)

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

(7)

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

(8)

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

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

2 模型評價

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

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

定義為：

(9)

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

定義為：

(10)

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

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

根據(jù)酵母16條染色體的Hi-c數(shù)據(jù)建立特定的分布函數(shù)，在此基礎(chǔ)上構(gòu)建目標(biāo)函數(shù)，然后利用梯度上升算法對每條染色體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á)到收斂，說明該目標(biāo)函數(shù)模型是有效可行的。文中使用梯度上升優(yōu)化算法對酵母16條染色體數(shù)據(jù)進(jìn)行三維空間結(jié)構(gòu)重構(gòu)，其模型的評價指標(biāo)如表3所示。

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

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

4 結(jié)論與展望

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