O6-甲基鳥嘌呤(O6-CMG)是DNA中的一種高致突變烷基化產物,它會導致生命體罹患胃腸道腫瘤。現有的研究主要是利用恥垢分枝桿菌膜蛋白(MspA)納米孔技術,借助枯草芽孢桿菌噬菌體Phi29 DNA多聚酶(Phi29 DNA polymerase)對突變進行精確定位。近年來,機器學習技術被廣泛應用于納米孔測序數據的分析,但是機器學習往往需要大量的數據標記,這給研究者們帶來了額外的工作負擔,大大影響了其實用性。因此,本文提出了一種納米無監督深度學習(nano-UDL)方法,該方法能自動識別含有突變段的納米孔數據。nano-UDL方法利用深度自動編碼器從納米孔數據中提取特征,然后通過均值漂移(MeanShift)聚類算法對特征數據進行分類。此外,該方法還聯合優化了聚類損失和重構損失,從而提取最優的特征用于聚類。實驗結果表明,nano-UDL方法在O6-CMG數據集上具有較高的識別精度,能準確識別出所有包含O6-CMG的序列段。為了進一步驗證nano-UDL方法的魯棒性,本文進行了超參數敏感性驗證和消融實驗。利用nano-UDL方法分析納米孔數據不但可以有效降低人工分析數據帶來的額外成本,而且對包括基因組測序在內的諸多生物研究具有重要意義。
引用本文: 關曉宇, 王宇, 張金月, 邵偉, 黃碩, 張道強. 基于無監督深度學習的納米孔測序O6-甲基鳥嘌呤識別. 生物醫學工程學雜志, 2022, 39(1): 139-148. doi: 10.7507/1001-5515.202104068 復制
引言
基因組測序有助于提高人類對疾病的認識,如何辨別與人類疾病相關的遺傳危險因素至關重要。O6-甲基鳥嘌呤(O6-carboxymethyl guanine,O6-CMG)在一定情況下會觸發突變,引起DNA復制時的編碼錯誤,這種錯誤往往與胃腸道腫瘤相關。近年來,納米孔測序技術已經成為了一種新興的大分子感知識別技術,可用于DNA測序[1-3]。在之前的研究工作中,Wang等[4]使用枯草芽孢桿菌噬菌體Phi29 DNA多聚酶 (Phi29 DNA polymerase)輔助的恥垢分枝桿菌膜蛋白(mycobacterium smegmatis porin A,MspA)納米孔作為識別O6-CMG的測序工具。納米孔裝置由兩個充滿液體的儲層組成,由納米級的孔道相連接[5]。儲層上下兩側的帶電粒子產生正負電位差,帶電粒子由于電位差的驅動而通過中間的孔道,因而產生了可以反映分子特性的電流信號[6-10]。
由于產生的這種電流信號主要是時間序列,這正是機器學習方法擅于處理的數據形式,因而最近越來越多的研究將機器學習方法用于處理納米孔數據。例如,Carral等[11]提出了一個無監督機器學習模型,可以通過分析基于二維納米孔原始數據來識別不同的分子事件。Farshad等[12]設計了一個具有列文伯格-馬夸爾特(Levenberg-Marquardt,LM)傳遞函數的兩層神經網絡來分析氮化硅(Si3N4)納米孔數據。對于納米孔甲基化數據, Jia等[13]訓練了自適應增強(adaptive boosting,AdaBoost)分類器,從牛津納米孔技術(oxford nanopore technology,ONT)原始數據中檢測DNA甲基化事件。Schreiber等[14]提出了一個隱馬爾可夫模型用來分割和整合納米孔數據,實驗結果表明,與以往的人工操作方法相比,該方法能顯著降低錯誤率。此外,Liu等[15]訓練了支持向量機(support vector machine,SVM)分類器可以實現高精度檢測納米孔N6-甲基腺嘌呤(N6-methyladenosine,m6A)RNA修飾。Ni等[16]提出了深度信號(DeepSignal)方法,該方法可以檢測基因組位點的5-甲基胞嘧啶(5-methylcytosine,5mC)和m6A甲基化信號,在基因組水平的精度(accuracy,ACC)為0.9。Simpson等[10]訓練了一個隱馬爾可夫模型來區分5mC和非5mC。Stoiber等[17]提出了一種新的算法納米源(Nanoraw),該算法可以有效識別4-甲基胞嘧啶(4-methylcytosine,4mC)、5mC和m6A三種不同的甲基化標記。
但是,上述大多數研究都是基于監督學習框架展開[18-20],而監督學習需要大量精確標注的樣本用于訓練模型,這就需要領域內的專家耗費大量的時間和精力對樣本進行精確標注,因此將監督學習應用到納米孔領域成本很高[21]。為此,Xu等[22]提出了一個方案,利用無監督聚類方法可以找出納米孔數據中的潛在關系,由于其學習過程是不需要任何標簽的,這將大大減少專家的工作量。基于上述研究,本文提出了一種納米無監督深度學習方法 (nano-unsupervised-deep-learning,nano-UDL)用于納米孔數據分析,并在O6-CMG納米孔數據集上驗證了方法的有效性。本文提出的nano-UDL方法可以自動識別突變序列,且無需人工標注標簽。nano-UDL方法首先采用自動編碼器從原始納米孔數據中提取序列特征;然后,對提取的特征進行聚類;最后,實現了對O6-CMG的高精度識別。本文擬定的對比方法有K均值聚類(K-Means)[23]、基于密度的噪聲應用空間聚類(density-based spatial clustering of applications with noise,DBSCAN)[24]、均值漂移(MeanShift)[25]、凝聚層次聚類(Agglomerative)[26]、譜聚類(SpectralClustering)[27]等。期望通過與上述方法進行對比研究,能夠證實本文方法的有效性,為今后納米孔突變檢測研究奠定理論研究基礎,為更多類似的生物信息相關研究問題提供研究思路。
本文的主要創新思路如下:
(1) nano-UDL方法是一種無監督學習方法,使用深度自動編碼器提取特征[28-29],隨后使用MeanShift聚類算法對提取的特征進行分類。
(2)對于O6-CMG納米孔數據集,本文擬通過超參數敏感性驗證和消融實驗證實nano-UDL方法的穩定性,并通過對比實驗的結果驗證了nano-UDL方法可以達到提升精度的目的。
1 方法
nano-UDL方法用于精確定位O6-CMG納米孔數據流程示意圖如圖1所示,整體的流程分為以下三個階段:

(a)數據生成過程;(b)nano-UDL方法的構建過程;(c)對數據進行分類并精準定位O6-CMG位點
Figure1. Diagram of nano-UDL for accurate location of O6-CMG nanopore data(a)The process of the data generate; (b) The construction process of nano-UDL method; (c) The data were classified and the O6-CMG locus are accurately located
(1)首先,對納米孔采集的原始序列執行分割算法,切割出若干子序列,如圖1(a)所示;
(2)其次,nano-UDL方法模型的構建階段,包含了自動編碼器提取特征用于聚類模型MeanShift聚類,然后聯合優化聚類損失與重構損失以實現自動編碼器參數調優,如圖1(b)所示;
(3)最后,將調優后的nano-UDL模型用于O6-CMG數據集聚類,分別形成了包含O6-CMG位點的序列簇(定義為類別A)和不包含O6-CMG位點的序列簇(定義為類別B)。只需對類別A中的序列進行計算即可實現對O6-CMG位點的精準識別,如圖1(c)所示。
1.1 數據生成
首先需要做的是切割原始的納米孔數據(本文所用數據來自Wang等[4]的工作,得到了作者的授權,如需引用此數據請聯系原作者)。O6-CMG的原始數據是長序列,包含所有突變序列和非突變序列。測序波形往往開始于一個非常高的位置,這個位置所對應的電流稱為開孔電流。當樣品通過納米孔時,電流下降。當樣品完全通過時,電流上升到開孔電流位置。在此期間產生的電流就是一個有效序列。換句話說,每一個事件都可以通過開孔電流的幅值來區分,低于開孔電流的序列段就是有效序列。每個小時間序列段就反映了分子的具體特性。O6-CMG數據的具體的分割圖如圖2所示,圖2中橫坐標為時間,單位是s,縱坐標為電流,單位是pA,實體三角為O6-CMG位點。本文基于切割信號基線對長原始序列進行切割,生成1 010個子序列,從而形成了O6-CMG納米孔數據集,其中包括存在突變段序列55個(類別A),非突變等噪聲序列955個(類別B)。

