陳 城，林 劼

(福建師范大學數學與統(tǒng)計學院，福建 福州 350117)

生物測序技術的發(fā)展產生大量的生物序列數據，對序列進行聚類分析，把功能相近的序列聚為一類，可以幫助科研人員快速了解生物序列的功能，為后續(xù)的研究奠定堅實的基礎.聚類在生物信息學中并不是一個新話題.基因表達數據的分析常對微陣列數據進行聚類，并認為在同一個簇中的基因具有相同的功能.學者們提出了幾種先進的微陣列聚類算法[1]，包括層次聚類[2]和集成聚類[3].為了避免微陣列數據中的噪聲，雙聚類[4]的方法應運而生，它能夠同時進行特征選擇和樣本選擇.單細胞測序數據的工作原理類似于基因表達數據，根據不同的細胞類型[5-6]對不同的基因表達進行聚類.然而，這些方法關注的是基因表達數據，而不是生物序列數據.

生物序列聚類是聚類算法在生物學中的應用，使用的數據為：核酸序列(如DNA、RNA序列)或氨基酸序列(如蛋白質序列).其中，核酸序列由4種堿基組成：A(腺嘌呤)、G(鳥嘌呤)、C(胞嘧啶)、T(胸腺嘧啶)/U(尿嘧啶)，蛋白質序列則是由20種氨基酸排列組合而成.序列聚類需要將序列轉化為數值向量.雖然有一些工具可以將DNA/RNA/蛋白質序列轉化為數值向量[7-10]，但對序列進行表示學習仍然是聚類過程中重要的一步.生物序列聚類研究構建衡量序列間的相似性或差異性的數學模型，旨在根據統(tǒng)計結果將序列劃分為幾個簇，同一簇中的序列具有相似的功能，這樣可以由已知序列的功能推測未知序列的功能[11].

生物序列聚類有利于分析序列潛在的結構、功能等信息，進一步推演序列在進化過程中發(fā)生的先后關系[12].迄今為止，學者們已經提出眾多優(yōu)秀的聚類算法，但是由于領域的差異、解釋性差等問題，不能直接應用在生物序列聚類上.目前，生物序列聚類的困難之處在于生物序列數據的特征難以正確提取，數據量大，計算復雜度高，內存需求高，不能保證一定能夠找到最優(yōu)解且很難從生物學角度進行結果解釋等.

本文提出一種基于多重解碼器的自編碼器模型，用于學習生物序列數據的表示，然后使用 k-means算法對序列的表示進行聚類. 試驗結果驗證本文提出的方法在 DNA 序列數據集上具有良好的性能.

1 相關工作

根據聚類方法的不同，生物序列上的聚類研究大致可分為：啟發(fā)式生物聚類算法、層次生物聚類算法和其他生物聚類算法.

1.1 啟發(fā)式生物聚類算法

基于啟發(fā)式的生物聚類算法采用一種簡單貪婪策略，始于一個種子(seed)，基于一定的搜索技術，合并、擴張，完成序列的聚類過程.其代表性算法為FastGroup[13]、CD-HIT[14]、UCLUST[15]等.

2001年，Seguritan首次提出一種基于啟發(fā)式的生物序列聚類算法FastGroup.該算法主要有3個步驟：首先將數據集中的所有序列相互比較；然后將相似的序列分組；最后從每組中輸出一個代表序列.CD-HIT使用貪婪的增量聚類算法，該算法通過短詞過濾策略和并行化技術得到提高，有效地對大規(guī)模序列數據集進行聚類.短詞過濾策略的整體思想為：若序列與代表序列的相似性低于預先設置的閾值，則該序列不進行序列比對.該策略可以使算法的運行速度加快，其復雜度為O(N)，N為序列數量.

USEARCH算法在進行比對時會在所有相似度達到閾值的序列中尋找合適的比對位點，基于這個特點，UCLUST提出一種改進方法，僅需要尋找一個或幾個合適的比對位點，減少了比對位點的搜索數量，提高了序列比對速度.與CD-HIT相比，這種方法的運算速度快、內存需求低，且靈敏度較高.

