融合卷積神經(jīng)網(wǎng)絡(luò)和最短路徑計(jì)算的染色體三維重構(gòu)

趙鳳嬌,鐘 誠(chéng)

(廣西大學(xué) 計(jì)算機(jī)與電子信息學(xué)院,南寧 530004)(廣西高校并行分布式計(jì)算技術(shù)重點(diǎn)實(shí)驗(yàn)室,南寧 530004) E-mail：1913392073@st.gxu.edu.cn

1 引 言

越來(lái)越多的研究表明,基因組的結(jié)構(gòu)和功能是相互統(tǒng)一的,染色體的三維結(jié)構(gòu)對(duì)基因表達(dá)調(diào)控等多種細(xì)胞功能非常重要,重建染色體的三維結(jié)構(gòu)對(duì)研究基因活動(dòng)的潛在機(jī)制具有重要意義[1].染色體構(gòu)象捕獲(3C)技術(shù)[2],特別是與新一代測(cè)序技術(shù)結(jié)合的高通量染色體構(gòu)象捕獲(Hi-C)技術(shù)[3],為染色體三維結(jié)構(gòu)的研究提供了全基因組的基因位點(diǎn)相互作用頻率(IF)數(shù)據(jù).

目前已有不少方法使用Hi-C數(shù)據(jù)進(jìn)行染色體三維重構(gòu),其中一些是基于距離約束實(shí)現(xiàn)的方法[4].文獻(xiàn)[5]建立了一個(gè)概率模型,由Hi-C數(shù)據(jù)確定空間距離,然后使用馬爾可夫鏈蒙特卡羅方法從數(shù)據(jù)中生成一個(gè)代表性樣本,以重建染色體三維結(jié)構(gòu).文獻(xiàn)[6]將染色體三維重構(gòu)表述為半定規(guī)劃問題,并利用黃金分割法確定出將IF值轉(zhuǎn)換為空間距離的合適參數(shù),以實(shí)現(xiàn)重建染色體的三維結(jié)構(gòu).文獻(xiàn)[7]通過泊松分布計(jì)算IF與空間距離之間的關(guān)系,采用多模型預(yù)測(cè)控制方法推斷染色體三維結(jié)構(gòu).以上幾種方法可以重建出質(zhì)量相對(duì)較好的染色體三維結(jié)構(gòu),然而由于算法模型的復(fù)雜性,它們只能有效處理幾百個(gè)基因位點(diǎn)大小的IF矩陣,無(wú)法處理更高分辨率的IF數(shù)據(jù).文獻(xiàn)[8]使用基于構(gòu)象能量和流形學(xué)習(xí)的方法來(lái)預(yù)測(cè)染色體的三維結(jié)構(gòu).文獻(xiàn)[9]將多維縮放(MDS)方法[10]與高斯公式推斷的相鄰基因位點(diǎn)局部相關(guān)性相結(jié)合,從染色體的IF矩陣中推斷出染色體三維結(jié)構(gòu).文獻(xiàn)[11]利用貝葉斯理論[12]定義染色體三維結(jié)構(gòu)空間的條件概率分布,使用分層優(yōu)化策略找出滿足給定約束條件下概率最高的染色體三維結(jié)構(gòu).ShRec3D方法[13]使用最短路徑算法預(yù)處理填充稀疏的IF矩陣,然后使用MDS方法重建染色體的三維結(jié)構(gòu).MDSGA方法[14]改進(jìn)了ShRec3D算法,使用遺傳算法優(yōu)化最短路徑推斷的距離,然后使用MDS方法重構(gòu)染色體三維結(jié)構(gòu).雖然與ShRec3D算法相比,MDSGA方法提高了重構(gòu)的三維結(jié)構(gòu)的精度,然而大大增加了算法所需的時(shí)間.

卷積神經(jīng)網(wǎng)絡(luò)(CNN)是深度學(xué)習(xí)的主流模型之一,已在圖像處理、計(jì)算機(jī)視覺和自然語(yǔ)言處理等領(lǐng)域得到廣泛應(yīng)用[15].為了在高分辨率染色體數(shù)據(jù)上有效重構(gòu)出精度更高的染色體三維結(jié)構(gòu),本文研究提出一種融合CNN和最短路徑計(jì)算的染色體三維重構(gòu)方法.余下的第2節(jié)詳細(xì)闡述重構(gòu)算法;第3節(jié)給出本文算法和同類具有代表性算法在Hi-C數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果;第4節(jié)總結(jié)本文的工作并給出下一步研究方向.