1.2 nano-UDL方法的構建
nano-UDL方法的思路是自編碼器提取特征之后用聚類算法進行聚類。核心功能是聚類,聚類問題可以視為將n個點的{xi∈X}被聚成k個簇,并且每個簇被表示成為μj,j=1, , k。
非線性映射 應用于數據變換,而不是直接在數據空間X中進行聚類,其中θ是可學習的參數,Z是潛在特征空間。
通過同時學習特征空間Z中的k個聚類中心 和將數據點映射到Z的深度神經網絡(deep neural network,DNN)參數來對數據進行聚類。nano-UDL方法有兩個階段:
(1)對深度自動編碼器進行參數初始化,提取特征用于聚類模型聚類,給定初始參數為θ和。本文應用MeanShift聚類算法對O6-CMG數據進行聚類[25]。
(2)參數優化,其中聚類損失以最小化庫爾貝克·萊布勒(Kullback-Leibler,KL)散度逼近原始分布[30]。解碼器的重構損失用均方誤差(mean squared error,MSE)來度量。聯合優化聚類損失和重構損失實現對參數的優化。
1.2.1 自動編碼器提取特征
本文選擇堆疊式自動編碼器來初始化網絡參數θ和簇質心{μj},因為最近的研究表明,它們的表現要明顯優于其他自動編碼器[28-29, 31-33],已經普遍應用于諸多真實場景。
堆疊式自動編碼器網絡的運行是逐層進行的,每一層都經過訓練,可以在隨機損壞后重建上一層的輸出[29]。自動編碼器是兩層神經網絡,如式(1)所示:
![]() |
其中,x是模型輸入,y是模型輸出,f(·)是隨機映射[34],它隨機將一部分輸入的維度更改為0,g1是編碼器層的激活函數,g2是解碼器層的激活函數,θ = {W1,b1,W2,b2}是網絡參數。為了訓練模型,將損失函數設置為最小二乘損失 ,將其最小化以實現收斂。當訓練完一層后,將當前層輸出作為下一層的輸入,以訓練下一層。在本文中所有采用的編碼器/解碼器網絡中都有整流線性單位(rectified linear units,ReLU)[35]。進行逐層訓練之后,將所有編碼器層與所有解碼器層連接在一起,反向順序訓練以構建深度自動編碼器,以最大程度地減少重構損失來對其進行微調。為了初始化聚類中心,把數據輸入到初始化的自動編碼器中以獲取嵌入的初始數據點,并在特征空間Z中執行MeanShift聚類算法以獲得k個初始質心
。
1.2.2 參數優化
給定初始估計的非線性映射fθ和初始簇,為了提高聚類性能,首先通過t分布來計算嵌入點與聚類質心之間的軟分配。接下來更新非線性映射fθ,并通過最小化軟分配分布和目標分布之間的KL散度來微調聚類中心。所有過程都是迭代的,直到收斂到最小值。用t分布作為核心來度量嵌入點zi與聚類中心μj相似性,如式(2)所示:
![]() |
其中, 對應嵌入后的xi∈X ,α為學生t分布的自由度,qij可以解釋為將樣本i分配到j簇的概率(即軟分配),
是聚類形成簇的標號,
是
所對應的聚類中心。由于在無監督設置下對驗證集執行交叉驗證是不合理的,所以在所有實驗中設置α = 1。
本文所建立的模型是通過匹配目標分布的軟分配來訓練的。為此,將KL散度作為目標損失Lc,散度計算的對象是軟分配qij與輔助分布pij,如式(3)所示:
![]() |
其中,如何選擇目標分布P是決定方法性能的關鍵,目標就是讓分布更加靠近聚類中心,所以本文需要目標分布P盡量接近于分布Q,以期待最終樣本有著更好的聚類結果。由于qij是軟分配,使用軟概率目標更自然、更靈活。在實驗中,首先將qij提高到2次方,然后通過每一簇的概率歸一化來計算pij,如式(4)所示:
![]() |
其中, 是聚類形成簇的標號,
是對應的針對于
的軟分配,
是針對于qij軟聚類概率,同理
是針對于
軟聚類概率。
解碼器重構損失Lr用MSE來度量,如式(5)所示:
![]() |
其中,,xi是輸入,
和
分別為編碼器映射和解碼器映射。自動編碼器可以保留數據生成分布的局部結構[36-37]。本文將聚類損失和重構損失合并為nano-UDL模型的損失函數L如式(6)所示:
![]() |
其中,γ是調節兩個損失的超參數調節因子,在這種情況下,使用整合后的損失不會對嵌入空間造成破壞。在本文2.6節將討論γ的取值。
為了計算最優損失函數,本文使用小批量隨機梯度下降(stochastic gradient decent,SGD)來優化聚類損失和重構損失。特別地,有兩種參數需要優化:自動編碼器的權值和聚類中心。Lc相對于嵌入點zi和聚類中心 μj的梯度如式(7)~式(8)所示:
![]() |
![]() |
其中,K是最大的聚類簇數,n為樣本總數,則對于小批量樣本m和學習率λ,更新后的 μj如式(9)所示:
![]() |
解碼器的權重 由下式進行更新,如式(10)所示:
![]() |
因此,編碼器的權重W通過下式進行更新,如式(11)所示:
![]() |
1.3 O6-CMG位點的定位
上述策略通過自動編碼器和聚類算法的聯合調優,構建了nano-UDL模型。接下來則需要nano-UDL模型對O6-CMG數據集進行聚類以實現對包含突變位點序列的篩選。由于O6-CMG數據集包含了兩個類別(類別A和類別B),因此O6-CMG數據集被聚類成為兩個簇分別對應類別A和類別B,如圖1(c)所示。
nano-UDL方法初步實現了對O6-CMG子序列的分類,而為了精確定位子序列中O6-CMG的位置,進一步需要對分類好的類別A中的序列執行分段計算。取類別A所在簇的其中一個序列S,將其分成n段等長小段序列[s1, s2, , sn],對這n個小段計算幅值變化率,取幅值變化最快的段即為O6-CMG位點所在的位置。具體計算方法如式(12)所示:
![]() |
其中,grad(·)是計算幅值變化率函數,本文統一取n = 100,通過上式計算最大的梯度所在的段即為當前序列S的O6-CMG位點(loc)所在的位置。
通過以上三個步驟(數據生成、nano-UDL方法的構建、O6-CMG位點的定位)可以實現在不需要人工干預的情況下,完成對O6-CMG數據集中的O6-CMG突變位點的精準定位。
2 實驗
本文在O6-CMG數據集上評估nano-UDL方法的性能,并與其他傳統的聚類算法進行比較。在本節中,不僅通過四個評價指標來展示nano-UDL方法的性能,還進行超參數敏感度驗證和消融實驗,以此驗證模型的魯棒性。
2.1 評價指標
評價指標是評價算法性能的指標。由于本文主要探討的是無監督學習,因此本文給出了四種標準的無監督評價指標,用于與其他算法進行評估和比較。基于數據集的類型,聚類簇的數量設置為2,本文用無監督聚類精度來評估性能,即準確率(accuracy,ACC),其計算如式(13)所示:
![]() |
其中,yi 為真實標簽,ci為算法產生的聚類分配,m為聚類與標簽之間可能的映射,n為樣本總數,1{·} 是非0即1判斷器,如果映射后的樣本標簽和真實標簽一致則為1,否則為0。
歸一化互信息(normalized mutual information,NMI)(符號為:NMI)用以計算來自同一數據集的兩個標簽之間相似度的歸一化度量,定義如式(14)所示:
![]() |
其中,為真實標簽l與預測聚類標簽c之間的互信息,H為其信息熵。NMI的結果不因簇(類)的排列而改變,它們被歸一化到[0, 1],0表示不相關,1表示完全相關。
另外一種評價指標是調整后的蘭德指數(adjusted rand index,ARI)[38],它經常用于聚類驗證,因為它是兩個類別之間的一致性的度量。
卡帕(Kappa)系數是新的聚類性能度量指標[39]。Kappa系數是用于一致性檢驗的指標,也可以用來衡量分級效果。因為對于分類問題來說,所謂一致性就是模型預測結果是否與實際分類結果一致。Kappa系數的計算是基于混淆矩陣的,混淆矩陣的范圍在?1~1之間,通常大于0。
基于混淆矩陣的Kappa系數計算公式如式(15)所示:
![]() |
其中,K表示為Kappa系數;po為觀測到的符合率,即ACC;pe表示隨機匹配率,即對應于所有類別的“實際數量與預測數量的乘積”之和除以“樣本總數的平方”。
2.2 實驗細節
一般監督學習使用驗證集上的交叉驗證來確定超參數,而非監督學習則不需要設置驗證集,因此本文使用默認超參數,盡量避免調優。本文將O6-CMG數據集的網絡尺寸依次設置為256、256、512和4,所有層都是全連接。
在預訓練階段,從均值為0和標準差為0.01的高斯分布中選取隨機數作為初始值。通過1 000次迭代對每個層進行預訓練。為了對整個自動編碼器進行參數調優,網絡執行2 000次的迭代。全部層都進行預訓練和調優,批次大小設置為256,開始學習率設置為0.1,每200次迭代衰減10%。參數γ將在第2.6節中進行討論。通過設置上述參數,使重構損失達到較好的效果。
為了初始化聚類中心,執行Meanshift聚類算法 20次,然后選擇最佳的聚類結果。為了使KL散度最小,設置學習率為0.01,收斂閾值設置為0.001。對于所有的基線算法,本文對參數進行調整,并選擇效果最好的結果進行比較。
本文主要目標是能夠最大限度地將突變段序列從整個分割后的序列中找出來,同時盡量降低噪聲序列被判斷為突變序列的情況,所以在最終的評價指標中,選擇使用最基本的二分類指標ACC作為初步評價指標,進一步增加了包括NMI、ARI、Kappa三種應用于衡量聚類效果的度量指標。
2.3 實驗結果
本文對不同的算法的性能進行了定量和定性的評估,如圖3所示,其中前5種對比方法的輸入特征是通過人工手動獲取,提取的特征包括諸如均值、方差等統計信息特征。對比方法的輸入特征是相同的,而本文所提出的nano-UDL方法輸入的特征是通過如圖1所示的自動編碼器優化獲取。通過觀察發現,在NMI、ARI、Kappa、ACC等所有評價指標上,nano-UDL方法均明顯高于其他方法,并在O6-CMG數據集上表現出良好的聚類性能。這是因為nano-UDL方法使用多層自動編碼器作為特征提取器,能夠捕捉到非深度學習模型無法獲得的局部特征。此外,圖3中前5種方法都是非深度學習模型,nano-UDL方法是基于深度學習模型的。由此可以觀察到基于深度學習模型的方法明顯比非深度學習的方法能發揮更好的性能。