Ghodsi等[16]提出了新的詞過濾方法DNACLUST，避免了序列間的兩兩比對，是一種貪婪算法.該算法在精確模型下的運算速度優(yōu)于CD-HIT與UCLUST，但在近似模型下的運算速度與UCLUST差不多.SEED[17]則使用開放哈希技術和一種特殊類型的稱為塊間隔的種子將輸入序列聚類.

UPARSE[18]是來自USEARCH的最新從頭聚類方法，它通過質量過濾算法過濾read，將其修整為固定長度.該方法使用UPARSE-OTU進行聚類，這是一種新的貪婪算法，可同時執(zhí)行嵌合過濾和OTU聚類.由于質量過濾算法的嚴格，生物物種豐富度和多樣性會被顯著低估，最后產生的生物學結果存在較多錯誤.大規(guī)模過濾序列能夠減少運行時間，但是過濾參數需要針對數據的不同進行改變，且不能自動對過濾參數進行選擇，需要多次試驗后人為進行選擇，而這個過程的時間成本可能很高.

在進行聚類的過程中，啟發(fā)式算法不需要計算距離矩陣，因此降低了存儲空間的需求和計算復雜度.缺點是不能保證一定能夠找到最優(yōu)解.

1.2 層次生物聚類算法

基于層次的生物序列聚類算法通過序列比對，根據一定的相似性或距離度量方式獲得距離矩陣，再采用貪婪層次聚類算法完成生物序列聚類，是目前最常用的生物序列聚類方法，其代表性算法為DOTUR[19]、Mothur[20]、ESPRIT[21]、mBKM[22]等.

2005年，Schloss提出一種無種子(seed)的生物序列聚類算法DOTUR.DOTUR使用序列之間的遺傳距離將序列分配給OTU.它通過使用最遠、平均或最近鄰居算法為OTU分配序列，并估計一個簇的豐富性和多樣性.

在DOTUR的基礎上，Schloss又開發(fā)了Mothur.Mothur被用于修剪、篩選和排列序列，計算距離，將序列分配給OTU，Alpha和Beta多樣性計算，序列比對，序列聚類注釋.距離矩陣對于Mothur聚集OTUs很重要.矩陣可以反映每個序列與其他序列的相似性或距離.計算矩陣和聚類具有較高的時間復雜度.基于此，Mothur不適合處理大數據集[23].該算法存在假陽性率高、噪音信號強等缺點，且很難從生物學角度進行結果解釋.

Sun采用ESPRIT算法對生物序列進行聚類劃分.該算法由4個模塊組成：(1)使用各種標準去除低質量read；(2)計算read的成對距離；(3)將read分組到不同差異水平的OTU中;(4)執(zhí)行統(tǒng)計推斷來估計物種豐富度.該算法使用了Needleman-Wunsch雙序列比對算法，以k-mers的形式過濾無需比對的序列，利用在線學習開發(fā)了一種名為Hduster的方法進行層次聚類，在一定程度上降低了計算機的內存需求，但計算復雜度為O(N2)，當N較大時，算法的時間成本高.

Cai[24]基于ESPRIT算法提出了一種新的在線學習的算法ESPRIT-Tree.基本思想是使用偽度量構造分區(qū)樹，利用分區(qū)樹將序列空間劃分為一組子空間，然后遞歸地細化這些子空間中的簇結構.該技術依賴于快速最近對搜索和一種動態(tài)插入和刪除樹結點的方法.為了避免窮舉計算簇之間的成對距離，該方法將序列的每個簇表示為概率序列，并進行一系列操作來比對這些概率序列并計算它們之間的遺傳距離.ESPRIT-Tree是啟發(fā)式生物序列聚類算法的一種，該算法同時解決了計算復雜度高和計算機內存需求多的問題，其計算復雜度幾乎與生物序列數目呈線性關系.

