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

1、基于深度學(xué)習(xí)的染色質(zhì)交互作用預(yù)測摘要：染色質(zhì)交互作用一般發(fā)生在生物基因表達(dá)過程中，可以對染色質(zhì)高維空間如染色質(zhì)環(huán)、拓?fù)浣Y(jié)構(gòu)域等產(chǎn)生影響，通過 特定實驗的高通量測序技術(shù)如5C、Hi-C等可獲得染色質(zhì)交互作用信息,預(yù)測染色質(zhì)交互作用信息具有重要價值,可 減少測序技術(shù)高昂的代價。提出基于組蛋白修飾數(shù)據(jù)來預(yù)測Hi-C數(shù)據(jù)的深度學(xué)習(xí)方法，從不同角度對結(jié)果進(jìn)行分 析，表明該方法可以較好預(yù)測原始染色質(zhì)交互作用情況，為組蛋白修飾預(yù)測染色質(zhì)交互作用數(shù)據(jù)提供一種可行路線。 關(guān)鍵詞：組蛋白修飾；Hi-C測序技術(shù)；深度學(xué)習(xí)；染色質(zhì)交互作用0引言染色質(zhì)是由DNA、組蛋白、非組蛋白等多種物質(zhì)組 成的遺傳物質(zhì)，其結(jié)構(gòu)復(fù)雜，

2、難以直接觀察，但對細(xì)胞 遺傳過程的基因表達(dá)有重要影響。自3C技術(shù)問世以 來，眾多方法被陸續(xù)報道用于捕獲染色質(zhì)構(gòu)象，其中 Hi-C技術(shù)是捕獲染色質(zhì)相互作用頻次的最新最常用 方法o Hi-C原始互作數(shù)據(jù)可以通過交互頻次的讀取 序列映射到對稱矩陣中,并且利用這種矩陣熱圖可以 表示并構(gòu)造為染色質(zhì)的高級結(jié)構(gòu)TAD121、隔間和染色質(zhì) 環(huán)等。染色質(zhì)的高級結(jié)構(gòu)與其功能密切相關(guān),對基因 表達(dá)和生物遺傳有重要影響131，如三維結(jié)構(gòu)變化可能誘 導(dǎo)腫瘤發(fā)育產(chǎn)生1 41 o目前對染色質(zhì)結(jié)構(gòu)的研究非常豐富，由于Hi-C數(shù) 據(jù)測量非常昂貴且耗時，但對與基因表達(dá)、轉(zhuǎn)錄和疾病 狀態(tài)相關(guān)的染色質(zhì)結(jié)構(gòu)有重要意義151，所以有很多

3、研究 開始關(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)問。最新 的三維結(jié)構(gòu)構(gòu)建方法還可以分析其基因功能，7|o而二 維結(jié)構(gòu)分析同樣基于Hi-C數(shù)據(jù)區(qū)域分割域DNA序 列預(yù)測，91來構(gòu)建TAD、染色質(zhì)環(huán)等結(jié)構(gòu)，識別其區(qū)域 性。染色質(zhì)測序技術(shù)的發(fā)展,還對表觀遺傳學(xué)方面的 實驗應(yīng)用非常重要1 101,表觀遺傳學(xué)包括組蛋白修飾等方 面,對基因表達(dá)調(diào)控和染色質(zhì)重塑有重要影響111。分析 組蛋白修飾的功能作用121，研究染色質(zhì)結(jié)構(gòu)功能和基因 表達(dá)131的影響，例如對染色質(zhì)開

4、放性141和染色質(zhì)狀態(tài)的 識別151，具有重要研究價值。利用深度學(xué)習(xí)預(yù)測染色質(zhì)交互作用能夠有效獲取 輸入數(shù)據(jù)的前后關(guān)聯(lián)信息和局部特征，目前對人類細(xì) 胞類型GM128781 161和果蠅細(xì)胞類型171都有相關(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

5、-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ì) 胞類

6、型的所有表觀遺傳修飾數(shù)據(jù)， HYPERLINK https:/egg2.wust https:/egg2.wust!. edu/roadmap/data/byFileType/signal/consolidated/ macs2signal/fO!dChange/。不同的細(xì)胞格類型對應(yīng)有不同的組蛋白修飾。對于 IMR90細(xì)胞類型，可以下載以下多種修飾因素：H3K23ac、H3K79Me1、H3K27Ac、H3K79me2、H3K27me2、 H3K79me3、H3K27me3、H3K9acH3K3K36me1、H3K9me1、 HP4、RPD3、H1、H3K36me2、H3K36me3、H3K

7、9me3、H3K- 4me 1、H4、H3、H3、3、3me 1、H4K3、3k4me 1、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)存:128GBPython版本和依賴包環(huán)境：Python 3.6,基于Ten- sorFlow 的 Keras。1.3模型與評估本文通過多種的深度學(xué)習(xí)方法評估訓(xùn)練結(jié)

