變分自編碼器對甲基化缺失數據的填補

王新峰，黃 偉

1.吉首大學 軟件學院，湖南 吉首 416000

2.中山大學 計算機學院，廣州 510000

3.中南大學 計算機學院，長沙 410000

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

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

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

1 填補實驗材料

1.1 數據集

基準實驗數據集選自The Cancer Genome Atlas（TCGA）平臺[31]，包含信息有：DNA甲基化數據（JHU-USC Human Methylation 450類型）、臨床信息和后續隨訪信息。數據集中共有33種不同的癌癥樣本9 756個，每個樣本具有485 577個甲基化位點。DNA甲基化水平通過計算甲基化和未甲基化等位基因之間的強度比來確定，稱為β值（范圍從0到1）。數據樣本中的部分位點在所有樣本上都缺失，這樣的缺失位點沒有填補價值在填補之前先刪除。

TCGA平臺上不同癌癥樣本分布極為不均從幾十到幾百不等，考慮到深度學習需要的樣本量和癌癥的代表性，最終從平臺中挑選了發病率很高的肺癌（lung adenocarcinoma，LUAD）作為填補實驗樣本，LUAD含有507個樣本。

1.2 模擬缺失值

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

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

2 填補實驗

填補實驗主要分為三個部分：（1）基于變分自編碼器原理建立用于填補DNA甲基化缺失值的模型VAEMethImp。（2）將VAE-MethImp模型與現有的方法進行對比，通過均方根誤差（root mean squared error，RMSE）和擬合度（R-squared，R2）對填補值的精度進行評估。（3）通過生存分析實驗來驗證填補值的真實有效性。

2.1 變分自編碼器模型結構

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

圖2 自編碼器結構示意圖Fig.2 Structure diagram of auto-encoder

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

圖3 VAE-MethImp模型結構圖Fig.3 VAE-MethImp model structure diagram

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

（2）隱含變量Z：Z是通過編碼層輸出的均值和方差計算出的輸入數據的專屬正態分布。

（3）解碼層：Z通過解碼層生成重構后的數據X?。

編碼層用來推斷Z的均值和方差，再從正態分布中取采樣ε，通過計算得到Z，如公式（1）所示：

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

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

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

超參數選擇是訓練深度神經網絡模型中重要的環節，選擇一組最優超參數可以提高學習的性能和效果。通過網格搜索對最常見的7種超參數進行嘗試和優化，最終選擇的超參數集如表1所示。

表1 測試的超參數集與選擇Table 1 Tested and selected hyperparameters

2.2 對比方法的實施

對現有的填補方法進行了分析與比較，確定了4種常見的填補方法，分別是Mean、KNN、imputePCA和SVD。隨機森林和多重填補法（MI）由于處理[507，10 000]大小的矩陣所耗費的時間（超48 h未結束）無法接受被舍棄。在實現時imputePCA直接調用missMDA包實現，SVD實現由fancyimpute包提供。所有方法均使用默認參數實現，同時進行了一些關鍵參數的優化。具體說來，在KNN方法中鄰居數K取值為30（從10、20、30和50中挑選），SVD的最大特征值設置為10（從5、10、20和30中挑選）。

2.3 生存分析

RMSE和R2只能用來評估填補值的平均準確性，無法衡量填補值的實際有效性。填補值的實際有效性可以用嶺回歸正則化的COX比例風險模型（cox proportional-hazards model）來預測。該模型使用R軟件包中適用于高維數據擬合的glmnet包構建[32]。通過一致性指數（C-Index）和p值（p-value）兩個指標對預測結果進行評估。C-Index是由Harrell教授提出用于計算生存分析中的COX模型預測值與真實之間的區分度[33]，C-Index值為0.5即為隨機猜測，數值越高則區分度越好。p值衡量將患者分為高風險組和低風險組的概率，值越低區分度越好。

3 結果與討論

統計了VAE-MethImp和現有的4種方法（Mean、KNN、imputePCA、SVD）在LUAD仿真缺失數據集上的填補實驗結果。通過RMSE、R2和C-Index三個指標對各方法填補后的數值精度和有效性進行了評估。

3.1 填補精度評估

圖4中顯示了VAE-MethImp方法與其他方法的填補結果對比，結果顯示VAE-MethImp方法優于其他方法。具體來看，Mean和KNN擁有較高誤差，KNN的RMSE隨著缺失率的增大而迅速增加。當缺失率超過40%時，KNN的結果差于Mean，因為在高缺失率下，KNN通過僅存的部分信息獲取的鄰居之間相似度差。imputePCA與SVD原理相同，但imputePCA基于更少的特征向量導致RMSE比SVD高。VAE-MethImp方法對比SVD將誤差平均縮小了4.8%，缺失率越低VAEMethImp的優勢越大，在80%缺失率時兩者接近原因是缺失數據太多可提取的特征變少。R2指標的趨勢與RMSE相同，VAE-MethImp填補后的數值相關性最好。

圖4 5種填補方法的RMSE和R2結果對比Fig.4 Comparisons of five methods by RMSE and R2

除了填補性能更好外，VAE-MethImp方法在模型訓練完后，對數據的填補花費時間也是屬于最優行列。在對比新樣本的填補時間時未將模型訓練時間放入考慮是因為訓練模型時間為一次性花費時間。在Titan 1080 GPU上測試對樣本的填補時間與Mean方法接近，而KNN、imputePCA和SVD方法花費的時間取決于樣本的大小，測試實驗結果顯示這些方法花費的時間是VAE-MethImp的10倍以上。

3.2 生存分析結果

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

圖5 5種方法的C-Index結果對比Fig.5 Comparisons of five methods by C-Index

3.3 LUAD和BRCA真實缺失數據實驗

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

表2 實驗數據集信息Table 2 Experimental dataset information

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

圖7展示了經過VAE-MethImp填補后的數據與未填補數據的生存曲線對比，填補之后的數據集對高低風險組的區分比未填補的高低風險組的區分曲線間隔大，區分度更明顯。兩個癌癥上的真實填補實驗可以證明VAE-MethImp在填補性能上比其他方法更加優秀，填補之后的數據包含更多有價值信息。

圖7 LUAD和BRCA的填補前后生存曲線對比圖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 總結

在本項研究中，開發了一種新的基于變分自編碼器的用于DNA甲基化缺失數據的填補方法，即VAEMethImp。利用VAE極強的重構生成輸入數據的能力，有效地解決了DNA甲基化缺失的問題。通過在模擬的缺失數據上的仿真實驗，證明VAE-MethImp在RMSE和R2方面均優于現有方法。通過對填補數據的生存分析進一步證實了填補數值的有效性。在LUAD和BRCA癌癥的真實缺失數據集上的填補實驗也表明，VAEMethImp方法比其他方法填補的數據，擁有更高的CIndex值和區分度更明晰的生存曲線。

目前，本研究僅限于癌癥的DNA甲基化缺失數據的填補上，且只利用了單個癌癥樣本之間的相關性。后續可以考慮利用泛癌樣本間DNA甲基化的相關性來提高單個癌癥的數據填補精度。隨著測序技術的發展，該填補方法可以進一步應用于其他生物組學類型、不同疾病類型及其他任務上，例如年齡預測和細胞分類。