mBKM是一種基于新的距離度量DMk的非比對算法，用于聚類基因序列.該方法將DNA序列轉化為特征向量，其中包含DNA序列中k-mer的出現次數、位置和順序關系.然后，將層次聚類算法應用于DNA序列.研究表明[25-26]，當基于同質性標準對數據集進行劃分時，該方法得到了較好的聚類結果.

層次生物序列聚類算法的缺點是在擁有龐大的序列數據時，所需要的儲存空間大，計算復雜度高.

1.3 其他生物聚類算法

其他生物聚類算法，如CROP[27]算法利用等級機制將需要比對的序列劃分為若干個子集，然后基于貝葉斯理論，采用高斯混合模型對序列進行聚類.在生物序列中，CROP算法可以推斷出最優(yōu)的聚類結果，其中高斯模型的抗噪聲能力能夠克服由于測序誤差而導致對序列的高估、內存需求高及計算效率低的問題，在一定程度上能夠較好地實現生物序列聚類，具有較強的抗噪聲能力和魯棒性.

CBE[28]是一種基于最大熵原理的聚類方法.這種方法基于數據的先驗信息來探索數據所有可能的概率分布空間，得到熵最大的分布，當熵值達到最大時，聚類結束.先驗信息基于以下假設：根據某種統(tǒng)計方法，簇中元素彼此相似.基于此，滿足條件的那些高熵分布優(yōu)于其他分布.

coreClust[29]是一種基于檢測保守區(qū)域的非比對聚類方法.這些區(qū)域的檢測可用于功能注釋和分組蛋白質序列的區(qū)域.coreClust基于一種名為MinHash的技術，該技術是一種局部敏感哈希方法，用于識別集合中的相似元素.它主要依賴于哈希，因此該方法非常適合MapReduce并行處理平臺，從而實現可伸縮性.

DeLUCS[33]模型使用DNA序列的頻率混沌博弈表示(frequency chaos game representations，FCGR)，并通過優(yōu)化多個神經網絡生成模擬序列的FCGRs來學習數據的模式(基因組簽名).然后使用多數投票方案來確定每個序列的最終簇分配.ALFATClust[34]利用快速成對非比對的序列距離計算和社區(qū)檢測來生成簇.ALFATClust可以通過考慮簇的分離和簇內序列相似性來動態(tài)確定生成每個簇的閾值，而不是對每個生成的簇應用單個閾值.

2 方法

本文提出一種基于多重解碼器的自編碼器模型，用于生物序列數據的表示學習，然后使用k-means算法對序列的表示進行聚類.將整個過程分為兩個階段：表示學習階段和聚類階段.

2.1 表示學習

在進行聚類之前，需要對生物序列進行表示學習，得到序列所對應的表示，以便進行后續(xù)的聚類任務.本方法使用的表示學習模型如圖1所示.首先，對生物序列數據進行k-mer劃分并統(tǒng)計k-mer的頻率，得到序列k-mer的頻率向量作為模型的輸入；其次，在自編碼器模型的基礎上應用多重解碼器的形式，從多個不同的空間對特征進行重構；最后，將序列的k-mer頻率向量輸入到模型中使用均方差損失進行訓練得到序列的表示.

圖1 基于多重解碼器的自編碼器模型

2.1.1 k-mer的提取和頻率向量

k-mer是一段長度為k的子串，它是由序列剪切一部分得到的.k-mer的提取指的是將一條序列連續(xù)切割，按照堿基的排列順序依次劃動得到的一系列長度為k的子串.

給定生物序列集S={S1，S2，…，Si，…，SN}，序列Si的長度為M，在序列上通過一個長度為k(2≤k≤M-1)、步長為1的滑動窗口從序列中提取長度為k的子串，即k-mer.在長度為M的序列中存在M-k+1個k-mer.一條具有|Σ|個字母的序列最多有|Σ|k種不同的k-mer.一條DNA序列由字母表中的字母組成，字母表Σ={A，C，G，T}且|Σ|=4.

