王新峰，黃 偉

1.吉首大學(xué) 軟件學(xué)院，湖南 吉首 416000

2.中山大學(xué) 計(jì)算機(jī)學(xué)院，廣州 510000

3.中南大學(xué) 計(jì)算機(jī)學(xué)院，長(zhǎng)沙 410000

DNA甲基化（DNA methylation）是表觀遺傳基因組中重要組成部分，它能在不改變DNA序列的前提下改變遺傳表現(xiàn)，也是研究最深入的表觀遺傳調(diào)控機(jī)制之一。DNA甲基化是指DNA序列上特定的堿基在DNA甲基轉(zhuǎn)移酶（DNA methyltransferase，DNMT）的催化作用下，以S-腺苷甲硫氨酸作為甲基供體，通過(guò)共價(jià)鍵結(jié)合的方式獲得一個(gè)甲基基團(tuán)的化學(xué)修飾過(guò)程[1]。它在基因表達(dá)和調(diào)控中起著重要作用，并參與許多細(xì)胞過(guò)程，包括細(xì)胞分化、發(fā)育和腫瘤發(fā)生。最常見(jiàn)的DNA甲基化類(lèi)型為5-甲基胞嘧啶（5mC）、6-甲基腺嘌呤（6mA）和4-甲基胞嘧啶（4mC），其中5mC是研究最廣泛的類(lèi)型[2]。很多研究表明DNA甲基化的變化與多種疾病的發(fā)病機(jī)制有關(guān)。例如，在許多癌癥類(lèi)型中，腫瘤細(xì)胞通常表現(xiàn)出與健康細(xì)胞不同的DNA甲基化模式，而腫瘤抑制基因很可能因啟動(dòng)子區(qū)域的高甲基化而失活[3]。CpG島（CpG island）的高甲基化與腎上腺皮質(zhì)癌（Adrenocortical carcinoma，ACC）癌癥的低存活率有關(guān)[4]。近年來(lái)高通量基因測(cè)序技術(shù)雖取得了巨大進(jìn)步，但DNA甲基化測(cè)序的序列數(shù)據(jù)因?yàn)楦鞣N技術(shù)局限性還常包含很大一部分缺失值導(dǎo)致無(wú)法直接使用[5]。因此有必要對(duì)這些序列數(shù)據(jù)中的缺失值進(jìn)行估算填補(bǔ)。

缺失數(shù)據(jù)可以分為三種類(lèi)型[6]：（1）完全隨機(jī)缺失（missing completely at random，MCAR），缺失情況與任何變量無(wú)關(guān)；（2）隨機(jī)缺失（missing at random，MAR），缺失的可能性不取決于缺失值，而是取決于觀察值；（3）非隨機(jī)缺失（missing not at random，MNAR），缺失情況只發(fā)生在特定的缺失數(shù)據(jù)上。填補(bǔ)方法需要對(duì)缺失數(shù)據(jù)的類(lèi)型進(jìn)行假設(shè)，沒(méi)有研究明確指出DNA甲基化的數(shù)據(jù)缺失類(lèi)型。根據(jù)已有文獻(xiàn)的假設(shè)和實(shí)驗(yàn)經(jīng)驗(yàn)，假設(shè)DNA甲基化的缺失類(lèi)型為MCAR或MAR。已有一些基于統(tǒng)計(jì)和機(jī)器學(xué)習(xí)的方法被應(yīng)用于隨機(jī)缺失值的填補(bǔ)。如最近鄰KNN（K-nearest neighbor）方法通過(guò)選擇與樣本最相似的K個(gè)樣本，計(jì)算這些樣本的平均值用于填補(bǔ)該樣本的缺失部分。郝勝軒等[7]基于KNN使用所有非噪聲最近鄰對(duì)缺失數(shù)據(jù)進(jìn)行填補(bǔ)取得很好效果。主成分分析（principal component analysis，PCA）[8]和奇異值分解（singular value decomposition，SVD）方法[9]基于相同原理都是通過(guò)計(jì)算數(shù)據(jù)矩陣的協(xié)方差矩陣，獲得特征值的特征向量，通過(guò)特征向量可以生成新的矩陣用于缺失數(shù)據(jù)的填補(bǔ)。