2 方 法

本文將Hi-C技術(shù)產(chǎn)生的染色體相互作用頻率(IF)數(shù)據(jù)構(gòu)成一個(gè)n×n規(guī)模的矩陣A,其中n代表基因位點(diǎn)的數(shù)量,矩陣元素Aij保存的是第i個(gè)基因位點(diǎn)和第j個(gè)基因位點(diǎn)之間的IF值.Hi-C實(shí)驗(yàn)得到的IF數(shù)據(jù)組成的矩陣常常是稀疏的[16],為了提高染色體三維重構(gòu)精度,首先需要降低IF矩陣的稀疏程度.為此,本文首先通過CNN建模,將原始稀疏IF矩陣推斷形成相對(duì)稠密IF矩陣,然后將相對(duì)稠密IF矩陣轉(zhuǎn)換為距離矩陣,依據(jù)距離矩陣構(gòu)造加權(quán)無(wú)向圖,并計(jì)算圖中任意兩頂點(diǎn)之間的最短距離,在此基礎(chǔ)上運(yùn)用多維縮放方法在完整距離矩陣上實(shí)現(xiàn)染色體三維重構(gòu).下面詳細(xì)闡述本文給出的融合卷積神經(jīng)網(wǎng)絡(luò)和最短路徑計(jì)算的染色體三維重構(gòu)方法.

2.1 推斷相對(duì)稠密IF矩陣

為了從染色體的原始IF矩陣推斷出相對(duì)稠密IF矩陣,本文構(gòu)建的CNN模型在訓(xùn)練中將原始IF矩陣作為已知的相對(duì)稠密IF矩陣,將使用隨機(jī)下采樣方法處理原始IF矩陣得到的結(jié)果作為稀疏IF矩陣.

(1)

本文參考文獻(xiàn)[15]搭建CNN框架,所構(gòu)建的染色體三維重構(gòu)CNN框架由3種不同類型的卷積層組成,具體結(jié)構(gòu)如圖1所示.

圖1 染色體三維重構(gòu)的卷積神經(jīng)網(wǎng)絡(luò)模型框架Fig.1 Convolution neural network model framework for 3D chromosome reconstruction

圖1中CNN框架的第1層用于提取和表示染色體稀疏IF矩陣上的特征,第2層用于調(diào)整矩陣的大小,第3層用于將提取到的特征非線性映射到相對(duì)稠密IF矩陣上,然后輸出推斷結(jié)果.CNN框架中的修正線性單元(ReLU)用于非線性激活函數(shù),殘差單元用于防止隨著網(wǎng)絡(luò)深度的增加,出現(xiàn)退化現(xiàn)象.

本文使用Tensorflow實(shí)現(xiàn)CNN模型,權(quán)重參數(shù)根據(jù)Xavier初始化方法進(jìn)行初始化,使用隨機(jī)梯度下降法對(duì)卷積核權(quán)值進(jìn)行更新.

2.2 推斷完整距離矩陣

染色體IF矩陣的大部分有效數(shù)據(jù)都集中在矩陣對(duì)角線附近,而偏離對(duì)角線的遠(yuǎn)距離數(shù)據(jù)則幾乎沒有.雖然使用CNN建?？梢院芎玫赝茢喑鋈旧wIF矩陣中對(duì)角線附近的缺失值,但是偏離對(duì)角線的遠(yuǎn)距離位置的數(shù)據(jù)依然很稀疏.為此,本文將相對(duì)稠密IF矩陣轉(zhuǎn)換為距離矩陣,然后通過最短路徑算法計(jì)算得到染色體的完整距離矩陣.

IF矩陣中的相互作用頻率越大,兩個(gè)對(duì)應(yīng)的基因位點(diǎn)間的歐氏距離就越小.設(shè)Bn×n表示相對(duì)稠密IF矩陣,其中Bij代表第i個(gè)基因位點(diǎn)和第j個(gè)基因位點(diǎn)間的IF值;Cn×n表示矩陣B對(duì)應(yīng)的距離矩陣,其中Cij代表第i個(gè)基因位點(diǎn)和第j個(gè)基因位點(diǎn)間的歐氏距離;i=1,2,…,n,j=1,2,…,n.本文通過翻轉(zhuǎn)相對(duì)稠密IF矩陣得到距離矩陣：