8、果。訓(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)。N匕=1皮爾遜相關(guān)系數(shù)公式: ( YrealredYpred)SPCC是基于PCC的一種相關(guān)系數(shù)計算方法，給 定變量基礎(chǔ)上給出一個排序差異集合d,由兩個變量集 合中的每個元素計算，最后使用PCC公式來獲取排序 變量的結(jié)果。除序列評估外，還可使用計

9、算峰值信噪比(Peak Signal to Noise Ratio)和結(jié)構(gòu)相似度指數(shù)(Structural SIMilarity)來進(jìn)行評估，這兩者都常用于圖像處理和 去噪。PSNR表示圖像信噪比，此值越大代表失真越少， MAX=max Ypred min Ypred 。PSNR = 20*log10( MAX )4MSESSIM值的范圍從0到1,衡量兩幅圖的相似度，判 斷預(yù)測圖像是否接近于原始圖像：SSIM -QreidNpred + 1)(realpred + C2)(4 + 說red + 回(以 + 優(yōu)red + C2)2方法2.1數(shù)據(jù)預(yù)處理對于輸入數(shù)據(jù)和輸出數(shù)據(jù),本文分別使用不同方 法

10、進(jìn)行預(yù)處理，因為不同數(shù)據(jù)實驗測序方法不同，其表 達(dá)值彼此差異較大,難以直接定性分析。Hi-C數(shù)據(jù)可使用標(biāo)準(zhǔn)化函數(shù)將其歸一化到0,1 的 區(qū)間范圍，表示每個位置交互的可能性,其原始序列交 互數(shù)據(jù)可從Rao等人(GSE63525)1181所做的公開可用的 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ì)

11、 長度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交互作用熱圖圖1 Hi-C交互作用熱圖神經(jīng)網(wǎng)絡(luò)介紹：卷積神經(jīng)網(wǎng)絡(luò)（CNN）一般用于圖 像處理或自然語言處理等高維特征自動提取,可學(xué)習(xí) 到平移不變性等特征，在生物信息研究上也廣泛應(yīng) 用。利用CNN可以快速提取相應(yīng)染色質(zhì)交互的相關(guān) 因子特征序列，構(gòu)建交互作用概

12、率分布輸入。在全連 接層部分,使用循環(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ì)的交

13、互和高維結(jié)構(gòu)受到堿基段的上下游影 響，有明顯的區(qū)域性，對于輸入的組蛋白修飾序列數(shù) 據(jù),每個交互位點(diǎn)的上下游各一個堿基段作為輸入，共 三個堿基段長度。因此，每個樣本關(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,因此 基

14、因組距離為500kb,即每個堿基段包含500k （500X 103 ）的堿基。這樣可以觀察交互作用密集區(qū) 域,排除交互作用發(fā)生不明顯的區(qū)域和較遠(yuǎn)距離的稀 疏數(shù)據(jù),使得預(yù)測結(jié)果更有價值。Bini Input | Bin2 Input| ConvlD(k：rnel=32) Bini Input | Bin2 InputIBatchNormalizationI |BatchNormalization| : :| ConvlD (k：rne1=64廠| | ConvlD(kernel=64門 : : ； |ConvlD(kerne1=128) | |ConvlD(kernel=128)|ConvlD(

15、kernel=l)ConvlD(kernel=32)ConvlD(kernel=l)BatchNormalizationConvlD(kernel=l)ConvlD(kernel=32)ConvlD(kernel=l)ConvlD(kernel=64) BatchNormalization| ConvlD (kernel=128)|BatchNormalization| ConvlD (kernel=128)ConvlD (kernel=l)Dense(w*8)Dense(w*4)Dense(w*2)RNN(w*2)Dense(=w)Dense(=1)圖2染色質(zhì)交互作用預(yù)測神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)圖圖

16、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ù)測過程中是針對各 個位點(diǎn)進(jìn)行的，需要根據(jù)位置和對稱性重建矩陣結(jié)果 進(jìn)行對比評估，預(yù)測結(jié)果中可得到多段w大小的預(yù)測 數(shù)值排列

17、，根據(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: Mi, i: i+w = Input1:w; Mi: i+w, i=M i, i: i+welse: x = N - I; Mi, i: i+w = Input1:x; Mi: i+w, i=Mi, i: i+w3結(jié)果我們對模型進(jìn)行了多次訓(xùn)練，設(shè)定超參數(shù)為訓(xùn)練 輪次設(shè)定為30輪,批次大小每次100個樣本，優(yōu)化器 是RMSProp

18、。訓(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）染色體序號整體PCC整體SPCC平均SS