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

魏梓軒， 周家樂

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

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

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

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

1 編碼DNA 載體阻斷事件數(shù)據(jù)集描述

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

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

2.1 深度嵌入聚類

2.2 變分自編碼器

3 無監(jiān)督聚類框架

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

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

3.2 變分自編碼器的網(wǎng)絡(luò)結(jié)構(gòu)設(shè)計(jì)

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

其中：x 為各層的線性輸出；γ 為待優(yōu)化參數(shù)。

4 仿真及結(jié)果分析

4.1 數(shù)據(jù)預(yù)處理

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

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

圖2 聚類框架中的變分編碼器結(jié)構(gòu)Fig. 2 Structure of the encoder of variational auto encoder in the clustering framework

另外，從圖3 中還可以發(fā)現(xiàn)原始阻斷事件的持續(xù)時(shí)間并不一致，這是因?yàn)榉肿哟┻^納米孔道過程中受到初始速度、分子結(jié)構(gòu)、進(jìn)入角度等隨機(jī)因素影響。傳統(tǒng)分析方法通過統(tǒng)計(jì)分布來分析待測堿基鏈長度[12,26]，而本文將阻斷事件的信號(hào)輸入到變分自編碼器中。由于變分自編碼器為一個(gè)神經(jīng)網(wǎng)絡(luò)模型，輸入信號(hào)維度固定，因此，除了將信號(hào)的電流水平對(duì)齊，還需要對(duì)信號(hào)進(jìn)行長度補(bǔ)全。選取每個(gè)信號(hào)的前50 個(gè)數(shù)據(jù)點(diǎn)（開孔電流信號(hào)）分別計(jì)算其標(biāo)準(zhǔn)差，得到平均值為0.009 5。進(jìn)而，選取均值μ = 0、標(biāo)準(zhǔn)差σ = 0.009 5 的高斯白噪聲，對(duì)信號(hào)的末尾補(bǔ)全至總共含700 個(gè)數(shù)據(jù)點(diǎn)。最終，不同編碼DNA 載體的阻斷事件數(shù)量如表1 所示。

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

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

4.2 仿真結(jié)果及分析

仿真平臺(tái)配置：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。本文選擇結(jié)合動(dòng)量和自適應(yīng)學(xué)習(xí)率的優(yōu)化方法(Adam)，對(duì)式(12)所示的損失函數(shù)進(jìn)行優(yōu)化。訓(xùn)練時(shí)，所選取的批大小為512，最大迭代次數(shù)為20 000，目標(biāo)分布的更新間隔為140 次，即T = 140。

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

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

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

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

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

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

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

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

將上述結(jié)果與表2 的結(jié)果進(jìn)行對(duì)比，可以發(fā)現(xiàn)所有聚類算法的無監(jiān)督聚類精度都得到了一定的提升，總體上升了12%。同時(shí)，本文提出的基于變分自編碼的深度嵌入聚類模型得到了較優(yōu)的聚類精度，最高提升了23%。這說明基于變分自編碼器隱變量模型對(duì)于信息感知和壓縮的有效性。本文設(shè)計(jì)的變分自編碼器模型采用了卷積神經(jīng)網(wǎng)絡(luò)的隱層結(jié)構(gòu)，可見卷積運(yùn)算在納米孔道信號(hào)濾波或特征提取中的實(shí)用性。

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

5 結(jié) 論

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

本文的無監(jiān)督聚類方法與卷積神經(jīng)網(wǎng)絡(luò)模型在DNA 載體阻斷事件中的高精度監(jiān)督分類效果，展示了卷積運(yùn)算在納米孔道數(shù)據(jù)處理中的研究前景。該方法不依賴于阻斷事件數(shù)據(jù)集的先驗(yàn)信息，可自主地發(fā)現(xiàn)數(shù)據(jù)集中的聚類簇，對(duì)于識(shí)別未知的信號(hào)種類或編碼特征具有較好的研究和應(yīng)用價(jià)值。