(2)

當(dāng)Bij為0時(shí),表示第i個(gè)基因位點(diǎn)和第j個(gè)基因位點(diǎn)之間的IF數(shù)據(jù)缺失,此時(shí)Bij對(duì)應(yīng)的Cij為未知距離,用-1表示.距離矩陣對(duì)角線元素的值表示某個(gè)基因位點(diǎn)和自身之間的距離,其值為0.

若距離矩陣中某一行或某一列除對(duì)角線外其余位置的元素值都為-1,則表示Hi-C實(shí)驗(yàn)沒有收集到該行或該列對(duì)應(yīng)基因位點(diǎn)的有效數(shù)據(jù),這樣的基因位點(diǎn)對(duì)染色體三維重構(gòu)沒有幫助,因此將它的數(shù)據(jù)從矩陣中移除.本文使用移除無(wú)效數(shù)據(jù)后得到的距離矩陣構(gòu)建加權(quán)無(wú)向圖,距離矩陣中的基因位點(diǎn)作為圖的頂點(diǎn),兩個(gè)基因位點(diǎn)之間的距離作為圖中對(duì)應(yīng)頂點(diǎn)之間邊的權(quán)值.若兩個(gè)基因位點(diǎn)之間的距離值為-1,則圖中對(duì)應(yīng)頂點(diǎn)間無(wú)邊.初始時(shí),若兩個(gè)頂點(diǎn)vi和vj之間有邊,則以邊的權(quán)值作為vi和vj之間的最短路徑長(zhǎng)度;若它們之間無(wú)邊,則其最短路徑長(zhǎng)度為∞.然后利用最短路徑算法求解更新任意一對(duì)頂點(diǎn)vi和vj之間的最短路徑,i=1,2,…,m,j=1,2,…,m,m≤n.本文將兩兩頂點(diǎn)之間的最短路徑長(zhǎng)度作為對(duì)應(yīng)基因位點(diǎn)在三維結(jié)構(gòu)中的歐氏距離,以推斷出染色體的完整距離矩陣.

2.3 推斷染色體三維結(jié)構(gòu)

獲得染色體的完整距離矩陣后,本文使用多維縮放(MDS)方法進(jìn)行染色體三維重構(gòu).MDS方法可以將數(shù)據(jù)樣本從高維空間轉(zhuǎn)換到低維空間,同時(shí)確保低維空間中任意兩個(gè)樣本之間的歐氏距離與對(duì)應(yīng)高維空間中的歐氏距離盡可能的相等.設(shè)染色體的完整距離矩陣為Dm×m,其中Dij代表第i個(gè)基因位點(diǎn)和第j個(gè)基因位點(diǎn)之間的距離,重構(gòu)的染色體結(jié)構(gòu)三維坐標(biāo)矩陣為F3×m,其中F·i代表第i個(gè)基因位點(diǎn)的三維坐標(biāo)(F1i,F2i,F3i);F·j代表第j個(gè)基因位點(diǎn)的三維坐標(biāo)(F1j,F2j,F3j),則矩陣F和矩陣D中的元素值應(yīng)盡可能滿足以下約束條件：

Dij=‖F(xiàn)·i-F·j‖

(3)

‖F(xiàn)·i-F·j‖表示F·i和F·j之間的歐氏距離.

設(shè)S為F的內(nèi)積矩陣,即S=FTF,由式(3)可以推導(dǎo)得到矩陣S和矩陣D滿足如下關(guān)系[10]：

(4)

求出Sij之后,i=1,2,…,m,j=1,2,…,m,進(jìn)而求解S的3個(gè)最大的特征值λ1、λ2、λ3以及對(duì)應(yīng)的3個(gè)特征向量X1、X2、X3,構(gòu)成對(duì)角矩陣Λ=diag(λ1,λ2,λ3)和特征向量矩陣X=(X1,X2,X3),則可以求出染色體結(jié)構(gòu)的三維坐標(biāo)矩陣F=Λ1/2XT.推斷出三維坐標(biāo)矩陣F后,即可使用繪圖工具繪制出染色體的三維結(jié)構(gòu).