深度學(xué)習(xí)是基于數(shù)據(jù)表示的機(jī)器學(xué)習(xí)方法，通過(guò)組合低級(jí)特征形成更加抽象的高級(jí)表示特征，以發(fā)現(xiàn)數(shù)據(jù)的分布式特征[10]。它強(qiáng)大的建模復(fù)雜非線性關(guān)系的能力在許多領(lǐng)域已顯示出極大的優(yōu)越性[11]，如語(yǔ)音識(shí)別[12]、圖像識(shí)別[13]和自然語(yǔ)言處理[14]等領(lǐng)域。近年來(lái)，深度學(xué)習(xí)的發(fā)展已被越來(lái)越多的應(yīng)用于生物信息領(lǐng)域，如蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)[15-17]、藥物設(shè)計(jì)[18-20]、基因組分析[21-23]和單細(xì)胞基因組分析[24-27]。自編碼器（Auto-Encoders，AE）是深度網(wǎng)絡(luò)中一種常見(jiàn)的模型，具備生成與訓(xùn)練數(shù)據(jù)相似數(shù)據(jù)樣本的功能，可應(yīng)用于數(shù)據(jù)增廣和填補(bǔ)。已有研究者將轉(zhuǎn)移學(xué)習(xí)技術(shù)與自編碼器相結(jié)合來(lái)插補(bǔ)缺失的RNA序列數(shù)據(jù)[28]。變分自編碼器[29]（variational auto-encoder，VAE）是深度隱含空間生成模型的一種，通過(guò)學(xué)習(xí)隱含編碼空間與數(shù)據(jù)生成空間的特征映射，進(jìn)而在輸出端重構(gòu)生成輸入數(shù)據(jù)。它是無(wú)監(jiān)督特征學(xué)習(xí)的重要工具，在圖像生成方向取得了極大的成功[30]。正是VAE這種強(qiáng)大的重構(gòu)輸入數(shù)據(jù)的能力，為DNA甲基化缺失數(shù)據(jù)的填補(bǔ)提供了新的機(jī)遇。

1 填補(bǔ)實(shí)驗(yàn)材料

1.1 數(shù)據(jù)集

基準(zhǔn)實(shí)驗(yàn)數(shù)據(jù)集選自The Cancer Genome Atlas（TCGA）平臺(tái)[31]，包含信息有：DNA甲基化數(shù)據(jù)（JHU-USC Human Methylation 450類(lèi)型）、臨床信息和后續(xù)隨訪信息。數(shù)據(jù)集中共有33種不同的癌癥樣本9 756個(gè)，每個(gè)樣本具有485 577個(gè)甲基化位點(diǎn)。DNA甲基化水平通過(guò)計(jì)算甲基化和未甲基化等位基因之間的強(qiáng)度比來(lái)確定，稱(chēng)為β值（范圍從0到1）。數(shù)據(jù)樣本中的部分位點(diǎn)在所有樣本上都缺失，這樣的缺失位點(diǎn)沒(méi)有填補(bǔ)價(jià)值在填補(bǔ)之前先刪除。

TCGA平臺(tái)上不同癌癥樣本分布極為不均從幾十到幾百不等，考慮到深度學(xué)習(xí)需要的樣本量和癌癥的代表性，最終從平臺(tái)中挑選了發(fā)病率很高的肺癌（lung adenocarcinoma，LUAD）作為填補(bǔ)實(shí)驗(yàn)樣本，LUAD含有507個(gè)樣本。

1.2 模擬缺失值

