基于流形嵌入的宏基因組疊連群分箱方法研究

何 翀 王美麗，2，3 景 旭*

1(西北農林科技大學信息工程學院 陜西 咸陽 712100)2(西北農林科技大學農業農村部農業物聯網重點實驗室 陜西 咸陽 712100)3(西北農林科技大學陜西省農業信息感知與智能服務重點實驗室 陜西 咸陽 712100)

0 引 言

微生物是自然界和人類不可或缺的一類生物[1]。近來研究表明，微生物不僅在自然環境中扮演著重要角色，人體內腸道微生物和多種疾病也存在著某種聯系[2]，如：肝硬化[3]、肥胖癥[4]、糖尿病[5]、腫瘤[6]，甚至神經行為類疾病和免疫類疾病[7-8]。科學家們一直在研究如何更好地培養和研究微生物，但是因為技術條件的限制，99%以上的微生物仍然是不可培養的[9]。這導致對于大多數的微生物，人們無法使用傳統方式研究它們。

高通量測序技術的不斷成熟和成本不斷的降低，產生了許多不同視角的研究方法。宏基因組學使用高通量測序技術直接對特定環境進行測序分析，不需要對微生物進行培養，極大地擴充了微生物學的研究對象。但是由于測序技術的限制，當前測序得到的讀長往往缺失測序對象的來源信息[10]，而且宏基因組讀長的組裝往往也無法得到完整的基因組，只能得到長片段的疊連群，這會對下游分析產生較大的影響。因此，眾多研究人員提出宏基因組學疊連群分箱方法以將宏基因組疊連群分發到操作分類單元(Operational Taxonomic Units,OTU)中，為下游分析提供更為明確的研究對象和更方便的操作對象。在疊連群分箱問題中，基于單拷貝基因確定初始化分箱數的效率和分箱器的效率較低。所以，高效準確的疊連群分箱成為了宏基因組學研究中的重要問題之一。

1 相關工作

宏基因組學疊連群分箱方法可分為依賴其他分箱方法的集成方法和不依賴其他分箱方法的獨立方法。獨立方法是基于某種數學模型來構建分箱器，比如：高斯混合模型(Gaussian Mixture Model,GMM[11])、期望最大化(Expectation Maximization,EM[12])、k-medoid[13]、仿射傳播算法(Affinity Propagation,AP)[14]和非負矩陣分解(Non-negative Matrix Factorization,NMF)[15]等。CONCOCT[11]使用疊連群的覆蓋度和四核苷酸頻率組成來提取聚類特征，接著使用高斯混合模型來對疊連群進行聚類。作為進一步的拓展，CONCOCT[11]可通過疊連群邊緣讀長的連接信息來進一步地整合前面聚類后的簇。混合高斯模型確實可以聚類得出一些OTU單元，然而混合高斯分布這一假設并非在所有數據上都成立。MaxBin[12]、MetaBAT[13]使用疊連群覆蓋度和組成信息聯合現存微生物的基因組數據學習得到四核苷酸頻率概率距離(Tetranucleotide frequency probability distance,TDP) 和覆蓋度距離概率(Abundance Distance Probability,ADP)，再使用這兩者來構建實驗數據中疊連群之間的距離，最后使用改進的k-medoid方法進行分箱。雖然利用現存的監督數據學習判別距離可以極大地提高分箱器的學習效率，但是監督模型受限于現有的基因組數據，模型學到的四核苷酸頻率距離的判別性能也受限(并不能分離出四核苷酸頻率相近的分類單元)。BinSanity[14]使用覆蓋度信息和仿射傳播算法，迭代地將每一條疊連群當作聚類中心，按照疊連群之間的相似度進行分箱。對AP算法的聚類結果中高冗余或低完整度的部分，利用四核苷酸頻率和GC含量等信息進行提純等后期處理。AP算法提高了算法的可適應性，減少了用戶的輸入參數數目，但是由于AP算法的時間復雜度很高，可拓展性較差，在實際應用中的可用性較低。文獻[16]和文獻[17]中使用關聯性特征進行分裝，但是對于豐度不均勻的數據效果較差。COCACOLA[15]使用非負矩陣分解的方法對疊連群進行分箱。矩陣分解提供了一個不同的分箱視角，使用序列組成、讀長的覆蓋度、疊連群比對結果和雙端讀長連接信息來構建特征矩陣，利用NMF算法對特征矩陣分解，最后得到每個疊連群的聚類簇。COCACOLA提出了一個新的可以融合較多信息的宏基因組學分箱框架。BMC3C[18]使用序列組成、覆蓋度信息和密碼子信息作為疊連群特征。通過多輪KMeans聚類得到的結果構建圖模型，使用圖分割算法進行劃分，將劃分結果中的子圖作為分箱的結果。SolidBin[10]使用半監督學習的方式將約束信息引入譜聚類方法中，根據不同的約束信息拓展出了幾種變體，使用歸一化切割方法分箱。雖然不同的變體適應于不同的應用場景，然而如何選取恰當的分箱參數α、β、K對于用戶較為困難。