2.4 算法描述

算法1形式描述了本文給出的融合卷積神經(jīng)網(wǎng)絡(luò)和最短路徑計(jì)算的染色體三維重構(gòu)算法(簡(jiǎn)記為CS3dC算法).

算法1.CS3dC

輸入：Hi-C數(shù)據(jù)矩陣Hl×l;染色體原始IF矩陣An×n,l、n為矩陣中染色體基因位點(diǎn)的數(shù)量

輸出：染色體三維坐標(biāo)矩陣F3×m,m為有矩陣中染色體效基因位點(diǎn)的數(shù)量

Begin

1.從染色體Hi-C數(shù)據(jù)集中學(xué)習(xí)染色體IF矩陣的特征,訓(xùn)練出具有利用IF矩陣元素值推斷出其鄰近區(qū)域元素值能力的CNN模型;

2.使用CNN模型根據(jù)染色體原始IF稀疏矩陣An×n推斷形成相對(duì)稠密矩陣Bn×n;

3. fori=0 ton-1 do //從矩陣B計(jì)算距離矩陣C

4. forj=0 ton-1 do

5. ifi=jthen

6.C[i][j]←0; //距離矩陣對(duì)角線值為0

7. else ifB[i][j]>0 then

8.C[i][j]←1/B[i][j];

9. else

10.C[i][j]←(-1);

11. end if

12. end for

13.end for

14.fori=0 ton-1 do //移除距離矩陣C中無(wú)效的行和列

15. is_empty←True;

16. forj=0 ton-1 do

17. ifC[i][j]>0 then

18. is_empty←False;

19. break;

20. end if

21. end for

22. if is_empty=True then

23. 移除矩陣C的第i行和第i列;

24. end if

25.end for

26.fork=0 tom-1 do //推斷出完整距離矩陣D

27. fori=0 tom-1 do

28. forj=0 tom-1 do

29. ifD[i][j]>C[i][k]+C[k][j]then

30.D[i][j]←C[i][k]+C[k][j];

31. end if

32. end for

33. end for

34.end for

38.fori=0 tom-1 do //計(jì)算內(nèi)積矩陣S的值

39. forj=0 tom-1 do

41. end for

42.end for

43.解方程|S-λE|=0,獲得矩陣S的3個(gè)最大的特征值λ1、λ2和λ3,構(gòu)成對(duì)角矩陣Λ=diag(λ1,λ2,λ3);

44.將λ1、λ2和λ3代入線性方程(Sx-λx)=0,解得特征向量X1、X2和X3,組成特征向量矩陣X=(X1,X2,X3);

45.計(jì)算三維坐標(biāo)矩陣F=Λ1/2XT;

End

算法CS3dC的步1用于訓(xùn)練CNN模型,所需時(shí)間為O(n2);步2使用CNN模型推斷相對(duì)稠密矩陣Bn×n,所需時(shí)間為O(n2);步3～步13使用相對(duì)稠密矩陣B計(jì)算出距離矩陣Cn×n,所需時(shí)間為O(n2);步14～步25用于找出并移除距離矩陣C中無(wú)效的行和列,所需時(shí)間為O(n2);步26～步34執(zhí)行從距離矩陣Cm×m推斷得到完整距離矩陣Dm×m,所需時(shí)間為O(m3);步35所需時(shí)間為O(m2),步36和步37所需時(shí)間均為O(m);步38～步42計(jì)算內(nèi)積矩陣S的值,所需時(shí)間為O(m2);步43和步44所需時(shí)間均為O(m2);步45計(jì)算三維坐標(biāo)矩陣F,所需時(shí)間為O(m2).因此,算法CS3dC的時(shí)間復(fù)雜度為O(n2+m3)=O(n3),m≤n.算法CS3dC運(yùn)行過程中使用染色體Hi-C數(shù)據(jù)矩陣Hl×l、原始IF矩陣An×n、相對(duì)稠密矩陣Bn×n、距離矩陣Cn×n、完整距離矩陣Dm×m、內(nèi)積矩陣Sm×m以及三維坐標(biāo)矩陣F3×m,因此算法所需的空間復(fù)雜度為O(l2+n2+m2+m)=O(l2+n2),m≤n.