由于癌癥樣本數(shù)據(jù)中缺失位置對(duì)應(yīng)的真實(shí)值無(wú)法獲取，則不能對(duì)填補(bǔ)值的準(zhǔn)確性進(jìn)行直觀評(píng)估。為了更好地評(píng)估不同填補(bǔ)方法的填補(bǔ)性能，需要在模擬缺失值上進(jìn)行仿真實(shí)驗(yàn)。主要流程為：（1）從LUAD數(shù)據(jù)集中選擇方差最大的前10 000個(gè)無(wú)缺失位點(diǎn)形成仿真數(shù)據(jù)集D∈R507×10000（507為樣本數(shù)量，10 000為每個(gè)樣本的甲基化位點(diǎn)數(shù)量）。（2）在仿真數(shù)據(jù)集中隨機(jī)引入不同的缺失比例（取20%、40%、60%、80%），并用相同大小的矩陣記錄下缺失位置信息。（3）對(duì)引入缺失值的數(shù)據(jù)進(jìn)行填補(bǔ)。圖1為在D∈R5×5矩陣的每列隨機(jī)引入20%缺失值的流程。原始矩陣大小為[5，5]，每列引入1個(gè)（5×0.2）缺失值，將該位置數(shù)值置為0。同時(shí)，在標(biāo)記矩陣中的相同位置置為1代表引入的缺失值位置。填補(bǔ)操作完成后，填補(bǔ)矩陣中X?ij為填補(bǔ)后的值。通過(guò)標(biāo)記矩陣中1的位置，在原始矩陣和填補(bǔ)后矩陣中找到填補(bǔ)前后的值就可進(jìn)行準(zhǔn)確性評(píng)估。

圖1 模擬缺失值引入流程Fig.1 Simulate process of introducing missing values

2 填補(bǔ)實(shí)驗(yàn)

填補(bǔ)實(shí)驗(yàn)主要分為三個(gè)部分：（1）基于變分自編碼器原理建立用于填補(bǔ)DNA甲基化缺失值的模型VAEMethImp。（2）將VAE-MethImp模型與現(xiàn)有的方法進(jìn)行對(duì)比，通過(guò)均方根誤差（root mean squared error，RMSE）和擬合度（R-squared，R2）對(duì)填補(bǔ)值的精度進(jìn)行評(píng)估。（3）通過(guò)生存分析實(shí)驗(yàn)來(lái)驗(yàn)證填補(bǔ)值的真實(shí)有效性。

2.1 變分自編碼器模型結(jié)構(gòu)

自編碼器模型主要有兩個(gè)部分組成：（1）編碼器（Encoder）：學(xué)習(xí)輸入數(shù)據(jù)的隱含特征的空間表示；（2）解碼器（Decoder）：從學(xué)習(xí)到的低維特征中重構(gòu)出原始的輸入數(shù)據(jù)。如圖2所示，編碼器h將原始輸入X編碼運(yùn)算得到特征X＇,X＇=h(X)；解碼器f將特征X＇解碼生成X?,即：X?=f(X＇)=f(h(X))。X和X?越接近意味著提取的X＇特征代表性越好。

圖2 自編碼器結(jié)構(gòu)示意圖Fig.2 Structure diagram of auto-encoder

變分自編碼器除了擁有與圖2中自編碼器結(jié)構(gòu)相似的編碼器和解碼器外，還新增了隱含編碼，可以看作是神經(jīng)網(wǎng)絡(luò)和貝葉斯網(wǎng)絡(luò)的混合體。在VAE中，隱藏編碼對(duì)應(yīng)的結(jié)點(diǎn)可以看成是隨機(jī)變量，其他結(jié)點(diǎn)視為普通神經(jīng)元。編碼器功能是一個(gè)變分推斷網(wǎng)絡(luò)，用來(lái)推薦均值和方差，而解碼器可以看作是將隱變量映射到觀測(cè)變量的生成網(wǎng)絡(luò)。VAE-MethImp模型如圖3所示，模型共分為3部分：

圖3 VAE-MethImp模型結(jié)構(gòu)圖Fig.3 VAE-MethImp model structure diagram

（1）編碼層：編碼層從輸入X推斷出分布的均值和方差，功能與自編碼器中的Encoder功能一致，Z_mean,Z_log_var=h(X)。

（2）隱含變量Z：Z是通過(guò)編碼層輸出的均值和方差計(jì)算出的輸入數(shù)據(jù)的專(zhuān)屬正態(tài)分布。

（3）解碼層：Z通過(guò)解碼層生成重構(gòu)后的數(shù)據(jù)X?。