使用組合策略的集成方法，利用其他軟件的聚類結果，對這些結果進行進一步的優化，比如挑選出其中質量或完整度較高的結果，剔除污染度較高的部分，或使用組裝前的讀長數據對這些分箱好的簇進行讀長的募集，以圖提高聚類結果的完整度。DASTools[19]使用基于單拷貝基因的評估函數對其他分箱軟件的分箱結果進行評分，然后迭代地從中選取最好的結果，最后集成這些結果形成分箱單元。Binning_refiner[20]的思想是將三種方法的分箱結果合并在一起，然后對合并文件進行比對，對比對結果設定一些閾值，篩選出三種方法共有的即提純后的簇。metaWRAP[21]在分箱上的想法是，比對以上兩種方法(DASTools,Binning_refiner)，在融合集成其他聚類結果的基礎上，對融合成的簇進行讀長募集，進一步地提高簇完整度。metaWRAP擴大了生物學分類單元的范圍，因為募集的過程可能會把近同源的讀長募集進簇，另外方法過程中涉及很多比對過程，耗時較長。

對于分箱的初始化，常用的方法有下面幾種：變分貝葉斯自動推斷[11](Variational Bayesian estimation of a Gaussian mixture,VBG)，使用功能單拷貝標記基因來確定初始聚類數或者使用集成聚類的方法確定最佳聚類數[18]。雖然以上分箱器使用這幾種方法，但是這些方法存在著一些問題，如CONCOCT[11]中使用的變分貝葉斯方法并不總能夠得到理想的結果，基于單拷貝標記基因的方法使用了外部程序對原始的疊連群集進行掃描、翻譯和識別，耗時較長，同樣地，使用多輪聚類方法進行初始化的方法也存在這樣的問題。

針對以上分箱器初始化聚類數效率低下，且未充分利用序列特征流形結構的問題，本文提出基于流形嵌入的分箱方法。流形嵌入是一種用于高維數據的非線性降維方法，降維后的數據不但可以反映出數據內部的隱式結構，為機器學習提供更有效的學習特征，同時也可以用于高維數據的可視化。本文的流形分箱方法可以自適應并高效地估計出初始化分箱數，同時高效地完成宏基因組學疊連群的分箱，如圖1所示。方法描述如下：

圖1 基于流形嵌入的疊連群分箱方法

(1) 對疊連群數據進行四核苷酸頻率和樣本覆蓋度特征提取，進行歸一化后合并得到原數據X。

(2) 對歸一化后的數據進行流形嵌入。學習得到原始數據嵌入的流形Y。

(3) 對嵌入流形上的數據Y進行VBG聚類，得到初始化的聚類數k。

(4) 從k開始增量迭代初始化聚類數，以輪廓系數為聚類性能度量指標，得到最佳的聚類數K。

(5) 以步驟(4)中參數K初始化KMeans聚類算法并將流形嵌入得到的結果Y進行聚類，得到分箱結果。

2 方 法

2.1 預處理

為了更有效地挖掘疊連群中的信息，同時也為得到每一條疊連群標準的特征描述，需要對疊連群數據進行預處理：提取疊連群的序列組成(即四核苷酸頻率)和疊連群的樣本覆蓋度。

2.1.1序列組成

(1)

2.1.2疊連群覆蓋度

(2)

(3)

2.2 流形嵌入

均勻流形近似和投影(Uniform Manifold Approximation and Projection,UMAP)[22]是2018年被提出的一種新的流形學習降維方法。定義d:X×X→R≥0為一種度量，給定一個k，在d下計算xi的近鄰集{xi1,xi2,…,xik},對于每一個xi，定義ρi和σi:

ρi=min{d(xi,xij)|1≤j≤k,d(xi,xij)>0}

(4)

(5)

(6)

B=A+AT-A°AT

(7)

式中：°為Hadamard積。可有對應B的無權圖G。對G的邊施加吸引力(Fattactive)，對邊上的頂點施加排斥力(Frepulsive)，使用模擬退火算法不斷減小吸引力和排斥力，最終算法收斂。吸引力和排斥力定義如下：

