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

豐繼華， 牟 錦， 郭亞茹

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

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

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

1 理論與方法

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

(1)

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

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

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

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

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

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

(2)

1.2 染色體三維模型建立

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

(3)

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

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

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

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

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

(4)

式(4)中t是迭代索引指標，S(t)是迭代索引指標t的結構坐標，λ(t)是在t處的學習速率，隨著迭代的進行可能發(fā)生變化，L(S(t))是結構中坐標的似然偏導數(shù)。式(5)表示的是在時間步長t處參數(shù)Si的似然梯度。因此，式(4)中的隨機梯度上升可以用式(6)表示。Si是S中的參數(shù)。

gt,i=(Si(t))

(5)

(6)

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

(7)

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

(8)

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

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

2 模型評價

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

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

定義為：

(9)

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

定義為：

(10)

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

3 實驗結果

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

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

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

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

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

4 結論與展望

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