編碼層用來(lái)推斷Z的均值和方差，再?gòu)恼龖B(tài)分布中取采樣ε，通過(guò)計(jì)算得到Z，如公式（1）所示：

其中zi∈Z，mi∈Z_mean,σi∈Z_log_var。

VAE-MethImp模型訓(xùn)練的目標(biāo)有兩個(gè)：編碼層需要將推斷出的正態(tài)分布與標(biāo)準(zhǔn)正態(tài)分布的KL散度KL(N(m,σ2)||N(0,I))作為額外的Loss，最小化損失函數(shù)如公式（2）所示：

其中xi∈X,x?i∈X?。

超參數(shù)選擇是訓(xùn)練深度神經(jīng)網(wǎng)絡(luò)模型中重要的環(huán)節(jié)，選擇一組最優(yōu)超參數(shù)可以提高學(xué)習(xí)的性能和效果。通過(guò)網(wǎng)格搜索對(duì)最常見(jiàn)的7種超參數(shù)進(jìn)行嘗試和優(yōu)化，最終選擇的超參數(shù)集如表1所示。

表1 測(cè)試的超參數(shù)集與選擇Table 1 Tested and selected hyperparameters

2.2 對(duì)比方法的實(shí)施

對(duì)現(xiàn)有的填補(bǔ)方法進(jìn)行了分析與比較，確定了4種常見(jiàn)的填補(bǔ)方法，分別是Mean、KNN、imputePCA和SVD。隨機(jī)森林和多重填補(bǔ)法（MI）由于處理[507，10 000]大小的矩陣所耗費(fèi)的時(shí)間（超48 h未結(jié)束）無(wú)法接受被舍棄。在實(shí)現(xiàn)時(shí)imputePCA直接調(diào)用missMDA包實(shí)現(xiàn)，SVD實(shí)現(xiàn)由fancyimpute包提供。所有方法均使用默認(rèn)參數(shù)實(shí)現(xiàn)，同時(shí)進(jìn)行了一些關(guān)鍵參數(shù)的優(yōu)化。具體說(shuō)來(lái)，在KNN方法中鄰居數(shù)K取值為30（從10、20、30和50中挑選），SVD的最大特征值設(shè)置為10（從5、10、20和30中挑選）。

2.3 生存分析

RMSE和R2只能用來(lái)評(píng)估填補(bǔ)值的平均準(zhǔn)確性，無(wú)法衡量填補(bǔ)值的實(shí)際有效性。填補(bǔ)值的實(shí)際有效性可以用嶺回歸正則化的COX比例風(fēng)險(xiǎn)模型（cox proportional-hazards model）來(lái)預(yù)測(cè)。該模型使用R軟件包中適用于高維數(shù)據(jù)擬合的glmnet包構(gòu)建[32]。通過(guò)一致性指數(shù)（C-Index）和p值（p-value）兩個(gè)指標(biāo)對(duì)預(yù)測(cè)結(jié)果進(jìn)行評(píng)估。C-Index是由Harrell教授提出用于計(jì)算生存分析中的COX模型預(yù)測(cè)值與真實(shí)之間的區(qū)分度[33]，C-Index值為0.5即為隨機(jī)猜測(cè)，數(shù)值越高則區(qū)分度越好。p值衡量將患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組的概率，值越低區(qū)分度越好。

3 結(jié)果與討論

統(tǒng)計(jì)了VAE-MethImp和現(xiàn)有的4種方法（Mean、KNN、imputePCA、SVD）在LUAD仿真缺失數(shù)據(jù)集上的填補(bǔ)實(shí)驗(yàn)結(jié)果。通過(guò)RMSE、R2和C-Index三個(gè)指標(biāo)對(duì)各方法填補(bǔ)后的數(shù)值精度和有效性進(jìn)行了評(píng)估。

3.1 填補(bǔ)精度評(píng)估