本文算法CS3dC通過融合卷積神經(jīng)網(wǎng)絡(luò)和最短路徑計(jì)算從染色體原始IF矩陣推斷出完整距離矩陣,降低了染色體矩陣的稀疏程度,減少了使用MDS算法重構(gòu)染色體三維過程中不確定值的數(shù)量,從而提高了重構(gòu)結(jié)果的精度.

3 實(shí) 驗(yàn)

實(shí)驗(yàn)在廣西大學(xué)Sugon 7000A超級(jí)并行計(jì)算機(jī)系統(tǒng)(1)https：//hpc.gxu.edu.cn/上進(jìn)行,使用CX50-G30計(jì)算節(jié)點(diǎn),其CPU為2×Intel Xeon Gold 6230,共40核,配置192GB DDR4內(nèi)存,操作系統(tǒng)為64位版本的CentOS 7.4,采用Python3.7編程實(shí)現(xiàn)算法.

本文使用的實(shí)驗(yàn)數(shù)據(jù)來(lái)自高通量基因表達(dá)數(shù)據(jù)庫(kù)(2)https：//www.ncbi.nlm.nih.gov/geo/中登記號(hào)為GSE63525的數(shù)據(jù)集,該數(shù)據(jù)集中包含多種不同類型細(xì)胞的Hi-C相互作用頻率矩陣數(shù)據(jù)[17].本文選用規(guī)模為342.46MB的人類B淋巴細(xì)胞數(shù)據(jù)集GM12878中的1～17號(hào)染色體數(shù)據(jù)作為CNN建模的訓(xùn)練集,18號(hào)染色體數(shù)據(jù)作為驗(yàn)證集,19～22號(hào)染色體數(shù)據(jù)作為測(cè)試集;選用規(guī)模為342.43MB的人類紅白血病細(xì)胞K562數(shù)據(jù)集、規(guī)模為342.46MB的人類胚肺成纖維細(xì)胞IMR90數(shù)據(jù)集和規(guī)模為260.55MB的小鼠B淋巴細(xì)胞CH12-LX數(shù)據(jù)集測(cè)試算法的染色體三維重構(gòu)效果.以上所有數(shù)據(jù)集分辨率均為100Kb,比對(duì)質(zhì)量(MAPQ)≥30.

實(shí)驗(yàn)測(cè)試了本文算法CS3dC與具有代表性的同類算法ShRec3D[13]、MDSGA[14]的均方根誤差值(RMSE)、Pearson相關(guān)系數(shù)值(PCC)和所需的運(yùn)行時(shí)間,其中RMSE值越小、PCC值越大,則說明重建得到的染色體三維結(jié)構(gòu)與真實(shí)結(jié)構(gòu)越接近,重建效果越好.

3.1 訓(xùn)練卷積神經(jīng)網(wǎng)絡(luò)的實(shí)驗(yàn)

本文參照文獻(xiàn)[18]的設(shè)置方式,將CNN訓(xùn)練的隨機(jī)下采樣比率設(shè)置為1/16、稀疏IF矩陣子矩陣規(guī)模的參數(shù)p設(shè)置為40、相對(duì)稠密IF矩陣子矩陣規(guī)模的參數(shù)q設(shè)置為28.本文設(shè)置CNN訓(xùn)練的學(xué)習(xí)率為0.001,每次處理的樣本大小為256,最大訓(xùn)練次數(shù)為200,當(dāng)達(dá)到最大訓(xùn)練次數(shù)或者驗(yàn)證集上的損失不再發(fā)生變化時(shí)停止訓(xùn)練.

本文使用Pearson相關(guān)系數(shù)值(PCC)作為CNN模型訓(xùn)練效果的測(cè)試指標(biāo),CNN訓(xùn)練過程中驗(yàn)證集上的的損失值(Loss)和PCC值變化情況如圖2所示.

圖2 CNN在驗(yàn)證集上的Loss值和PCC值Fig.2 Lossand PCCof CNN on the validation set