使用單個字母的頻率作為特征表示樣本時，所有的關聯信息會被丟失.k-mer頻率向量通過描繪序列局部信息將每條樣本轉化成|Σ|k維向量，從而彌補這一缺陷.因此，使用k-mer頻率向量來表示樣本序列.

具有|Σ|k維的k-mer頻率向量可以表示為：

Xi= [Xi1，Xi2，…，Xij，…，Xi|Σ|k]T,

(1)

其中，Xij是第i條序列中第j種k-mer出現的頻率，Xij可由公式(2)計算得到：

(2)

其中，Nij是在第i條序列提取的所有的k-mer中第j種k-mer出現的次數，Ni是第i條序列中提取的所有k-mer的數量.

2.1.2 編碼器

對生物序列Si提取k-mer得到頻率向量Xi后，使用編碼器E(·)對k-mer頻率向量Xi進行特征提取，如公式(3)所示：

Zi=E(Xi)，i=1，2，…，N,

(3)

其中，Zi表示序列Si的特征向量.在本文中，使用編碼器E(·)對序列數據進行映射，得到序列的表示.它需要與解碼器Dl(·)一起使用均方差損失進行訓練.訓練完成后，不再使用解碼器的部分，而是單獨使用編碼器.由編碼器得到的表示可以用于下游任務，本文中，應用學習到的表示進行聚類.

2.1.3 多重解碼器

與傳統(tǒng)的自編碼器模型不同，本文采用多重解碼器的形式，具有提高編碼器表征能力的優(yōu)點.給定L個解碼器Dl(·)，l=1，2，…，L，對于序列Si的特征向量Zi，使用解碼器對其進行解碼，重構序列Si的k-mer頻率向量Xi，如公式(4)所示：

(4)

2.1.4 學習策略

綜上所述，模型是由一個編碼器和多個解碼器組成.訓練時，分別使用多個解碼器的輸出和序列的k-mer頻率向量計算均方差損失，解碼器Dl(·)的損失函數Ll如公式(5)所示：

(5)

最后，同時最小化解碼器的均方差損失，總損失函數如公式(6)所示：

(6)

模型訓練完畢后，使用編碼器的輸出Zi作為序列的表示進行后續(xù)的任務.

2.2 聚類

為了推測未知生物序列的功能，分析基因間進化的先后順序關系，對生物序列進行聚類有助于科研人員快速了解生物序列，加快推動相應研究的進展.本文使用的聚類算法是k-means算法.該算法具有收斂速度快的優(yōu)點，并且能夠達到較好的聚類效果.

k-means算法是最常用的聚類算法，主要思想是在給定k值和k個初始簇中心點的情況下，把每個點(即樣本)分到離其最近的簇中心點所代表的簇中，所有點分配完畢后，根據一個簇內的所有點重新計算該簇的中心點(取平均值)，然后再迭代進行分配點和更新簇中心點的步驟，直至簇中心點的變化很小，或者達到指定的迭代次數.

給定樣本的表示Z(即為編碼器的輸出)，包含了n個樣本的表示Z={Z1，Z2，Z3，…，Zn}，其中每個表示都具有m個維度的特征.對于k-means算法，首先需要初始化k個簇中心{C1，C2，C3，…，Ck}，1

(7)

其中，Zi表示第i個對象，1≤i≤n，Cj表示第j個簇中心，1≤j≤k，Zit表示第i個對象的第t個分量，1≤t≤m，Cjt表示第j個簇中心的第t個分量.

依次比較每一個樣本到每一個簇中心的距離，將對象分配到距離最近的簇中心的簇中，得到k個簇{S1，S2，S3，…，Sk}.

聚類時，目標是最小化每個樣本到其對應的簇中心的距離，如公式(8)所示：

(8)