自動編碼器和不同聚類算法結合的實驗結果對比如圖4所示。圖4中前4個柱狀圖都是自動編碼器與對應聚類算法結合后的結果,結果表明,nano-UDL方法優于其他組合。本文選用的MeanShift聚類算法針對于O6-CMG數據集有著良好的聚類效果,其得益于O6-CMG數據在通過自動編碼器特征提取后的特征分布更適合這種基于密度類型的聚類方法。

2.4 超參數敏感性
為了研究nano-UDL方法在O6-CMG數據集上的參數敏感性,本文對不同參數下的ACC進行柱狀圖直觀對比,如圖5所示。結果表明,nano-UDL方法在大多數參數組合下都保持了良好的結果,并且相對穩定。進一步地,對其他幾種度量指標也同樣做了類似的實驗,實驗結果同圖5類似,同樣證實了nano-UDL穩定性。

2.5 消融實驗
本文方法所采用的損失函數是聚類損失和解碼器重構損失的線性加權和。因此本文主要是通過對聚類損失和解碼器重構損失進行的聯合優化,使nano-UDL方法的自動編碼器模型得到最優的聚類結果。為驗證重構損失添加之前和之后對實驗結果的影響,本文執行了一個消融實驗。其中在“無重構損失”的情況下,ACC為98.12%,聚類時間為1.735 s,相對而言在“有重構損失”的情況下ACC為98.32%,聚類時間為2.228 s。實驗結果表示,“有重構損失”比“無重構損失”發揮出了更好的性能。對于聚類時間而言,“有重構損失”聚類時間略長于“無重構損失”,原因是計算解碼器重構損失需要額外時間。總的來說,將聚類損失和重構損失相結合時nano-UDL方法性能最好。
2.6 參數γ的確定
為了解如式(6)所示的聚類損失參數γ對nano-UDL方法性能的影響,本文在 [0.01, 100] 范圍中采樣,選取13個值。優化器設置為隨機梯度下降法(stochastic gradient descent, SGD),動量β=0.9,學習速率λ設為0.1。結果如圖6所示,可以觀察到當γ值過小時,聚類損失項的積極作用被消除,而值過大則會導致潛在特征空間的扭曲從而影響性能,其中橫坐標值做了對數拉伸處理,以便橫坐標可以等比例顯示。當γ接近0.5時,聚類結果可以達到最優。因此,本文在所有的實驗中設置γ=0.5。

2.7 與監督方法對比
在本節中,將nano-UDL方法與其他監督學習算法進行比較,驗證的算法有分類回歸樹[40](classification and regression trees,CART)、SVM[41]、AdaBoost[42]、隨機森林(random forest,RF)[43]、K近鄰(k-nearest neighbor,KNN)和卷積神經網絡(convolutional neural networks,CNN)。這些監督算法均具有不同的監督率,為此,本文在常規監督學習算法上使用0.1~0.9比例的訓練集用于模型訓練。除此之外,本文利用CNN模型來驗證傳統深度學習模型的性能。實驗結果如圖7所示,隨著監督率的增加,模型的分類效果也在逐步提高。雖然在完全監督情況下,nano-UDL方法不如某些全監督算法。但是,在沒有監督或者監督信息很少的情況下,本文所提出的nano-UDL方法具有顯著的優勢。除CNN和nano-UDL方法這兩種深度學習模型外,其余模型輸入采用相同的特征,從圖7中可以看出CART、SVM、AdaBoost與KNN在全監督情況下不如nano-UDL方法。也就是說,手工提取的特征用上述這四種模型訓練并不能達到預期的效果,但例外的是RF模型在同樣的實驗設置下卻可以表現得較好。這也從另一個側面驗證了深度學習相比于手工提取特征技術,可以針對性地對特定任務進行特征調整。與傳統機器學習相比不同的是,在某些特征情況下,模型間可能存在較大的性能偏差。而本文提出的nano-UDL方法在無監督的背景條件下,也能達到部分有監督信息訓練算法的性能,這也證實了無監督深度學習在納米孔數據分析領域具有一定的應用價值。