圖4中顯示了VAE-MethImp方法與其他方法的填補(bǔ)結(jié)果對(duì)比，結(jié)果顯示VAE-MethImp方法優(yōu)于其他方法。具體來(lái)看，Mean和KNN擁有較高誤差，KNN的RMSE隨著缺失率的增大而迅速增加。當(dāng)缺失率超過(guò)40%時(shí)，KNN的結(jié)果差于Mean，因?yàn)樵诟呷笔氏?，KNN通過(guò)僅存的部分信息獲取的鄰居之間相似度差。imputePCA與SVD原理相同，但imputePCA基于更少的特征向量導(dǎo)致RMSE比SVD高。VAE-MethImp方法對(duì)比SVD將誤差平均縮小了4.8%，缺失率越低VAEMethImp的優(yōu)勢(shì)越大，在80%缺失率時(shí)兩者接近原因是缺失數(shù)據(jù)太多可提取的特征變少。R2指標(biāo)的趨勢(shì)與RMSE相同，VAE-MethImp填補(bǔ)后的數(shù)值相關(guān)性最好。

圖4 5種填補(bǔ)方法的RMSE和R2結(jié)果對(duì)比Fig.4 Comparisons of five methods by RMSE and R2

除了填補(bǔ)性能更好外，VAE-MethImp方法在模型訓(xùn)練完后，對(duì)數(shù)據(jù)的填補(bǔ)花費(fèi)時(shí)間也是屬于最優(yōu)行列。在對(duì)比新樣本的填補(bǔ)時(shí)間時(shí)未將模型訓(xùn)練時(shí)間放入考慮是因?yàn)橛?xùn)練模型時(shí)間為一次性花費(fèi)時(shí)間。在Titan 1080 GPU上測(cè)試對(duì)樣本的填補(bǔ)時(shí)間與Mean方法接近，而KNN、imputePCA和SVD方法花費(fèi)的時(shí)間取決于樣本的大小，測(cè)試實(shí)驗(yàn)結(jié)果顯示這些方法花費(fèi)的時(shí)間是VAE-MethImp的10倍以上。

3.2 生存分析結(jié)果

為了驗(yàn)證填補(bǔ)數(shù)值的真實(shí)有效性，通過(guò)嶺回歸正則化Cox模型對(duì)不同填補(bǔ)方法的填補(bǔ)數(shù)值進(jìn)行了生存分析。如圖5所示，VAE-MethImp在4種缺失率下始終表現(xiàn)最佳。在缺失率為20%時(shí)，VAE-MethImp的C-Index為0.618，比Mean、KNN、imputePCA和SVD的結(jié)果分別提升10.5%、8%、4.7%和4%。很明顯，所有填補(bǔ)方法的C-Index值都會(huì)隨著缺失率的增加而減小，而填補(bǔ)方法之間的差異也會(huì)擴(kuò)大。在缺失率為80%時(shí)，VAEMethImp比SVD方法的優(yōu)勢(shì)有所減少，與圖4中的精度對(duì)比結(jié)果趨勢(shì)相同。

圖5 5種方法的C-Index結(jié)果對(duì)比Fig.5 Comparisons of five methods by C-Index

3.3 LUAD和BRCA真實(shí)缺失數(shù)據(jù)實(shí)驗(yàn)

在模擬缺失數(shù)據(jù)集上進(jìn)行的仿真實(shí)驗(yàn)表明，VAEMethImp模型在RMSE、R2和C-Index指標(biāo)上都優(yōu)于其他方法。為了進(jìn)一步證明模型的魯棒性，需要在真實(shí)缺失數(shù)據(jù)集上進(jìn)行驗(yàn)證。除了已有的LUAD癌癥，還增加了乳腺癌（Breast adenocarcinoma，BRCA）數(shù)據(jù)集進(jìn)行真實(shí)填補(bǔ)實(shí)驗(yàn)，實(shí)驗(yàn)中包含的數(shù)據(jù)集如表2所示。首先，從LUAD和BRCA數(shù)據(jù)集中挑選缺失率最高的前10 000個(gè)位點(diǎn)，組成真實(shí)缺失數(shù)據(jù)集。然后，將5種方法分別對(duì)真實(shí)數(shù)據(jù)集進(jìn)行填補(bǔ)，由于真實(shí)缺失位點(diǎn)對(duì)應(yīng)的真實(shí)值無(wú)法獲得，故RMSE和R2指標(biāo)不能使用。最后，通過(guò)C-Index和生存曲線來(lái)查看不同方法的填補(bǔ)性能。圖6所示為5種填補(bǔ)方法填補(bǔ)后的數(shù)據(jù)集和原始未填補(bǔ)數(shù)據(jù)集（No impute）的C-Index。圖中顯示，未填補(bǔ)數(shù)據(jù)集的C-Index接近0.5，表明原始數(shù)據(jù)集由于缺失數(shù)據(jù)太多變成無(wú)效數(shù)據(jù)。而經(jīng)過(guò)不同方法填補(bǔ)后的C-Index有了提升，證明填補(bǔ)后的數(shù)據(jù)是有價(jià)值的，其中VAE-MethImp比SVD有了平均2%的提升，比未補(bǔ)齊的數(shù)據(jù)提升了12%。

