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

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

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

0 引言

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

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

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

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

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

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

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

不同的細(xì)胞格類型對(duì)應(yīng)有不同的組蛋白修飾。對(duì)于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 項(xiàng)目中公開(kāi)訪問(wèn)，首先對(duì)組蛋白修飾的數(shù)據(jù)進(jìn)行預(yù)處理，保持與Hi-C 數(shù)據(jù)同樣分辨率大?。ɡ绱笮?0kb）。

1.2 實(shí)驗(yàn)環(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 模型與評(píng)估

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

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

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

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

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

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

SSIM 值的范圍從0 到1，衡量?jī)煞鶊D的相似度，判斷預(yù)測(cè)圖像是否接近于原始圖像：

2 方法

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

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

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

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

圖1 Hi-C交互作用熱圖

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

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

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

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

圖3 兩段bins輸入

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

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

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

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

4 結(jié)語(yǔ)

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

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

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

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

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