基于VAE 的編碼DNA 載體阻斷事件聚類分析與研究

魏梓軒， 周家樂

（華東理工大學信息科學與工程學院，上海 200237）

納米孔道檢測技術作為一種可用于單分子檢測的有力工具，具有巨大的研究價值和應用潛力[1-4]。目前，納米孔道已經被成功應用于DNA 測序[5]、蛋白質檢測[6]、蛋白質折疊[7]等研究中。納米孔道檢測裝置通常由納米孔道及其連通的兩部分電解質溶液組成，在外加電場的驅動下，待測物分子穿越納米孔道時，使得流經納米孔道的離子電流發生變化，進而形成特征阻斷電流[8]。因此，待測物分子的穿越行為對離子電流產生的信號進行了調制，這些信號中包含了分子電性、尺寸和結構等特異性信息[3,9-10]。

待測物分子穿越納米孔道產生的瞬時電流阻斷為一次阻斷事件。為識別這些事件，可以利用電流閾值的方法來實現[11]。然而，由于納米孔道信號的信噪比較低，增加了后續特征提取和分析的困難。傳統的數據分析方法主要提取事件持續時間和阻斷幅值兩個特征，通過研究其散點分布對待測物分子進行區分[12-13]。該方法往往依賴于納米孔道對不同分析物的分辨率，具有一定的隨機性。即使利用具有較高分子空間結構分辨率的納米孔道，也難以完全對混合待測物散點分布中的每一個事件進行區分，因此，該數據分析方法對納米孔道阻斷事件的分辨率較低[14]。

文獻[15]利用DNA 分子鏈作為蛋白質分子的載體，并在載體DNA 鏈上的不同位置加以啞鈴型束夾修飾。鑲嵌啞鈴型束夾后，載體DNA 分子鑲嵌部位的體積增大，在穿越孔道過程中會產生二級阻斷事件。通過在預先設定的若干位置進行束夾鑲嵌，形成編碼結構，進一步通過分析二級阻斷事件特征，便可對編碼DNA 載體上的蛋白質分子進行分類。為實現對編碼DNA 載體的納米孔道信號進行分析，需精確識別二級阻斷事件的特征。然而，由于納米孔道檢測存在的局限，如幅值變化、信噪比低、信號持續時間的分布范圍大、重疊事件、DNA 分子鏈的纏繞和折疊等問題，使得自動化的分析方法難以適用于各種復雜的情況，需要依賴于手動調整參數[4,16-18]。Misiunas 等[10]利用卷積神經網絡的自動特征提取特性提出了QuipuNet 網絡結構，對編碼DNA 載體種類和有無蛋白質綁定進行分類，取得了較高精度。然而，QuipuNet 作為一種監督分類方法，依賴于人工對阻斷事件進行精確標注以用于模型的訓練，這需要耗費大量的時間，效率低。無監督學習方法不依賴于對阻斷事件進行標注，可自動發現信號中的類簇[19-20]。基于這一優勢，本文建立了基于深度嵌入特征[21]和變分自編碼器[22]的聚類模型，該模型實現了對不同持續時間的阻斷事件的低維嵌入，并在嵌入特征上進行聚類。為驗證此模型的聚類效果，將其與經典的聚類方法進行了對比分析。在無監督聚類精度的評價標準下，本文模型實現了對納米孔道數據的精確分析。

1 編碼DNA 載體阻斷事件數據集描述

Misiunas 等[10]將編碼DNA 載體的阻斷事件的表格信息整理到HDF5 文件格式中，可用Pandas 提供的數據接口進行讀取。該數據集提供了32 個納米孔道產生的58 178 條數據，包含了“000”～“111”8 種編碼DNA 載體有無蛋白質綁定的阻斷事件信息。每條數據包含編碼、納米孔道編號、阻斷事件信號數據點、有無蛋白質綁定、展開幅值等信息。本文采用該數據集，使用阻斷事件信號數據點作為模型輸入，輸出得到對應信號的8 種編碼聚類簇，并利用已有的標注信息對算法的效果進行分析。在對信號進行聚類前，還需對阻斷事件進行預處理，利用載體DNA 分子鏈的展開幅值（Unfolded level）去除不同納米孔道對信號幅值產生的影響，以避免模型對納米孔道編號過擬合。

2 深度嵌入聚類與變分自編碼器

2.1 深度嵌入聚類

2.2 變分自編碼器

3 無監督聚類框架

3.1 基于變分自編碼器的深度嵌入聚類

圖1 基于變分自編碼器的深度嵌入聚類Fig. 1 Deep embedded clustering based on the variational auto encoder

3.2 變分自編碼器的網絡結構設計

