基于深度學(xué)習(xí)的染色質(zhì)交互作用預(yù)測

任景瑞，李川，張振毓，鄧凱

（1.四川大學(xué)計算機(jī)學(xué)院，四川610065；2.四川大學(xué)網(wǎng)絡(luò)空間安全學(xué)院，四川610065）

0 引言

染色質(zhì)是由DNA、組蛋白、非組蛋白等多種物質(zhì)組成的遺傳物質(zhì)，其結(jié)構(gòu)復(fù)雜，難以直接觀察，但對細(xì)胞遺傳過程的基因表達(dá)有重要影響。自3C 技術(shù)問世以來，眾多方法被陸續(xù)報道用于捕獲染色質(zhì)構(gòu)象，其中Hi-C 技術(shù)是捕獲染色質(zhì)相互作用頻次的最新最常用方法[1]。Hi-C 原始互作數(shù)據(jù)可以通過交互頻次的讀取序列映射到對稱矩陣中，并且利用這種矩陣熱圖可以表示并構(gòu)造為染色質(zhì)的高級結(jié)構(gòu)TAD[2]、隔間和染色質(zhì)環(huán)等。染色質(zhì)的高級結(jié)構(gòu)與其功能密切相關(guān)，對基因表達(dá)和生物遺傳有重要影響[3]，如三維結(jié)構(gòu)變化可能誘導(dǎo)腫瘤發(fā)育產(chǎn)生[4]。

目前對染色質(zhì)結(jié)構(gòu)的研究非常豐富，由于Hi-C 數(shù)據(jù)測量非常昂貴且耗時，但對與基因表達(dá)、轉(zhuǎn)錄和疾病狀態(tài)相關(guān)的染色質(zhì)結(jié)構(gòu)有重要意義[5]，所以有很多研究開始關(guān)注減少實驗進(jìn)行預(yù)測，對染色質(zhì)結(jié)構(gòu)研究存在三維結(jié)構(gòu)構(gòu)建，二維結(jié)構(gòu)預(yù)測以及基因表達(dá)和功能性研究。對染色質(zhì)三維構(gòu)建目前存在多種方法，分別使用了多種距離模型算法來構(gòu)建三維空間結(jié)構(gòu)[6]。最新的三維結(jié)構(gòu)構(gòu)建方法還可以分析其基因功能[7]。而二維結(jié)構(gòu)分析同樣基于Hi-C 數(shù)據(jù)區(qū)域分割[8]或DNA 序列預(yù)測[9]來構(gòu)建TAD、染色質(zhì)環(huán)等結(jié)構(gòu)，識別其區(qū)域性。染色質(zhì)測序技術(shù)的發(fā)展，還對表觀遺傳學(xué)方面的實驗應(yīng)用非常重要[10]，表觀遺傳學(xué)包括組蛋白修飾等方面，對基因表達(dá)調(diào)控和染色質(zhì)重塑有重要影響[11]。分析組蛋白修飾的功能作用[12]，研究染色質(zhì)結(jié)構(gòu)功能和基因表達(dá)[13]的影響，例如對染色質(zhì)開放性[14]和染色質(zhì)狀態(tài)的識別[15]，具有重要研究價值。

利用深度學(xué)習(xí)預(yù)測染色質(zhì)交互作用能夠有效獲取輸入數(shù)據(jù)的前后關(guān)聯(lián)信息和局部特征，目前對人類細(xì)胞類型GM12878[16]和果蠅細(xì)胞類型[17]都有相關(guān)報道，使用的數(shù)據(jù)一般圍繞DNA 序列和表觀遺傳學(xué)數(shù)據(jù)。為了深入研究組蛋白修飾和染色質(zhì)交互作用兩種數(shù)據(jù)相關(guān)性，本文提出了一種從組蛋白修飾數(shù)據(jù)中預(yù)測人類基因組中Hi-C 數(shù)據(jù)的方法，基于卷積神經(jīng)網(wǎng)絡(luò)，針對常見人類細(xì)胞類型IMR90 建立了深度學(xué)習(xí)模型比較預(yù)測，通過線性相關(guān)系數(shù)皮爾遜系數(shù)以及圖相似性系數(shù)等評估，并最終在預(yù)測結(jié)果與原始結(jié)果之間表現(xiàn)出高相關(guān)性。

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)預(yù)處理