由圖2可知,隨著訓(xùn)練次數(shù)的增加,驗(yàn)證集上的Loss值先是急劇下降,然后下降趨勢(shì)變緩,在下降過程中有降有升,但總體趨勢(shì)是下降;驗(yàn)證集上的PCC值先是快速增大,然后增速變緩,最終趨近于1.這表明在訓(xùn)練過程中,開始時(shí)CNN的收斂速度很快,模型的訓(xùn)練效果也有很大的提升,然后CNN的收斂速度降低,模型的訓(xùn)練效果提升也變得緩慢,最終CNN停止收斂,此時(shí)CNN模型的訓(xùn)練效果也達(dá)到最佳.

本文將訓(xùn)練好的CNN模型應(yīng)用到染色體三維重構(gòu)中,用于從染色體稀疏的原始IF矩陣中推斷相對(duì)稠密IF矩陣.

3.2 人類細(xì)胞數(shù)據(jù)集上的實(shí)驗(yàn)

表1 算法在K562數(shù)據(jù)集上運(yùn)行得到的RMSE值和PCC值Table 1 RMSEand PCCobtained by running three algorithms on K562 data set

表2 算法在K562數(shù)據(jù)集上運(yùn)行所需的時(shí)間(s)Table 2 Required time(s) of running three algorithms on K562 data set

首先測(cè)試算法ShRec3D、MDSGA和本文算法CS3dC在人類紅白血病細(xì)胞K562數(shù)據(jù)集上的RMSE值、PCC值和所需的運(yùn)行時(shí)間,實(shí)驗(yàn)結(jié)果分別如表1和表2所示.

表1表明,對(duì)于K562數(shù)據(jù)集的1～22號(hào)染色體,本文算法CS3dC獲得的RMSE值最低,算法MDSGA的RMSE值次之,算法ShRec3D的RMSE值最高;對(duì)于K562數(shù)據(jù)集的1～22號(hào)染色體,本文算法CS3dC的PCC值最大,算法ShRec3D的次之,算法MDSGA的最小.這是因?yàn)樗惴⊿hRec3D僅使用最短路徑方法推算距離矩陣中的缺失值,算法MDSGA使用遺傳算法進(jìn)一步優(yōu)化距離矩陣中的最短路徑距離,而本文算法CS3dC融合卷積神經(jīng)網(wǎng)絡(luò)和最短路徑計(jì)算推斷完整距離矩陣,從而使得完整距離矩陣中的數(shù)據(jù)有更高的精度,進(jìn)而在染色體三維重構(gòu)過程中始終能夠獲得更低的RMSE值和更大的PCC值.

由表2可知,在3種算法中,對(duì)于K562數(shù)據(jù)集的1～20號(hào)染色體,本文算法CS3dC的運(yùn)行時(shí)間最少;對(duì)于20～22號(hào)染色體,算法ShRec3D運(yùn)行時(shí)間最少,其次是本文算法CS3dC的運(yùn)行時(shí)間,算法MDSGA的運(yùn)行時(shí)間最多.從整體上看,本文算法CS3dC的運(yùn)行時(shí)間平均值最少,為50.0秒;算法ShRec3D的運(yùn)行時(shí)間平均值次之,為67.5秒;算法MDSGA的運(yùn)行時(shí)間平均值最多,達(dá)到8686.8秒.

為了進(jìn)一步測(cè)試算法在其他人類細(xì)胞數(shù)據(jù)集上運(yùn)行的性能,實(shí)驗(yàn)比對(duì)了算法ShRec3D、MDSGA和本文算法CS3dC在人類胚肺成纖維細(xì)胞IMR90數(shù)據(jù)集上的RMSE值、PCC值和所需的運(yùn)行時(shí)間,結(jié)果如表3和表4所示.

從表3可以看到,本文算法CS3dC的RMSE值在所有染色體上都顯著低于其他兩種算法MDSGA和ShRec3D;與算法MDSGA和ShRec3D相比,本文算法CS3dC的PCC值更大.這表明在IMR90數(shù)據(jù)集上仍然是本文算法CS3dC的染色體三維重構(gòu)效果最好,與真實(shí)情況最接近.

表3 算法在在IMR90數(shù)據(jù)集上運(yùn)行得到的RMSE值和PCC值Table 3 RMSEand PCCobtained by running three algorithms on IMR90 data set

表4 算法在IMR90數(shù)據(jù)集上所需的運(yùn)行時(shí)間(s)Table 4 Required time(s) of running three algorithms on MR90 data set