將卷積網絡結構設計用于聚類框架的變分自編碼器。如圖2 所示，該網絡輸入維度為700，卷積層輸出按兩組相同維度的特征圖為一個元胞進行鏈接，其維度按照自編碼器的通用設計方法采用逐級降低方式進行設定。因此，設定卷積層輸出的維度為{350, 350, 128, 128, 87, 87}，那么卷積核的大小依次為{7, 7, 5, 5, 3, 3}。對于越深層的卷積層，選擇更多數目的特征圖來提取信號中的特征，其數目依次為{16, 16, 32, 32, 64, 64}。將卷積層輸出連接至維度為1 024 的全連接層，并分別連接至均值μ 和標準差σ，其維度為32。再通過式(8)所示的重參數技巧，利用均值μ 和標準差σ 進行高斯分布的采樣，得到隱變量表示z，至此構建成編碼器結構。解碼器按照與編碼器鏡像的方式進行構建，卷積算子利用反卷積算子進行代替，以還原出信號空間的維度。該網絡中的激活函數全部采用LeakyRelu 函數：

其中：x 為各層的線性輸出；γ 為待優化參數。

4 仿真及結果分析

4.1 數據預處理

納米孔道檢測技術存在的局限使得編碼DNA 載體阻斷事件中包含著與具體編碼類別無關的干擾，因此，在應用本文的聚類方法前，有必要對數據進行預處理，以降低聚類模型對噪聲特征的過擬合風險。由于納米孔道的直徑具有固定尺寸，編碼DNA 載體自身直徑和加之以靶位點后的直徑之間存在較大的差距，因此會在孔道中產生不同的占位作用，從而產生不同大小的二級阻斷電流幅值。圖3示出了8 種不同編碼DNA 載體對應的阻斷電流信號，圖中紅色虛線表示載體DNA 分子鏈的展開幅值，黑色虛線表示阻斷事件的開孔電流幅值，箭頭表示開始位點至結束位點的方向。

根據圖3 中箭頭指向可以判斷，編碼DNA 載體兩端進入納米孔道的順序是不固定的，但通過判斷靶位點尖峰的偏向可以判斷載體兩端的進入順序，進而從高位到低位讀出編碼。同一種類、不同批次的納米孔道，由于制備條件的差異，會在尺寸和結構上存在偏差，尤其是固體納米孔道。如圖3 中紅色虛線表示了編碼DNA 載體進入納米孔道時去折疊過程中的主要電流水平，是對電流數據設置3 個聚類中心進行聚類分析的結果。在將開孔電流對齊后，可看到紅色虛線所示的電流水平存在明顯差異，因此，需要對這些信號的去折疊電流水平進行對齊，以防止模型聚焦于孔差異而帶來錯誤的聚類結果。

圖2 聚類框架中的變分編碼器結構Fig. 2 Structure of the encoder of variational auto encoder in the clustering framework

另外，從圖3 中還可以發現原始阻斷事件的持續時間并不一致，這是因為分子穿過納米孔道過程中受到初始速度、分子結構、進入角度等隨機因素影響。傳統分析方法通過統計分布來分析待測堿基鏈長度[12,26]，而本文將阻斷事件的信號輸入到變分自編碼器中。由于變分自編碼器為一個神經網絡模型，輸入信號維度固定，因此，除了將信號的電流水平對齊，還需要對信號進行長度補全。選取每個信號的前50 個數據點（開孔電流信號）分別計算其標準差，得到平均值為0.009 5。進而，選取均值μ = 0、標準差σ = 0.009 5 的高斯白噪聲，對信號的末尾補全至總共含700 個數據點。最終，不同編碼DNA 載體的阻斷事件數量如表1 所示。

圖3 8 種編碼DNA 載體的阻斷電流信號Fig. 3 Blockade signals of eight encoded DNA carriers

表1 編碼DNA 載體的阻斷事件數量Table 1 Numbers of blockade events produced by encoded DNA carriers

4.2 仿真結果及分析

仿真平臺配置：Intel?Xeon?CPU E5-2650 v4 @2.20 GHz, 252 GB RAM, 64-bit GNU/Linux Centos 7,Python 3.6.8, Keras 2.2.4, Tensorflow-mkl 1.12.0。本文選擇結合動量和自適應學習率的優化方法(Adam)，對式(12)所示的損失函數進行優化。訓練時，所選取的批大小為512，最大迭代次數為20 000，目標分布的更新間隔為140 次，即T = 140。

利用2.3 節中的優化方法和參數，對基于變分自編碼器深度聚類框架進行優化。隨機選取了4 個聚類中心，并分別繪制出離聚類中心最近的5 個阻斷事件樣本，如圖4 所示。結合表1 中的數據樣本數量可以發現，離聚類中心距離較近的樣本中出現頻率較高的“011”、“101”等編碼，同時也是樣本數據集中占比較高的編碼。表1 中的數據存在嚴重不平衡的現象，表明聚類過程中占比較高的樣本會完全淹沒占比較小的樣本，從而產生與編碼無關的聚類模型。

對表1 中數據進行隨機降采樣，使得8 種編碼的分布保持一致，最終得到每個編碼對應樣本數量為963 的數據集。在降采樣后的數據集上，利用2.3 節的方法對上述模型重新訓練。選取4 類聚類中心，并分別列出5 個最相似樣本，如圖5 所示。對比圖5 中的結果可以發現，降采樣后的聚類效果得到提升，每類聚類中心周圍的樣本純度較高。因此，對原始的數據樣本進行降采樣后，訓練得到的模型在聚類效果上得到了提升。