Hi-C 數(shù)據(jù)：在GEO 公開數(shù)據(jù)集上，可以通過訪問代碼GSE63525 下載IMR90 細(xì)胞類型的數(shù)據(jù)。我們從實驗原始觀測的Hi-C 序列交互數(shù)據(jù)生成原始矩陣數(shù)據(jù)，并根據(jù)分辨率確定每段基因的長度，計算對應(yīng)段位置（例如本文采用10kb 分辨率，k 為一千數(shù)量，b 指代堿基段即一段堿基序列，就會將每個交互位置數(shù)據(jù)除以104，獲得其bin 序號，每條染色體按細(xì)胞類型和染色體不同有上百萬或上億堿基長度）。本文關(guān)注研究染色質(zhì)內(nèi)交互作用，即同一序號染色體間的交互作用（同染色質(zhì)間交互），最終獲得22 條染色體的Hi-C 原始交互作用熱圖（去除性染色體影響）。

組蛋白修飾數(shù)據(jù)：從Roadmap 上可以下載每種細(xì)胞類型的所有表觀遺傳修飾數(shù)據(jù)，https://egg2.wustl.edu/roadmap/data/byFileType/signal/consolidated/macs2signal/foldChange/。

不同的細(xì)胞格類型對應(yīng)有不同的組蛋白修飾。對于IMR90 細(xì)胞類型，可以下載以下多種修飾因素：

H3K23ac、H3K79Me1、H3K27Ac、H3K79me2、H3K27me2、H3K79me3、H3K27me3、H3K9acH3K3K36me1、H3K9me1、HP4、RPD3、H1、H3K36me2、H3K36me3、H3K9me3、H3K-4me1、H4、H3、H3、3、3me1、H4K3、3k4me1、H4K16ac

以上數(shù)據(jù)可在ENCODE 項目中公開訪問，首先對組蛋白修飾的數(shù)據(jù)進(jìn)行預(yù)處理，保持與Hi-C 數(shù)據(jù)同樣分辨率大小（例如大小為10kb）。

1.2 實驗環(huán)境

GPU：NVIDIA TU102[GeForce RTX 2080 Ti Rev.A]（rev a1）

CPU：48 英特爾至強(qiáng)CPU E5-2650 v4@2.20GHz

內(nèi)存：128GB

Python 版本和依賴包環(huán)境：Python 3.6，基于TensorFlow 的Keras。

1.3 模型與評估

本文通過多種的深度學(xué)習(xí)方法評估訓(xùn)練結(jié)果。訓(xùn)練神經(jīng)網(wǎng)絡(luò)并使其擬合的過程重，使用MSE（Mean Square Error）作為損失函數(shù)，使用MAE（Mean Absolute Error）作為目標(biāo)函數(shù)。

線性相關(guān)性上使用皮爾遜相關(guān)系數(shù)（PCC）和斯皮爾曼相關(guān)系數(shù)（SPCC）來評估預(yù)測結(jié)果，結(jié)果數(shù)值將顯示預(yù)測結(jié)果與原始數(shù)據(jù)之間的線性相關(guān)性，相關(guān)系數(shù)的絕對值越接近1，其相關(guān)性就越強(qiáng)。

皮爾遜相關(guān)系數(shù)公式：

SPCC 是基于PCC 的一種相關(guān)系數(shù)計算方法，給定變量基礎(chǔ)上給出一個排序差異集合d，由兩個變量集合中的每個元素計算，最后使用PCC 公式來獲取排序變量的結(jié)果。

除序列評估外，還可使用計算峰值信噪比（Peak Signal to Noise Ratio）和結(jié)構(gòu)相似度指數(shù)（Structural SIMilarity）來進(jìn)行評估，這兩者都常用于圖像處理和去噪。

PSNR 表示圖像信噪比，此值越大代表失真越少，MAX=max{Ypred}-min{Ypred}。

SSIM 值的范圍從0 到1，衡量兩幅圖的相似度，判斷預(yù)測圖像是否接近于原始圖像：

2 方法

2.1 數(shù)據(jù)預(yù)處理

對于輸入數(shù)據(jù)和輸出數(shù)據(jù)，本文分別使用不同方法進(jìn)行預(yù)處理，因為不同數(shù)據(jù)實驗測序方法不同，其表達(dá)值彼此差異較大，難以直接定性分析。

Hi-C 數(shù)據(jù)可使用標(biāo)準(zhǔn)化函數(shù)將其歸一化到[0,1]的區(qū)間范圍，表示每個位置交互的可能性，其原始序列交互數(shù)據(jù)可從Rao 等人（GSE63525）[18]所做的公開可用的Hi-C 實驗中獲得。原始數(shù)據(jù)為每個染色質(zhì)每段堿基部分間的交互作用強(qiáng)度，在10k 堿基分辨率下，每個作用強(qiáng)度數(shù)據(jù)就表示染色質(zhì)上按順序排列的兩段10k 長度堿基的交互頻率。由于每個染色質(zhì)長度彼此不同，使用字母標(biāo)識i 和j 表示兩個段堿基序號，它們之間的計數(shù)nij 表示染色質(zhì)段上相互作用次數(shù)，通過堿基位置對應(yīng)關(guān)系可構(gòu)成對稱矩陣，矩陣大小為N 表示染色質(zhì)長度L/分辨率R（如圖1 所示）。