表4的結(jié)果表明,對(duì)于IMR90數(shù)據(jù)集的1～17號(hào)染色體,本文算法CS3dC所需的運(yùn)行時(shí)間最少,算法ShRec3D的次之,MDSGA的最多;對(duì)于18～22號(hào)染色體,算法ShRec3D所需的運(yùn)行時(shí)間最少,本文算法CS3dC的次之,算法MDSGA的最多.對(duì)于1～22號(hào)染色體,本文算法CS3dC的運(yùn)行時(shí)間平均值最低,為49.1秒;算法ShRec3D的運(yùn)行時(shí)間平均值次之,為58.9秒;算法MDSGA的運(yùn)行時(shí)間平均值最多,高達(dá)8708.5秒.

綜上所述,無(wú)論是在K562數(shù)據(jù)集上還是在IMR90數(shù)據(jù)集上,本文算法CS3dC都能夠在所需運(yùn)行時(shí)間平均值最少的情況下,獲得比算法ShRec3D和MDSGA更好的染色體三維重構(gòu)效果.

3.3 小鼠B淋巴細(xì)胞數(shù)據(jù)集上的實(shí)驗(yàn)

為了測(cè)試算法在其他不同物種細(xì)胞數(shù)據(jù)集上的染色體三維結(jié)構(gòu)重構(gòu)效果,本文使用小鼠B淋巴細(xì)胞CH12-LX數(shù)據(jù)集測(cè)試算法ShRec3D、MDSGA和本文算法CS3dC的RMSE值、PCC值和所需的運(yùn)行時(shí)間,實(shí)驗(yàn)結(jié)果分別如表5和表6所示.

由表5可知,對(duì)于CH12-LX數(shù)據(jù)集的1～19號(hào)染色體,本文算法CS3dC獲得的RMSE值均是3種算法中最低的,且所有值均小于1,而算法MDSGA和ShRec3D的值均大于1;對(duì)于CH12-LX數(shù)據(jù)集的1～19號(hào)染色體,本文算法CS3dC的PCC值一直是3種算法中最高的.這表明,本文算法CS3dC在小鼠物種數(shù)據(jù)上的染色體三維重構(gòu)效果明顯優(yōu)于其他2種算法ShRec3D和MDSGA.

表6的結(jié)果表明,對(duì)于1～15號(hào)染色體和17、18號(hào)染色體,本文算法CS3dC所需的運(yùn)行時(shí)間最少,對(duì)于16、19號(hào)染色體,算法ShRec3D所需的運(yùn)行時(shí)間最少.在整個(gè)數(shù)據(jù)集上,本文算法CS3dC的運(yùn)行時(shí)間平均值為35.4秒,算法ShRec3D的運(yùn)行時(shí)間平均值為44.5秒,算法MDSGA的運(yùn)行時(shí)間平均值為2131.6秒.因此,從整體上看,本文算法CS3dC在小鼠B淋巴細(xì)胞CH12-LX數(shù)據(jù)集上所需的運(yùn)行時(shí)間是最少的.

上述結(jié)果表明,本文算法CS3dC可以在較短的時(shí)間內(nèi)重構(gòu)出精度更高的小鼠B淋巴細(xì)胞染色體三維結(jié)構(gòu).

表5 算法在CH12-LX數(shù)據(jù)集上運(yùn)行得到的RMSE值和PCC值Table 5 RMSEand PCCobtained by running three algorithgs on CH12-LX data set

表6 算法在CH12-LX數(shù)據(jù)集上所需的運(yùn)行時(shí)間(s)Table 6 Required time(s) of running three algorithms on CH12-LX data set

圖3 算法重構(gòu)的染色體三維結(jié)構(gòu)效果Fig.3 Effect of chromosomes 3D structure reconstructed by three algorithms

綜合3.2節(jié)和3.3節(jié)的實(shí)驗(yàn)結(jié)果可知,與算法ShRec3D和MDSGA相比,本文算法CS3dC不僅可以在人類細(xì)胞數(shù)據(jù)集上也可以在小鼠細(xì)胞數(shù)據(jù)集上,快速地重構(gòu)出更準(zhǔn)確的染色體三維結(jié)構(gòu).這說明本文提出的融合卷積神經(jīng)網(wǎng)絡(luò)和最短路徑計(jì)算的染色體三維重構(gòu)算法能夠應(yīng)用于多種物種,具有通用性.