表2 實(shí)驗(yàn)數(shù)據(jù)集信息Table 2 Experimental dataset information

圖6 不同填補(bǔ)方法在LUAD和BRCA癌癥的生存分析結(jié)果Fig.6 C-Index value of survival prediction of original data（no imputation）and imputed methylation values on LUAD and BRCA

圖7展示了經(jīng)過(guò)VAE-MethImp填補(bǔ)后的數(shù)據(jù)與未填補(bǔ)數(shù)據(jù)的生存曲線對(duì)比，填補(bǔ)之后的數(shù)據(jù)集對(duì)高低風(fēng)險(xiǎn)組的區(qū)分比未填補(bǔ)的高低風(fēng)險(xiǎn)組的區(qū)分曲線間隔大，區(qū)分度更明顯。兩個(gè)癌癥上的真實(shí)填補(bǔ)實(shí)驗(yàn)可以證明VAE-MethImp在填補(bǔ)性能上比其他方法更加優(yōu)秀，填補(bǔ)之后的數(shù)據(jù)包含更多有價(jià)值信息。

圖7 LUAD和BRCA的填補(bǔ)前后生存曲線對(duì)比圖Fig.7 Kaplan-Meier plots of predicted high-and low-risk patients by original and predicted values by VAE-MethImp for LUAD and BRCA

4 總結(jié)

在本項(xiàng)研究中，開(kāi)發(fā)了一種新的基于變分自編碼器的用于DNA甲基化缺失數(shù)據(jù)的填補(bǔ)方法，即VAEMethImp。利用VAE極強(qiáng)的重構(gòu)生成輸入數(shù)據(jù)的能力，有效地解決了DNA甲基化缺失的問(wèn)題。通過(guò)在模擬的缺失數(shù)據(jù)上的仿真實(shí)驗(yàn)，證明VAE-MethImp在RMSE和R2方面均優(yōu)于現(xiàn)有方法。通過(guò)對(duì)填補(bǔ)數(shù)據(jù)的生存分析進(jìn)一步證實(shí)了填補(bǔ)數(shù)值的有效性。在LUAD和BRCA癌癥的真實(shí)缺失數(shù)據(jù)集上的填補(bǔ)實(shí)驗(yàn)也表明，VAEMethImp方法比其他方法填補(bǔ)的數(shù)據(jù)，擁有更高的CIndex值和區(qū)分度更明晰的生存曲線。

目前，本研究?jī)H限于癌癥的DNA甲基化缺失數(shù)據(jù)的填補(bǔ)上，且只利用了單個(gè)癌癥樣本之間的相關(guān)性。后續(xù)可以考慮利用泛癌樣本間DNA甲基化的相關(guān)性來(lái)提高單個(gè)癌癥的數(shù)據(jù)填補(bǔ)精度。隨著測(cè)序技術(shù)的發(fā)展，該填補(bǔ)方法可以進(jìn)一步應(yīng)用于其他生物組學(xué)類(lèi)型、不同疾病類(lèi)型及其他任務(wù)上，例如年齡預(yù)測(cè)和細(xì)胞分類(lèi)。