(8)

(9)

式中：yi、yj是低維空間中的坐標；a、b為參數；為無窮小量。

最終得到低維空間下的坐標yi∈Rd，所有樣本構成的Y∈Rn×d即流形嵌入得到的結果。

2.3 流形分箱

2.3.1分箱初始化

從疊連群數據集StrainMock的主成分分析降維和流形嵌入降維的可視化結果(如圖2所示，圖中每個點代表一個疊連群)中，發現流形嵌入可以更好地反映出數據集中潛在分箱數。潛在簇的數量是進行聚類數搜索過程良好的初始化，這一點在文獻[22]中也有說明。對于流形嵌入得到的結果Y,使用變分貝葉斯高斯混合模型發現的聚類數k作為初始化分箱數。通過一組聚類簇數初始化KMeans聚類，并計算對應的sihouette系數以評估聚類性能，最終確定最小sihouette系數對應的簇數K為最終的聚類簇數。

(a) StrainMock數據集主成分分析

2.3.2分 箱

對流形嵌入下的結果Y進行分箱。為使分箱方法簡單、有效并易于解釋，以2.3.1節得到的最佳初始化分箱數K為參數，使用KMeans算法對Y進行分箱。KMeans算法將n條疊連群劃分成K個簇，以最小化目標函數：

(10)

作為拓展，可對不同流形嵌入維度下的結果分別使用KMeans算法進行聚類，然后使用集成聚類的一致函數將不同的聚類結果合并，形成最終的分箱結果。

3 實 驗

3.1 實驗環境和數據集

3.1.1實驗環境

實驗運行的環境為Ubuntu 16.04，48核512 GB內存,Intel(R) Xeon(R) CPU E5- 2650 v4,2.20 GHz。本文使用的軟件MaxBin的版本為2.2.7，COCACOLA使用Python版本，SolidBin軟件使用SolidBin-naive變體。

3.1.2數據集

(1) SpeciesMock數據集是基于人類微生物計劃(Human Microbiome Project,HMP)中16S rRNA數據集構建的，它包括96個樣本、101個不同的物種。數據集用來評估在物種水平上的分箱性能。數據集中有37 628條疊連群。

(2) StrainMock數據集是用來測試不同水平上的分箱效果，它包括64個樣本、20個不同的物種或菌株。數據集中有9 417條疊連群。

3.2 評 估

3.2.1分箱初始化

對于分箱初始化聚類數，使用了3.1.2節中兩個標準的分箱數據集進行測試。為了驗證流形嵌入對初始聚類數的影響，對每一個數據集分別使用原數據集、PCA降維結果和UMAP嵌入結果進行測試。然后對不同的變換分別使用XMeans[23]、Meanshift[24]和VBG[25]三種方法計算初始化分箱數，結果越接近于原始的聚類數，說明方法的性能越好。

3.2.2流形分箱

對流形分箱方法，使用3.1.2節中兩個標準的分箱數據集進行測試。分別使用COCACOLA、MaxBin和SolidBin進行分箱，使用標準的評估方法來評估聚類算法的性能：準確度(ACC)，歸一化互信息(NMI)和蘭德指數(ARI)。

準確度：在聚類中，準確度定義是預測標簽和真實標簽的最佳匹配。

(11)

式中：n是樣本數；ci和yi是第i個樣本的預測標簽和真實標簽；Ⅱ為指示函數；m是匹配函數。

歸一化互信息：歸一化互信息可以視作對互信息做歸一化處理的結果。

(12)

式中：y是真實標簽；c是聚類的預測標簽；H代表熵；I是真實標簽和預測標簽之間的互信息。

歸一化蘭德系數：

(13)

(14)

式中：a是在C、K中都在同一集合中的樣本對數；b是在C、K中都不在同一集合中的樣本對數；n為樣本數；RI為蘭德系數；E為期望。

3.3 實驗過程和結果

3.3.1分箱數初始化

實驗首先測試使用貝葉斯高斯混合模型在流形嵌入數據上估計初始化分箱數的效果。實驗中使用了原始數據、PCA變換和UMAP變換后的數據集，在數據變換的結果上使用XMeans、Meanshift和VBG來測試潛在的聚類數即初始化分箱數。

實驗結果標注：實驗項用“X_Y”來標注，其中下劃線前面的字母表示方法，三種方法簡寫如下：X(XMeans)，M(Meanshift)，V(VBG)。下劃線后面的字母表示數據變換：X(Raw data)，P(PCA transform)，U(UMAP embedding)。例如，M_X表示在原數據集上進行Meanshift算法的結果。考慮到UMAP變換中的隨機影響，5次運行的結果如表1和表2所示。