3.4 重構(gòu)的染色體三維結(jié)構(gòu)可視化

在Hi-C實(shí)驗(yàn)中,基因位點(diǎn)傾向于優(yōu)先與同一拓?fù)湎嚓P(guān)結(jié)構(gòu)域(TAD)內(nèi)的基因位點(diǎn)發(fā)生相互作用,這表明同一TAD內(nèi)的基因位點(diǎn)在染色體三維結(jié)構(gòu)中距離更近[19].本文依據(jù)在人類細(xì)胞數(shù)據(jù)集和小鼠細(xì)胞數(shù)據(jù)集上的實(shí)驗(yàn)數(shù)據(jù)結(jié)果,使用Python的繪圖庫(kù)Matplotlib繪制染色體的三維結(jié)構(gòu).以K562的3號(hào)染色體、IMR90的10號(hào)染色體和CH12-LX的18號(hào)染色體為例,算法ShRec3D、MDSGA和本文算法CS3dC重構(gòu)的三維結(jié)構(gòu)可視化結(jié)果如圖3所示.

從圖3(a)可以看到,K562數(shù)據(jù)集3號(hào)染色體的IF矩陣熱圖包含兩個(gè)明顯的TAD,本文算法CS3dC重構(gòu)的染色體三維結(jié)構(gòu)分為兩個(gè)獨(dú)立的區(qū)域(圖中黑線兩側(cè)所示);與IF矩陣熱圖一致,而算法ShRec3D和MDSGA重構(gòu)的三維結(jié)構(gòu)中的兩個(gè)部分區(qū)分比較模糊(圖中黑線兩側(cè)所示),沒有很好地體現(xiàn)出此染色體結(jié)構(gòu)的特點(diǎn).

圖3(b)表明,IMR90數(shù)據(jù)集10號(hào)染色體的IF矩陣熱圖由兩個(gè)TAD組成,本文算法CS3dC重構(gòu)的染色體三維結(jié)構(gòu)也可以劃分為兩個(gè)結(jié)構(gòu)(圖中黑線兩側(cè)所示),與IF矩陣熱圖具有對(duì)應(yīng)關(guān)系,而算法ShRec3D和MDSGA重構(gòu)的三維結(jié)構(gòu)區(qū)分度很低,可視化效果較差.

從圖3(c)也可以看出,CH12-LX數(shù)據(jù)集18號(hào)染色體的IF矩陣熱圖包含3個(gè)較大的TAD,本文算法CS3dC重構(gòu)的染色體三維結(jié)構(gòu)分為3個(gè)緊密的結(jié)構(gòu)(圖中黑線兩側(cè)所示),與染色體IF矩陣熱圖表現(xiàn)出的特征相似,而算法ShRec3D和MDSGA的重構(gòu)結(jié)果難以區(qū)分,未能有效表示此染色體的三維結(jié)構(gòu).

上述結(jié)果表明,相比于算法ShRec3D和MDSGA,本文算法CS3dC重構(gòu)得到的染色體三維結(jié)構(gòu)與染色體的IF矩陣熱圖具有更好的對(duì)應(yīng)關(guān)系,它與真實(shí)的染色體三維結(jié)構(gòu)具有更高的一致性.

4 結(jié) 語(yǔ)

本文設(shè)計(jì)的染色體三維結(jié)構(gòu)重構(gòu)算法的特色和新穎之處是：通過卷積神經(jīng)網(wǎng)絡(luò)建模處理染色體稀疏的原始相互作用頻率矩陣數(shù)據(jù),并融合最短路徑計(jì)算推斷出完整距離矩陣,降低了染色體數(shù)據(jù)的稀疏程度,然后運(yùn)用多維縮放方法實(shí)現(xiàn)推斷出精度更高的染色體三維結(jié)構(gòu),同時(shí)算法所需的運(yùn)行時(shí)間較短,而且能夠應(yīng)用于多種物種的染色體三維重構(gòu).下一步的工作方向是探索采用其他機(jī)器學(xué)習(xí)方法,例如強(qiáng)化學(xué)習(xí)和聯(lián)邦學(xué)習(xí)方法,建模處理染色體數(shù)據(jù)是否能獲得更好的三維結(jié)構(gòu)重構(gòu)效果.