對于組蛋白數(shù)據(jù)，首先通過bwtools 和指定bin 長度（分辨率）生成組蛋白修飾序列信號數(shù)據(jù)矩陣H，此分辨率與Hi-C 數(shù)據(jù)分辨率一致，且進(jìn)行截取使實驗數(shù)據(jù)長度相同，矩陣H 中每列為組蛋白修飾類型，共M列，使用最大最小值歸一化使數(shù)據(jù)在[0,1]區(qū)間，每列數(shù)據(jù)表示組蛋白修飾與染色質(zhì)產(chǎn)生作用的可能性。

圖1 Hi-C交互作用熱圖

2.2 Hi-C預(yù)測模型

神經(jīng)網(wǎng)絡(luò)介紹：卷積神經(jīng)網(wǎng)絡(luò)（CNN）一般用于圖像處理或自然語言處理等高維特征自動提取，可學(xué)習(xí)到平移不變性等特征，在生物信息研究上也廣泛應(yīng)用。利用CNN 可以快速提取相應(yīng)染色質(zhì)交互的相關(guān)因子特征序列，構(gòu)建交互作用概率分布輸入。在全連接層部分，使用循環(huán)神經(jīng)網(wǎng)絡(luò)考慮染色質(zhì)上下游數(shù)據(jù)帶來的影響，并使所有神經(jīng)元之間的參數(shù)共享，這些參數(shù)針對染色質(zhì)基相互作用進(jìn)行了優(yōu)化，可生成用于不同組蛋白修飾的最佳濾波器，均方誤差（MSE）作為損失函數(shù)，平均絕對誤差（MAE）作為目標(biāo)函數(shù)。對輸入輸出數(shù)據(jù)采用取對數(shù)值獲取[-0,1]區(qū)間值域范圍，并作為神經(jīng)網(wǎng)絡(luò)模型的標(biāo)簽和樣本，可理解為對每段染色質(zhì)交互作用的可能性預(yù)測。

本文提出的模型基于每對序列對應(yīng)的方式，對矩陣中數(shù)據(jù)一一預(yù)測，如圖3 所示，使用一個w 大小的窗口來獲取每個交互基因距離下的所有數(shù)據(jù)，由于Hi-C矩陣中的對稱特性，交互數(shù)據(jù)可以按列或行獲取。因為染色質(zhì)的交互和高維結(jié)構(gòu)受到堿基段的上下游影響，有明顯的區(qū)域性，對于輸入的組蛋白修飾序列數(shù)據(jù)，每個交互位點的上下游各一個堿基段作為輸入，共三個堿基段長度。因此，每個樣本關(guān)注w 大小的染色質(zhì)交互作用，使用x-1 到x+1 段（x 為交互作用發(fā)生的位置）的3 個長度的組蛋白修飾作為輸入來預(yù)測相應(yīng)的Hi-C 交互作用情況。本文對模型輸入部分進(jìn)行劃分，構(gòu)建一個輸入模塊獲取兩對不同位置的輸入數(shù)據(jù)，此模型將在Hi-C 矩陣中的對角線區(qū)域附近生成每個bin 的交互作用數(shù)據(jù)。最后結(jié)果用熱圖重建方法來還原預(yù)測矩陣。窗口大小w 的選取，可使用的Hi-C 矩陣為10kb（104）分辨率下的實驗數(shù)據(jù)，設(shè)定w=50，因此基因組距離為500kb，即每個堿基段包含500k（500×103）的堿基。這樣可以觀察交互作用密集區(qū)域，排除交互作用發(fā)生不明顯的區(qū)域和較遠(yuǎn)距離的稀疏數(shù)據(jù)，使得預(yù)測結(jié)果更有價值。

圖2 染色質(zhì)交互作用預(yù)測神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)圖

圖3 兩段bins輸入

模型預(yù)測過程包括三個階段，第一階段為輸入數(shù)據(jù)的卷積和特征聚合階段，對每兩對輸入的組蛋白修飾信號矩陣，通過多層一維卷積神經(jīng)網(wǎng)絡(luò)獲取其多種類型修飾的一維聚合特征序列，代表不同修飾類型共同作用的綜合權(quán)重分布。第二階段為拼接層，將兩對輸入構(gòu)成二維矩陣，再次利用卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行不同位置間的交互作用影響特征提取。最后一層為全連接層和循環(huán)神經(jīng)網(wǎng)絡(luò)層，利用堿基上下游序列的影響，轉(zhuǎn)化為時序序列的前后特征，來預(yù)測相應(yīng)染色質(zhì)堿基段的交互作用結(jié)果。由于神經(jīng)網(wǎng)絡(luò)預(yù)測過程中是針對各個位點進(jìn)行的，需要根據(jù)位置和對稱性重建矩陣結(jié)果進(jìn)行對比評估，預(yù)測結(jié)果中可得到多段w 大小的預(yù)測數(shù)值排列，根據(jù)提取過程的順序可以依次對應(yīng)到交互作用發(fā)生的位置上，并使用重構(gòu)算法將其恢復(fù)成數(shù)值矩陣。因此，最終結(jié)果也是對角線區(qū)域在指示基因組距離上的完整矩陣，并存在一定的拓?fù)浣Y(jié)構(gòu)區(qū)域。