表1 不同方法在StrainMock上的分箱數

表2 不同方法在SpeciesMock上的分箱數

可以看出，在StrainMock數據集中，三種方法在三種變換得到的初始聚類數中最接近真實分箱數(20)的方法就是“V_U”,即在流形嵌入結果下使用VBG進行聚類。雖然XMeans方法較為穩定，但是結果與真實分箱數相差較大,這也說明XMeans方法的自適應效果差。Meanshift方法依賴于帶寬參數的設置，對帶寬參數極為敏感。在Meanshift默認參數下得到的聚類數與真實的聚類數差距較大，在兩個數據集下均未得到理想的初始化分箱數。

在兩個數據集下使用VBG方法進行初始化分箱數估計的實驗結果中，StrainMock數據集下，原始數據集和PCA變換后數據集的結果均為50,這與真實的分箱數相差較大，而在UMAP變換下，VBG得到的分箱數為19～25，更接近于真實分箱數。在SpeciesMock數據集中結果類似。雖然原始數據集和PCA變換后數據集的結果較為穩定，但與真實分箱數(101)相差較大。實驗結果中可以發現，流形嵌入對于VBG估計出良好的初始化分箱數非常重要，雖然UMAP中有隨機過程的影響，但是作為初始化分箱數使用，仍比其他方法要好。

3.3.2流形分箱

如表3、表4所示，在分箱數據中,COCACOLA在原始StrainMock數據集上表現很好，準確率、歸一化互信息和歸一化蘭德指數分別達到了0.974 83、0.957 33和0.950 78。雖然COCACOLA方法較有競爭力，但是本文方法比COCACOLA運行所需時間更少，較MaxBin方法和SolidBin方法相比分箱性能有較大幅度提高，準確率、歸一化互信息和歸一化蘭德指數分別提高0.095 67、0.027 26和0.104 09。

表3 不同方法在StrainMock上的性能

表4 不同方法在SpeciesMock上的性能

在SpeciesMock數據集中，本文方法分箱性能最好，準確率、歸一化互信息和歸一化蘭德指數分別達到了0.998 64、0.998 13和0.997 23,對比COCACOLA的結果有較大程度的性能提高。SolidBin的性能較差，可能是該方法未充分利用四核苷酸頻率和疊連群樣本覆蓋度信息，過度依賴于其他監督信息的原因。

3.4 結果分析

由表5中可以看到，流形嵌入分箱方法占用內存情況。在StrainMock數據中的內存占用不是最少的，但較于SolidBin來說，本文方法內存占用與COCACOLA和MaxBin方法同在一個數量級。在SpeciesMock數據集，本文方法占用的內存是最少的，且橫向對比不同方法在兩個數據集上的內存消耗情況，本文方法的增長幅度不大，這說明本文方法具有很好的可拓展性。

表5 不同方法內存使用 單位：MB

不同方法的運行時間如表6所示，本文方法與其他分箱方法相比，運行時間最短，在StrainMock數據集上用了不到一分鐘的時間,在SpeciesMock數據集上也僅僅用了88.21 s。另外三種分箱方法最少運行時間也是本文方法的4.95倍。SolidBin使用的點對點的約束項是造成其內存消耗較大的原因，而MaxBin方法中EM過程中的迭代是其消耗時間較長的主要原因。從結果中發現，流形嵌入作為疊連群非線性降維方法，為分類器提供了更有效的分箱特征，極大地縮短了分箱器的運行時間，提高了分箱器的計算性能。

表6 不同方法運行時間 單位：s

4 結 語

本文將流形嵌入引入到了宏基因組分箱方法中，開發了基于流形嵌入的宏基因組疊連群分箱工具，基于數據流形嵌入，使用變分貝葉斯混合高斯模型可以更高效方便地估算出初始化分箱數，相比較之前使用基于單拷貝基因的方法，所需運行時間較少，運行效率更高，同時具有良好的可拓展性，可運行在個人電腦上。本文提出的基于流形分箱的方法，相比較MaxBin、COCACOLA、SolidBin，在SpeciesMock擁有更高的ACC、NMI、ARI。本方法尤其適用于菌種水平上、多樣本測序的數據。本方法的不足之處是對于樣本較少的數據效果較差，即對疊連群的覆蓋度信息敏感。在之后的研究里，我們將融入其他先驗信息，比如讀長連接信息到流形分箱方法中，進一步地提高本文方法的性能。