3 討論
為了證明本文所提出框架的有效性,本文采用了消融實驗、精度對比等機器學習模型驗證策略對nano-UDL方法進行了驗證。通過對實驗輸出結果進行分析,本文從以下兩個方面對nano-UDL方法進行討論:
可靠性:nano-UDL方法有效地保證了突變序列檢測的高精度和高置信度,這是構建可靠且穩定系統的必然要求。如圖3所示的nano-UDL方法的結果比一些傳統的聚類算法具有更高精度和置信度。從圖5的結果可以看出超參數對聚類結果的影響很小,充分說明了本文提出模型的穩定性。通過改變監督樣本的大小,從圖7可以看出,nano-UDL方法的結果相對較好而且相對穩定。綜上所述,實驗結果充分說明了nano-UDL方法的可靠性。
訓練代價:雖然nano-UDL方法相比于其他聚類方法有著較好的性能,但由于這種方法策略是基于多梯度進行計算的,其中聚類損失和重構損失,增加了訓練的時間復雜度。在訓練階段,本文使用少量的樣本對模型進行預訓練,以減少模型訓練的成本。通過使用預先訓練的模型初始化網絡,可以有效地保證模型在訓練過程中時刻處于最佳。由于本文的訓練數據集樣本較少(1 010例),所以消耗的系統內存、GPU資源等計算資源都在有效范圍內。
目前的納米孔數據分析技術都是基于統計學習方法或有監督的機器學習算法,不能保證時間的有效性,因為其需要大量的分析計算成本。本文提出的方法是一種可以降低學習代價的無監督學習算法,該方法對納米孔DNA突變的研究有較好的效果,將對今后更多納米孔DNA突變的研究提供指導思路。本文提出的方法是嘗試將無監督深度學習方法應用到納米孔數據的首次探索,希望這種在生物信息學應用中利用先進的機器學習技術的思想在未來能夠逐漸走向成熟。
4 結論
針對O6-CMG納米孔數據,本文提出了一種新的無監督深度聚類模型nano-UDL,該模型將數據點聚類到一個可以通過聯合優化的特征空間中,最終實現精準識別突變位點的目的。為了逼近目標分布,該方法通過最小化聚類損失和重構損失進行迭代訓練。nano-UDL方法可以看作是一種不需要任何人工干預的啟發式工具,它可以用來處理原始納米孔數據。通過實驗分析,nano-UDL方法在提高性能的同時,對超參數設置具有魯棒性。nano-UDL這種無監督深度學習策略可以在減少人工計算成本的情況直接檢測出O6-CMG突變位點。期望本文提出的方法將對未來更多類似的從DNA測序中精準定位突變位點的任務提供一個理論依據。
重要聲明
利益沖突聲明:本文全體作者均聲明不存在利益沖突。
作者貢獻聲明:關曉宇主要負責項目主持、算法程序設計以及論文編寫;王宇和張金月主要負責實驗數據的采集和分析; 邵偉, 黃碩, 張道強主要負責提供實驗指導,數據分析指導,論文審閱修訂。
引言
基因組測序有助于提高人類對疾病的認識,如何辨別與人類疾病相關的遺傳危險因素至關重要。O6-甲基鳥嘌呤(O6-carboxymethyl guanine,O6-CMG)在一定情況下會觸發突變,引起DNA復制時的編碼錯誤,這種錯誤往往與胃腸道腫瘤相關。近年來,納米孔測序技術已經成為了一種新興的大分子感知識別技術,可用于DNA測序[1-3]。在之前的研究工作中,Wang等[4]使用枯草芽孢桿菌噬菌體Phi29 DNA多聚酶 (Phi29 DNA polymerase)輔助的恥垢分枝桿菌膜蛋白(mycobacterium smegmatis porin A,MspA)納米孔作為識別O6-CMG的測序工具。納米孔裝置由兩個充滿液體的儲層組成,由納米級的孔道相連接[5]。儲層上下兩側的帶電粒子產生正負電位差,帶電粒子由于電位差的驅動而通過中間的孔道,因而產生了可以反映分子特性的電流信號[6-10]。
由于產生的這種電流信號主要是時間序列,這正是機器學習方法擅于處理的數據形式,因而最近越來越多的研究將機器學習方法用于處理納米孔數據。例如,Carral等[11]提出了一個無監督機器學習模型,可以通過分析基于二維納米孔原始數據來識別不同的分子事件。Farshad等[12]設計了一個具有列文伯格-馬夸爾特(Levenberg-Marquardt,LM)傳遞函數的兩層神經網絡來分析氮化硅(Si3N4)納米孔數據。對于納米孔甲基化數據, Jia等[13]訓練了自適應增強(adaptive boosting,AdaBoost)分類器,從牛津納米孔技術(oxford nanopore technology,ONT)原始數據中檢測DNA甲基化事件。Schreiber等[14]提出了一個隱馬爾可夫模型用來分割和整合納米孔數據,實驗結果表明,與以往的人工操作方法相比,該方法能顯著降低錯誤率。此外,Liu等[15]訓練了支持向量機(support vector machine,SVM)分類器可以實現高精度檢測納米孔N6-甲基腺嘌呤(N6-methyladenosine,m6A)RNA修飾。Ni等[16]提出了深度信號(DeepSignal)方法,該方法可以檢測基因組位點的5-甲基胞嘧啶(5-methylcytosine,5mC)和m6A甲基化信號,在基因組水平的精度(accuracy,ACC)為0.9。Simpson等[10]訓練了一個隱馬爾可夫模型來區分5mC和非5mC。Stoiber等[17]提出了一種新的算法納米源(Nanoraw),該算法可以有效識別4-甲基胞嘧啶(4-methylcytosine,4mC)、5mC和m6A三種不同的甲基化標記。
但是,上述大多數研究都是基于監督學習框架展開[18-20],而監督學習需要大量精確標注的樣本用于訓練模型,這就需要領域內的專家耗費大量的時間和精力對樣本進行精確標注,因此將監督學習應用到納米孔領域成本很高[21]。為此,Xu等[22]提出了一個方案,利用無監督聚類方法可以找出納米孔數據中的潛在關系,由于其學習過程是不需要任何標簽的,這將大大減少專家的工作量。基于上述研究,本文提出了一種納米無監督深度學習方法 (nano-unsupervised-deep-learning,nano-UDL)用于納米孔數據分析,并在O6-CMG納米孔數據集上驗證了方法的有效性。本文提出的nano-UDL方法可以自動識別突變序列,且無需人工標注標簽。nano-UDL方法首先采用自動編碼器從原始納米孔數據中提取序列特征;然后,對提取的特征進行聚類;最后,實現了對O6-CMG的高精度識別。本文擬定的對比方法有K均值聚類(K-Means)[23]、基于密度的噪聲應用空間聚類(density-based spatial clustering of applications with noise,DBSCAN)[24]、均值漂移(MeanShift)[25]、凝聚層次聚類(Agglomerative)[26]、譜聚類(SpectralClustering)[27]等。期望通過與上述方法進行對比研究,能夠證實本文方法的有效性,為今后納米孔突變檢測研究奠定理論研究基礎,為更多類似的生物信息相關研究問題提供研究思路。
本文的主要創新思路如下:
(1) nano-UDL方法是一種無監督學習方法,使用深度自動編碼器提取特征[28-29],隨后使用MeanShift聚類算法對提取的特征進行分類。
(2)對于O6-CMG納米孔數據集,本文擬通過超參數敏感性驗證和消融實驗證實nano-UDL方法的穩定性,并通過對比實驗的結果驗證了nano-UDL方法可以達到提升精度的目的。
1 方法
nano-UDL方法用于精確定位O6-CMG納米孔數據流程示意圖如圖1所示,整體的流程分為以下三個階段:

(a)數據生成過程;(b)nano-UDL方法的構建過程;(c)對數據進行分類并精準定位O6-CMG位點
Figure1. Diagram of nano-UDL for accurate location of O6-CMG nanopore data(a)The process of the data generate; (b) The construction process of nano-UDL method; (c) The data were classified and the O6-CMG locus are accurately located
(1)首先,對納米孔采集的原始序列執行分割算法,切割出若干子序列,如圖1(a)所示;
(2)其次,nano-UDL方法模型的構建階段,包含了自動編碼器提取特征用于聚類模型MeanShift聚類,然后聯合優化聚類損失與重構損失以實現自動編碼器參數調優,如圖1(b)所示;
(3)最后,將調優后的nano-UDL模型用于O6-CMG數據集聚類,分別形成了包含O6-CMG位點的序列簇(定義為類別A)和不包含O6-CMG位點的序列簇(定義為類別B)。只需對類別A中的序列進行計算即可實現對O6-CMG位點的精準識別,如圖1(c)所示。
1.1 數據生成
首先需要做的是切割原始的納米孔數據(本文所用數據來自Wang等[4]的工作,得到了作者的授權,如需引用此數據請聯系原作者)。O6-CMG的原始數據是長序列,包含所有突變序列和非突變序列。測序波形往往開始于一個非常高的位置,這個位置所對應的電流稱為開孔電流。當樣品通過納米孔時,電流下降。當樣品完全通過時,電流上升到開孔電流位置。在此期間產生的電流就是一個有效序列。換句話說,每一個事件都可以通過開孔電流的幅值來區分,低于開孔電流的序列段就是有效序列。每個小時間序列段就反映了分子的具體特性。O6-CMG數據的具體的分割圖如圖2所示,圖2中橫坐標為時間,單位是s,縱坐標為電流,單位是pA,實體三角為O6-CMG位點。本文基于切割信號基線對長原始序列進行切割,生成1 010個子序列,從而形成了O6-CMG納米孔數據集,其中包括存在突變段序列55個(類別A),非突變等噪聲序列955個(類別B)。