基于多重解碼器的自編碼器模型的生物序列聚類算法，如算法1所示.

算法1 基于多重解碼器的自編碼器模型的生物序列聚類算法輸入:生物序列集S,編碼器E(·),解碼器Dl,序列數量N,解碼器數量L輸出:生物序列的簇標簽1:for i=1 to Ndo2: Xi=從Si中提取kmer的頻率向量3: Zi = E(Xi)4: for l=1 to L do5: ^Xli=Dl(Zi)6: end for7:end for8: 計算損失 Ltotal如公式(6)所示9: 使用隨機梯度下降法計算梯度10:更新模型參數,直到收斂11:得到訓練完畢的編碼器 E(·)12:for i=1 to N do13: Zi= E(Xi)14:end for15: 使用k-means 對序列表示集 Z進行聚類,得到生物序列的簇標簽

3 數據試驗

3.1 數據集

HOGENOM是一個全序列生物同源基因家族的數據庫，其中包含來自真核生物、細菌和古細菌的182個完整基因組的同源基因家族.本選題使用從HOGENOM中隨機抽取的HOG100、HOG200、HOG300共3個DNA序列數據集，其對應的序列類別數分別為100、200和300，這些DNA數據由蛋白質生物學和化學研究所(IBCP)使用基于人口的增量學習(PBIL)進行收集，可以從ftp:∥pbil.univlyon1.fr/pub/hogenom/release—06/獲得相關數據.表1列出了這3個DNA數據集的詳細信息.

表1 DNA數據集的詳細信息

3.2 參數設置

本文選擇的k-mer的k為3.編碼器的網絡結構為3個全連接層，神經元參數分別設置為128、64和16，激活函數均為relu函數.解碼器D1(·)的網絡結構為4個全連接層，神經元參數分別設置為128、64、32和64，激活函數均為relu函數；解碼器D2(·)的網絡結構為5個全連接層，神經元參數分別設置為128、64、56、32和64，激活函數均為relu函數；解碼器D3(·)的網絡結構為5個全連接層，神經元參數分別設置為128、256、128、32和64，第一層和最后一層的激活函數為tanh函數，其余均為relu函數.

3.3 評價指標

本文采用的評價指標是歸一化互信息 (normalized mutual information，NMI)，NMI是聚類評價指標中的一種外部指標，用于度量聚類結果與真實標簽之間的關系，其范圍為[0，1].當聚類的結果越接近真實的情況，NMI的數值就越接近1.若兩者之間相互獨立，那么NMI就為0.NMI的定義如公式(9)所示：

(9)

其中，I表示互信息(mutual information)，H為熵.

(10)

其中，P(wk)、P(cj)、P(wk∩cj)可以分別看作樣本屬于聚類簇wk，屬于類別cj，同時屬于兩者的概率.第二個等價式子則是由概率的極大似然估計推導而來.

(11)

3.4 試驗結果

表2 聚類結果的比較

從試驗的結果可以看出，本文提出的方法在HOG100數據集上具有競爭性，在其他2個數據集上性能優(yōu)于另外2個方法.隨著數據集類別數量的增加，本方法的表現更加穩(wěn)定，其他2種方法的性能反而出現下滑的趨勢，說明本方法在種類數量大的序列數據集中一樣適用.

4 結論

本文提出了一種基于多重解碼器的自編碼器模型的聚類方法.在自編碼器的基礎上進行了改造，增加了解碼器的數量，從不同的空間對輸入數據進行了重構，從而加強了模型進行表示學習的能力；接著使用模型中編碼器的輸出作為k-means的輸入來進行聚類.試驗表明本方法獲得了良好的聚類效果，而且在當前數據集下隨著數據集的類別數量的增加，本方法的性能將趨于穩(wěn)定.

未來會進一步針對解碼器研究適用于不同任務的學習策略，繼續(xù)提高編碼器的表征能力，使得學習到的嵌入根據學習策略的不同能夠滿足一個或多個下游任務的需求.