重構(gòu)矩陣熱圖偽代碼：

M 為最后結(jié)果矩陣

For i in N:

if i ＜N-w+1:M[i,i:i+w]=Input[1:w];M[i:i+w,i]=M[i,i:i+w]

else:x=N-I;M[i,i:i+w]=Input[1:x];M[i:i+w,i]=M[i,i:i+w]

3 結(jié)果

我們對模型進(jìn)行了多次訓(xùn)練，設(shè)定超參數(shù)為訓(xùn)練輪次設(shè)定為30 輪，批次大小每次100 個樣本，優(yōu)化器是RMSProp。訓(xùn)練完畢后可以獲得一維濾波器的各項權(quán)重，表示每對固定組蛋白修飾序列的綜合作用特征。最終結(jié)果分布表示序列相關(guān)性和矩陣相似性。

最終結(jié)果分別使用線性相關(guān)和圖像相似性進(jìn)行分析。數(shù)據(jù)集使用GEO 數(shù)據(jù)庫中的GSE63525 訪問代碼獲取的10kb 分辨率的IMR90 原始交互作用數(shù)據(jù)，重構(gòu)成交互作用矩陣。其中訓(xùn)練集使用1-17 號染色體的Hi-C 樣本，測試集使用18-22 號染色體的Hi-C 樣本。最后結(jié)果顯示，在500kb 的基因組距離下，對Hi-C 樣本進(jìn)行預(yù)測結(jié)果分析。線性相關(guān)性分析中測試集PCC 最好達(dá)到0.85，SPCC 最好達(dá)到0.8，SSIM 的測試樣本中最好可達(dá)到0.98。

表1 細(xì)胞類型IMR90 染色質(zhì)編號的測試數(shù)據(jù)集為染色體18 號-22 號/基因組距離（500kb）

4 結(jié)語

本文對染色質(zhì)交互作用的二維結(jié)構(gòu)與表觀遺傳學(xué)數(shù)據(jù)組蛋白修飾信號進(jìn)行了相關(guān)性預(yù)測研究，提出了一種基于組蛋白修飾信號序列數(shù)據(jù)交叉預(yù)測染色質(zhì)交互作用的方法。結(jié)果表明，組蛋白修飾信號在染色質(zhì)相互作用中可起到重要作用，也為染色質(zhì)結(jié)構(gòu)預(yù)測在深度學(xué)習(xí)和不同數(shù)據(jù)上提供了一種可行路線。

對于人類細(xì)胞類型數(shù)據(jù)，具有數(shù)據(jù)量大，交互作用矩陣數(shù)據(jù)稀疏，結(jié)構(gòu)作用域難以預(yù)測等特點，本文提出的方法針對稀疏數(shù)據(jù)進(jìn)行了基因距離篩選，一定程度解決了數(shù)據(jù)稀疏性和不相關(guān)數(shù)據(jù)干擾的問題，并對每對數(shù)據(jù)進(jìn)行分別預(yù)測，減少了不同樣本差異帶來的訓(xùn)練過擬合影響。高維染色質(zhì)結(jié)構(gòu)與染色質(zhì)組成的內(nèi)部物質(zhì)高度相關(guān)，可以通過其他組成數(shù)據(jù)進(jìn)行推測。

圖4 chr22 16mb-18mb 上的交互

圖5 chr22 20mb-22mb 上的交互

本文的深度學(xué)習(xí)模型，并應(yīng)用了不同的數(shù)據(jù)集和復(fù)雜的神經(jīng)網(wǎng)絡(luò)層結(jié)構(gòu)，從組蛋白修飾推測染色質(zhì)數(shù)據(jù)。結(jié)果表明了模型的能力以及組蛋白修飾的方向，如何影響染色質(zhì)組織。但方法使用數(shù)據(jù)類型不豐富，染色質(zhì)結(jié)構(gòu)不止與組蛋白修飾信號相關(guān)，輸入樣本數(shù)據(jù)特征仍然可以增加，模型預(yù)測數(shù)據(jù)在高維結(jié)構(gòu)的留存性上仍有待提高的。