1.2 nano-UDL方法的構建
nano-UDL方法的思路是自編碼器提取特征之后用聚類算法進行聚類。核心功能是聚類,聚類問題可以視為將n個點的{xi∈X}被聚成k個簇,并且每個簇被表示成為μj,j=1, , k。
非線性映射 應用于數據變換,而不是直接在數據空間X中進行聚類,其中θ是可學習的參數,Z是潛在特征空間。
通過同時學習特征空間Z中的k個聚類中心 和將數據點映射到Z的深度神經網絡(deep neural network,DNN)參數來對數據進行聚類。nano-UDL方法有兩個階段:
(1)對深度自動編碼器進行參數初始化,提取特征用于聚類模型聚類,給定初始參數為θ和。本文應用MeanShift聚類算法對O6-CMG數據進行聚類[25]。
(2)參數優化,其中聚類損失以最小化庫爾貝克·萊布勒(Kullback-Leibler,KL)散度逼近原始分布[30]。解碼器的重構損失用均方誤差(mean squared error,MSE)來度量。聯合優化聚類損失和重構損失實現對參數的優化。
1.2.1 自動編碼器提取特征
本文選擇堆疊式自動編碼器來初始化網絡參數θ和簇質心{μj},因為最近的研究表明,它們的表現要明顯優于其他自動編碼器[28-29, 31-33],已經普遍應用于諸多真實場景。
堆疊式自動編碼器網絡的運行是逐層進行的,每一層都經過訓練,可以在隨機損壞后重建上一層的輸出[29]。自動編碼器是兩層神經網絡,如式(1)所示:
![]() |
其中,x是模型輸入,y是模型輸出,f(·)是隨機映射[34],它隨機將一部分輸入的維度更改為0,g1是編碼器層的激活函數,g2是解碼器層的激活函數,θ = {W1,b1,W2,b2}是網絡參數。為了訓練模型,將損失函數設置為最小二乘損失 ,將其最小化以實現收斂。當訓練完一層后,將當前層輸出作為下一層的輸入,以訓練下一層。在本文中所有采用的編碼器/解碼器網絡中都有整流線性單位(rectified linear units,ReLU)[35]。進行逐層訓練之后,將所有編碼器層與所有解碼器層連接在一起,反向順序訓練以構建深度自動編碼器,以最大程度地減少重構損失來對其進行微調。為了初始化聚類中心,把數據輸入到初始化的自動編碼器中以獲取嵌入的初始數據點,并在特征空間Z中執行MeanShift聚類算法以獲得k個初始質心
。
1.2.2 參數優化
給定初始估計的非線性映射fθ和初始簇,為了提高聚類性能,首先通過t分布來計算嵌入點與聚類質心之間的軟分配。接下來更新非線性映射fθ,并通過最小化軟分配分布和目標分布之間的KL散度來微調聚類中心。所有過程都是迭代的,直到收斂到最小值。用t分布作為核心來度量嵌入點zi與聚類中心μj相似性,如式(2)所示:
![]() |
其中, 對應嵌入后的xi∈X ,α為學生t分布的自由度,qij可以解釋為將樣本i分配到j簇的概率(即軟分配),
是聚類形成簇的標號,
是
所對應的聚類中心。由于在無監督設置下對驗證集執行交叉驗證是不合理的,所以在所有實驗中設置α = 1。
本文所建立的模型是通過匹配目標分布的軟分配來訓練的。為此,將KL散度作為目標損失Lc,散度計算的對象是軟分配qij與輔助分布pij,如式(3)所示:
![]() |
其中,如何選擇目標分布P是決定方法性能的關鍵,目標就是讓分布更加靠近聚類中心,所以本文需要目標分布P盡量接近于分布Q,以期待最終樣本有著更好的聚類結果。由于qij是軟分配,使用軟概率目標更自然、更靈活。在實驗中,首先將qij提高到2次方,然后通過每一簇的概率歸一化來計算pij,如式(4)所示:
![]() |
其中, 是聚類形成簇的標號,
是對應的針對于
的軟分配,
是針對于qij軟聚類概率,同理
是針對于
軟聚類概率。
解碼器重構損失Lr用MSE來度量,如式(5)所示:
![]() |
其中,,xi是輸入,
和
分別為編碼器映射和解碼器映射。自動編碼器可以保留數據生成分布的局部結構[36-37]。本文將聚類損失和重構損失合并為nano-UDL模型的損失函數L如式(6)所示:
![]() |
其中,γ是調節兩個損失的超參數調節因子,在這種情況下,使用整合后的損失不會對嵌入空間造成破壞。在本文2.6節將討論γ的取值。
為了計算最優損失函數,本文使用小批量隨機梯度下降(stochastic gradient decent,SGD)來優化聚類損失和重構損失。特別地,有兩種參數需要優化:自動編碼器的權值和聚類中心。Lc相對于嵌入點zi和聚類中心 μj的梯度如式(7)~式(8)所示:
![]() |
![]() |
其中,K是最大的聚類簇數,n為樣本總數,則對于小批量樣本m和學習率λ,更新后的 μj如式(9)所示:
![]() |
解碼器的權重 由下式進行更新,如式(10)所示:
![]() |
因此,編碼器的權重W通過下式進行更新,如式(11)所示:
![]() |
1.3 O6-CMG位點的定位
上述策略通過自動編碼器和聚類算法的聯合調優,構建了nano-UDL模型。接下來則需要nano-UDL模型對O6-CMG數據集進行聚類以實現對包含突變位點序列的篩選。由于O6-CMG數據集包含了兩個類別(類別A和類別B),因此O6-CMG數據集被聚類成為兩個簇分別對應類別A和類別B,如圖1(c)所示。
nano-UDL方法初步實現了對O6-CMG子序列的分類,而為了精確定位子序列中O6-CMG的位置,進一步需要對分類好的類別A中的序列執行分段計算。取類別A所在簇的其中一個序列S,將其分成n段等長小段序列[s1, s2, , sn],對這n個小段計算幅值變化率,取幅值變化最快的段即為O6-CMG位點所在的位置。具體計算方法如式(12)所示:
![]() |
其中,grad(·)是計算幅值變化率函數,本文統一取n = 100,通過上式計算最大的梯度所在的段即為當前序列S的O6-CMG位點(loc)所在的位置。
通過以上三個步驟(數據生成、nano-UDL方法的構建、O6-CMG位點的定位)可以實現在不需要人工干預的情況下,完成對O6-CMG數據集中的O6-CMG突變位點的精準定位。
2 實驗
本文在O6-CMG數據集上評估nano-UDL方法的性能,并與其他傳統的聚類算法進行比較。在本節中,不僅通過四個評價指標來展示nano-UDL方法的性能,還進行超參數敏感度驗證和消融實驗,以此驗證模型的魯棒性。
2.1 評價指標
評價指標是評價算法性能的指標。由于本文主要探討的是無監督學習,因此本文給出了四種標準的無監督評價指標,用于與其他算法進行評估和比較。基于數據集的類型,聚類簇的數量設置為2,本文用無監督聚類精度來評估性能,即準確率(accuracy,ACC),其計算如式(13)所示:
![]() |
其中,yi 為真實標簽,ci為算法產生的聚類分配,m為聚類與標簽之間可能的映射,n為樣本總數,1{·} 是非0即1判斷器,如果映射后的樣本標簽和真實標簽一致則為1,否則為0。
歸一化互信息(normalized mutual information,NMI)(符號為:NMI)用以計算來自同一數據集的兩個標簽之間相似度的歸一化度量,定義如式(14)所示:
![]() |
其中,為真實標簽l與預測聚類標簽c之間的互信息,H為其信息熵。NMI的結果不因簇(類)的排列而改變,它們被歸一化到[0, 1],0表示不相關,1表示完全相關。
另外一種評價指標是調整后的蘭德指數(adjusted rand index,ARI)[38],它經常用于聚類驗證,因為它是兩個類別之間的一致性的度量。
卡帕(Kappa)系數是新的聚類性能度量指標[39]。Kappa系數是用于一致性檢驗的指標,也可以用來衡量分級效果。因為對于分類問題來說,所謂一致性就是模型預測結果是否與實際分類結果一致。Kappa系數的計算是基于混淆矩陣的,混淆矩陣的范圍在?1~1之間,通常大于0。
基于混淆矩陣的Kappa系數計算公式如式(15)所示:
![]() |
其中,K表示為Kappa系數;po為觀測到的符合率,即ACC;pe表示隨機匹配率,即對應于所有類別的“實際數量與預測數量的乘積”之和除以“樣本總數的平方”。
2.2 實驗細節
一般監督學習使用驗證集上的交叉驗證來確定超參數,而非監督學習則不需要設置驗證集,因此本文使用默認超參數,盡量避免調優。本文將O6-CMG數據集的網絡尺寸依次設置為256、256、512和4,所有層都是全連接。
在預訓練階段,從均值為0和標準差為0.01的高斯分布中選取隨機數作為初始值。通過1 000次迭代對每個層進行預訓練。為了對整個自動編碼器進行參數調優,網絡執行2 000次的迭代。全部層都進行預訓練和調優,批次大小設置為256,開始學習率設置為0.1,每200次迭代衰減10%。參數γ將在第2.6節中進行討論。通過設置上述參數,使重構損失達到較好的效果。
為了初始化聚類中心,執行Meanshift聚類算法 20次,然后選擇最佳的聚類結果。為了使KL散度最小,設置學習率為0.01,收斂閾值設置為0.001。對于所有的基線算法,本文對參數進行調整,并選擇效果最好的結果進行比較。
本文主要目標是能夠最大限度地將突變段序列從整個分割后的序列中找出來,同時盡量降低噪聲序列被判斷為突變序列的情況,所以在最終的評價指標中,選擇使用最基本的二分類指標ACC作為初步評價指標,進一步增加了包括NMI、ARI、Kappa三種應用于衡量聚類效果的度量指標。
2.3 實驗結果
本文對不同的算法的性能進行了定量和定性的評估,如圖3所示,其中前5種對比方法的輸入特征是通過人工手動獲取,提取的特征包括諸如均值、方差等統計信息特征。對比方法的輸入特征是相同的,而本文所提出的nano-UDL方法輸入的特征是通過如圖1所示的自動編碼器優化獲取。通過觀察發現,在NMI、ARI、Kappa、ACC等所有評價指標上,nano-UDL方法均明顯高于其他方法,并在O6-CMG數據集上表現出良好的聚類性能。這是因為nano-UDL方法使用多層自動編碼器作為特征提取器,能夠捕捉到非深度學習模型無法獲得的局部特征。此外,圖3中前5種方法都是非深度學習模型,nano-UDL方法是基于深度學習模型的。由此可以觀察到基于深度學習模型的方法明顯比非深度學習的方法能發揮更好的性能。

自動編碼器和不同聚類算法結合的實驗結果對比如圖4所示。圖4中前4個柱狀圖都是自動編碼器與對應聚類算法結合后的結果,結果表明,nano-UDL方法優于其他組合。本文選用的MeanShift聚類算法針對于O6-CMG數據集有著良好的聚類效果,其得益于O6-CMG數據在通過自動編碼器特征提取后的特征分布更適合這種基于密度類型的聚類方法。