利用K-means、Auto-encoder + K-means (AE + Kmeans)、VAE + K-means、IDEC[24]、基于變分自編碼器的IDEC（VAE + IDEC）分別對降采樣后的數據進行聚類，利用式(13)對聚類后的結果進行評價，結果如表2 所示。不難發現，可聚類的變分自編碼器模型相較于其他聚類算法可顯著提升分析精度，提升量最高可達29%。K-means 聚類算法的結果較差，原因在于該算法要求原始數據嚴格對齊，更適用于表格型數據，并且聚類效果也容易受到噪聲的干擾。而編碼DNA 載體阻斷事件信號具有時序數據中變形、位移等特征，這對K-means 中計算樣本與聚類中心距離的方法提出了更高的要求。因此，變分自編碼器能更有效地提取出原始信號中與阻斷事件類型相關的特征，進而通過將變分自編碼器與聚類層聯合訓練的方式，使得編碼器產生的隱變量分布逐漸向聚類中心對齊。AE + K-means 等兩階段的聚類方法中，表示特征轉換的編碼層網絡的訓練與聚類過程是分離的，因此難以保證編碼器所轉換特征的對齊效果。在計算樣本與聚類中心距離時，引入了一定的噪聲，使得聚類效果較差。基于VAE 的聚類模型相比于基于AE 的模型都取得了較高的聚類精度，原因在于VAE 對數據的概率分布參數進行建模，而不僅僅是對輸入數據進行壓縮表示。概率分布參數反映了輸入數據分布的本質特征，基于此可以得到更優的聚類精度。

圖4 0 ～ 3 類中離聚類中心最近的5 個樣本Fig. 4 Five samples closest to the centroids of cluster 0 to 3

圖5 降采樣后的聚類結果Fig. 5 Clustering results after under-sampling

表2 K-means、AE + K-means、VAE + K-means、IDEC、VAE +IDEC 的聚類結果比較（8 個聚類中心）Table 2 Clustering results comparison of K-means, AE + Kmeans, VAE + K-means, IDEC, VAE + IDEC (8 clusters)

由于原始數據，即編碼DNA 載體的阻斷電流信號，存在不規則的信號特征，如持續時間分布的方差、極差較大，使得原本分布在時間域上的數據又可能分離出更多子類。因此，通過對表2 中各個聚類方法的比較，本文認為變分自編碼的隱變量特征最大化地保留了信號特征，通過對隱變量進行分布假設，使得原始信號與類別相關的特征被轉換到更低維的特征分布中。

從圖3 可以發現，相同編碼的DNA 載體進入納米孔道的方向是不確定的。因此，編碼001 和編碼100、編碼011 和編碼110 的DNA 載體穿越納米孔道會產生相同時序的阻斷事件。在該聚類問題中，無論是傳統聚類算法還是本文提出的基于變分自編碼器的聚類算法，都無法直接利用尖峰的偏向信息來建立較優的聚類模型,所以本文將編碼001 和編碼100、編碼011 和編碼110 的數據分別設定為同一類別，得到6 個中心的聚類模型。聚類結果如表3 所示。

將上述結果與表2 的結果進行對比，可以發現所有聚類算法的無監督聚類精度都得到了一定的提升，總體上升了12%。同時，本文提出的基于變分自編碼的深度嵌入聚類模型得到了較優的聚類精度，最高提升了23%。這說明基于變分自編碼器隱變量模型對于信息感知和壓縮的有效性。本文設計的變分自編碼器模型采用了卷積神經網絡的隱層結構，可見卷積運算在納米孔道信號濾波或特征提取中的實用性。

表3 K-means、AE + K-means、VAE + K-means、IDEC、VAE +IDEC 的聚類結果比較（6 個聚類中心） Table 3 Clustering results comparison of K-means, AE + Kmeans, VAE + K-means, IDEC, VAE + IDEC (6 clusters)

5 結 論

數字編碼的DNA 載體為檢測蛋白質分子提供了重要的途徑，然而傳統的數據分析方法往往依賴于手動進行大量的調參。本文提出了一種基于變分自編碼器的深度嵌入聚類方法。首先，利用阻斷事件數據集對變分自編碼器進行預訓練，對信號的隱變量空間進行建模。再將預訓練好的變分自編碼器加入到深度嵌入聚類的框架中，在隱變量空間中對信號的嵌入特征進行聚類。隱變量模型能夠保留信號中的主要特征，可以有效過濾阻斷事件信號的噪聲、變形等干擾。與傳統的以及多階段的聚類算法相比較，本文提出的聚類框架以及設計的變分自編碼器網絡實現了較高精度的聚類效果。

本文的無監督聚類方法與卷積神經網絡模型在DNA 載體阻斷事件中的高精度監督分類效果，展示了卷積運算在納米孔道數據處理中的研究前景。該方法不依賴于阻斷事件數據集的先驗信息，可自主地發現數據集中的聚類簇，對于識別未知的信號種類或編碼特征具有較好的研究和應用價值。