2.4 超參數敏感性
為了研究nano-UDL方法在O6-CMG數據集上的參數敏感性,本文對不同參數下的ACC進行柱狀圖直觀對比,如圖5所示。結果表明,nano-UDL方法在大多數參數組合下都保持了良好的結果,并且相對穩定。進一步地,對其他幾種度量指標也同樣做了類似的實驗,實驗結果同圖5類似,同樣證實了nano-UDL穩定性。

2.5 消融實驗
本文方法所采用的損失函數是聚類損失和解碼器重構損失的線性加權和。因此本文主要是通過對聚類損失和解碼器重構損失進行的聯合優化,使nano-UDL方法的自動編碼器模型得到最優的聚類結果。為驗證重構損失添加之前和之后對實驗結果的影響,本文執行了一個消融實驗。其中在“無重構損失”的情況下,ACC為98.12%,聚類時間為1.735 s,相對而言在“有重構損失”的情況下ACC為98.32%,聚類時間為2.228 s。實驗結果表示,“有重構損失”比“無重構損失”發揮出了更好的性能。對于聚類時間而言,“有重構損失”聚類時間略長于“無重構損失”,原因是計算解碼器重構損失需要額外時間。總的來說,將聚類損失和重構損失相結合時nano-UDL方法性能最好。
2.6 參數γ的確定
為了解如式(6)所示的聚類損失參數γ對nano-UDL方法性能的影響,本文在 [0.01, 100] 范圍中采樣,選取13個值。優化器設置為隨機梯度下降法(stochastic gradient descent, SGD),動量β=0.9,學習速率λ設為0.1。結果如圖6所示,可以觀察到當γ值過小時,聚類損失項的積極作用被消除,而值過大則會導致潛在特征空間的扭曲從而影響性能,其中橫坐標值做了對數拉伸處理,以便橫坐標可以等比例顯示。當γ接近0.5時,聚類結果可以達到最優。因此,本文在所有的實驗中設置γ=0.5。

2.7 與監督方法對比
在本節中,將nano-UDL方法與其他監督學習算法進行比較,驗證的算法有分類回歸樹[40](classification and regression trees,CART)、SVM[41]、AdaBoost[42]、隨機森林(random forest,RF)[43]、K近鄰(k-nearest neighbor,KNN)和卷積神經網絡(convolutional neural networks,CNN)。這些監督算法均具有不同的監督率,為此,本文在常規監督學習算法上使用0.1~0.9比例的訓練集用于模型訓練。除此之外,本文利用CNN模型來驗證傳統深度學習模型的性能。實驗結果如圖7所示,隨著監督率的增加,模型的分類效果也在逐步提高。雖然在完全監督情況下,nano-UDL方法不如某些全監督算法。但是,在沒有監督或者監督信息很少的情況下,本文所提出的nano-UDL方法具有顯著的優勢。除CNN和nano-UDL方法這兩種深度學習模型外,其余模型輸入采用相同的特征,從圖7中可以看出CART、SVM、AdaBoost與KNN在全監督情況下不如nano-UDL方法。也就是說,手工提取的特征用上述這四種模型訓練并不能達到預期的效果,但例外的是RF模型在同樣的實驗設置下卻可以表現得較好。這也從另一個側面驗證了深度學習相比于手工提取特征技術,可以針對性地對特定任務進行特征調整。與傳統機器學習相比不同的是,在某些特征情況下,模型間可能存在較大的性能偏差。而本文提出的nano-UDL方法在無監督的背景條件下,也能達到部分有監督信息訓練算法的性能,這也證實了無監督深度學習在納米孔數據分析領域具有一定的應用價值。

3 討論
為了證明本文所提出框架的有效性,本文采用了消融實驗、精度對比等機器學習模型驗證策略對nano-UDL方法進行了驗證。通過對實驗輸出結果進行分析,本文從以下兩個方面對nano-UDL方法進行討論:
可靠性:nano-UDL方法有效地保證了突變序列檢測的高精度和高置信度,這是構建可靠且穩定系統的必然要求。如圖3所示的nano-UDL方法的結果比一些傳統的聚類算法具有更高精度和置信度。從圖5的結果可以看出超參數對聚類結果的影響很小,充分說明了本文提出模型的穩定性。通過改變監督樣本的大小,從圖7可以看出,nano-UDL方法的結果相對較好而且相對穩定。綜上所述,實驗結果充分說明了nano-UDL方法的可靠性。
訓練代價:雖然nano-UDL方法相比于其他聚類方法有著較好的性能,但由于這種方法策略是基于多梯度進行計算的,其中聚類損失和重構損失,增加了訓練的時間復雜度。在訓練階段,本文使用少量的樣本對模型進行預訓練,以減少模型訓練的成本。通過使用預先訓練的模型初始化網絡,可以有效地保證模型在訓練過程中時刻處于最佳。由于本文的訓練數據集樣本較少(1 010例),所以消耗的系統內存、GPU資源等計算資源都在有效范圍內。
目前的納米孔數據分析技術都是基于統計學習方法或有監督的機器學習算法,不能保證時間的有效性,因為其需要大量的分析計算成本。本文提出的方法是一種可以降低學習代價的無監督學習算法,該方法對納米孔DNA突變的研究有較好的效果,將對今后更多納米孔DNA突變的研究提供指導思路。本文提出的方法是嘗試將無監督深度學習方法應用到納米孔數據的首次探索,希望這種在生物信息學應用中利用先進的機器學習技術的思想在未來能夠逐漸走向成熟。
4 結論
針對O6-CMG納米孔數據,本文提出了一種新的無監督深度聚類模型nano-UDL,該模型將數據點聚類到一個可以通過聯合優化的特征空間中,最終實現精準識別突變位點的目的。為了逼近目標分布,該方法通過最小化聚類損失和重構損失進行迭代訓練。nano-UDL方法可以看作是一種不需要任何人工干預的啟發式工具,它可以用來處理原始納米孔數據。通過實驗分析,nano-UDL方法在提高性能的同時,對超參數設置具有魯棒性。nano-UDL這種無監督深度學習策略可以在減少人工計算成本的情況直接檢測出O6-CMG突變位點。期望本文提出的方法將對未來更多類似的從DNA測序中精準定位突變位點的任務提供一個理論依據。
重要聲明
利益沖突聲明:本文全體作者均聲明不存在利益沖突。
作者貢獻聲明:關曉宇主要負責項目主持、算法程序設計以及論文編寫;王宇和張金月主要負責實驗數據的采集和分析; 邵偉, 黃碩, 張道強主要負責提供實驗指導,數據分析指導,論文審閱修訂。