關聯性特征在宏基因組分裝中的應用*

張倩倩，曹唱唱，丁 嘯，孫 嘯

(東南大學生物科學與醫學工程學院，南京210096)

宏基因組是指環境中全部微生物DNA的總和，最早是由Handelsman等人于98年在一篇研究土壤微生物的文章中提出［1］。宏基因組學，直接對混合的微生物群落樣本進行基因組提取，然后以多種微生物基因組的混合體為測序模版，對其進行高通量測序。傳統的研究方法是對單一微生物進行單純培養，然后再對其進行分離測序研究，這樣的方法停留在單一微生物物種的水平上。由于環境中99%的微生物都難以用常規的方法進行培養，所以傳統的方法在很大程度上受到了限制。

然而，隨著高通量測序技術的發展，當環境中的微生物樣本被測序之后，需要對宏基因組測序片段建立種族進化關系［2］，這樣一個步驟叫做分裝。分裝也是一種特殊的聚類方式，指對宏基因組測序片段進行重疊區域保守拼接后根據一定的規則進行聚類，并將其歸類到已知的系統發生關系中。歸類的物種層次精度不一樣，精確的可以歸類至種，而粗糙的只能歸類到界、門。歸類的精度取決于多個因素如分裝方法、群落結構和測序質量及深度［3］。微生物群落結構由研究樣本所決定，測序質量和深度則依靠測序技術和樣本尺寸。所以對于分裝主要的研究重點就是分裝方法。目前存在的分裝方法分為2類:一類是基于序列相似性比對的分裝方法;另一類是基于序列特征的分裝方法。

基于序列相似性比對的分裝方法其主要步驟是將宏基因組DNA片段進行處理之后與已知微生物基因組數據庫做比對，然后再用比對所獲得的信息對輸入的測序片段進行分類和標注。最早的此類分裝方法是由07年德國蒂賓根大學的Huson等人開發的MEGAN軟件［4］。其原理直截了當，只需要用戶利用序列比對工具BLAST將宏基因組測序數據與NCBI的微生物基因組數據庫進行比對，生成一個軟件可讀取的比對中間文件。最后通過計算機分析，基于最近公共祖先算法將宏基因組DNA序列歸類到物種分類樹的相關結點上，并輸出結果。但是這類算法太依賴數據庫，局限性很大，畢竟我們已知的參考基因組數量有限，而且只有不到1%的微生物基因組可以培養和測序得到。于是研究者們開始分析并利用序列本身固有的成分特征來進行分裝，從而提出了基于序列成分特征的分裝方法。

基于序列特征的分裝方法又包含2類:一類是基于物種分類標志物基因;另一類是基于堿基子串使用偏向性。物種分類標志物基因是指某一微生物在編碼區域所特有的一段DNA序列，用于區分不同物種的標志，如傳統的16S rRNA［5］，以及新加入的rec A和rpo B基因等。利用這類特征進行分裝具有以下局限性:首先，不是所有的宏基因組測序得到的DNA序列中都包含標志物基因的;其次，近年來科學家也發現，標志物基因也存在著種間平行轉移的情況;最后，這種方法仍然還是依賴于物種分類標志物基因數據庫，極大地限制了其應用范圍。

針對以上限制，有學者提出使用基因組序列特征進行分裝。目前主流的分裝方法都是利用序列中堿基子串的組成特征，這種方法的理論依據是由Karlin等人在上世紀90年代中期對多種微生物全基因組序列進行相關研究后得出的［6-7］。他們發現同一物種基因組的不同片段有著非常相似的堿基子串使用偏向性，而不同基因組之間使用偏向性就很大?；谶@樣一種特性近年來也有很多新穎的分裝方法，例如本文中用來做性能對比的由YANG等人提出的分裝軟件MetaCluster［8］，還有Chan等人提出的基于自組織生長圖(GSOM)算法［9-10］等。

基因組序列特征可分為序列組成性特征和關聯性特征，而目前分裝算法使用的都是組成性特征，如果能夠在分裝算法中引入關聯性特征，則可望提高分裝算法的性能。本文主要分析序列關聯性特征在宏基因組分裝上的應用。

1 基于序列關聯性特征的分裝方法

1.1 堿基對關聯性

BBC特征反映的是堿基對的關聯性(Base-Base Correlation)［11］，該特征由序列的互信息 MIF(Mutual Information Function)定義而來。序列的互信息計算公式如下:

其中Pi表示單個核苷酸ni∈{A，G，C，T}出現的頻率，Pij(k)表示一對被k個核苷酸分隔的核苷酸ni和nj出現的頻率。這樣I(k)表示，當識別到核苷酸X，得到相距k個核苷酸的核苷酸為Y時產生的信息量(以比特為單位)。舉個例子來直觀地說明I(k)的含義。

例子1 考慮一條隨機序列，組成序列的各個核苷酸是獨立無關的。直觀上就可以看出，我們不能從X中得到任何Y的信息，所以對于任意k，I(k)都為零。事實上，從式(1)也可以得到相同的結果。由于所有的核苷酸統計上相互獨立，所以由統計學公式可以得到:對所有的i、j和k，Pij(k)=PiPj。把這個式子代入到式(1)中，就可以得到對數中的因子為1，所以I(k)就為零。

由此我們可以定義序列的BBC特征如下:

其中Pi和Pij(l)的含義同上。Tij(k)表示不同間隔的二核苷酸組合在k+1長度上的平均相關性，反映了核苷酸序列的一種局部特征。

1.2 三聯體核苷酸關聯性

三聯體核苷酸關聯性特征，即序列中每一個三聯體核苷酸單堿基與堿基對之間的關聯性的量化。其定義為:宏基因組測序片段中三聯體核苷酸第一位堿基與后兩位二聯體核苷酸的關聯特性參數。這種特征一共包含了68個特征參數，分為2個部分。

第1部分是在確定第一位堿基的條件下后兩位二聯體核苷酸的出現頻率，一共包含64個特征參數。這部分特征值的計算公式為:

其中nij表示確定第一位堿基后，三聯體核苷酸的后兩位核苷酸取16種不同堿基對時在序列中分別出現的次數表示確定第一位堿基的所有三聯體核苷酸在序列中出現的次數。這部分特征值體現的是當第一位堿基確定之后，后兩位核苷酸作為一個整體使用的頻度。

第2部分是三聯體核苷酸中的4種第一位堿基與后兩位二聯體核苷酸的相互信息量，一共包含4個特征參數。計算公式為:

其中Pi表示第一位堿基在序列中出現4種類型的頻率，Pj表示三聯體核苷酸后兩位核苷酸取特定堿基對時在序列中出現的頻率，Pij表示指定三聯體核苷酸在序列中出現的頻率。1～4和1，2，…，16分別對應單堿基A、T、C、G和16種堿基對的編號。

1.3 機器學習算法

本文采用K均值分類法，K均值分類法原理是給定一個數據點集合和需要的聚類數目K，K均值算法根據某個距離函數反復把數據分入K個聚類中。其算法步驟為先隨機選取K個對象作為初始的聚類中心。然后計算每個對象與各個種子聚類中心之間的距離，把每個對象分配給距離它最近的聚類中心。一旦全部對象都被分配了，每個聚類的聚類中心會根據聚類中現有的對象被重新計算。這個過程將不斷重復直到滿足某個終止條件。K均值分類法是一種無監督分類法，并以其簡潔和效率而廣泛使用。本文使用的K均值分類軟件為Gene Cluster3.0［12］。

本文之所以選擇K均值分類算法是因為首先這種算法分別是無監督方法，其次因為與作者提取的特征值進行對比的MetaCluster3.0軟件里使用的分類算法就是K均值分類法，這樣使用同一種分類法進行測試效果才有公平性。

2 實驗結果及其討論

2.1 模擬宏基因組測序數據集

從IMG［13］系統中根據分類樹隨機提取了18種微生物，其中包括2種古生菌、4種真核微生物、12種細菌。在此基礎上根據不同的宏基因組數據復雜度創建了116組模擬數據集。首先根據分類學物種差異層次，數據集被分為同門異綱、同綱異目、同目異科、同科異屬、同屬異種5大組;然后我們在每個大組中，根據不同的基因組測序片段reads長度、測序錯誤率、相對豐度再進行模擬。當模擬數據集中的測序片段reads長度分別取500 bp、1 000 bp、2 000 bp、3 000 bp時，默認的測序錯誤率為1%，相對豐度為1∶1;當模擬數據集中的測序錯誤率分別取2%、3%、5%時，默認的 Reads長度為2 000 bp，相對豐度為1∶1;當模擬數據集中的相對豐度分別取1∶2、1 ∶4、1 ∶8時，默認的 Reads長度為 2 000 bp，測序錯誤率為1%。其中不同大組的數據集也會包含不同數量的微生物，其中可能有2種、3種、4種、8種或者10種微生物。

2.2 分類效果及其分析

對模擬數據集的分類結果分為2部分:第1部分是分別利用K均值分類法加入三聯體核苷酸關聯性特征參數和MetaCluster3.0軟件對模擬數據集進行分類;第2部分是利用SVM分別加入三聯體核苷酸關聯性特征參數、三聯體核苷酸使用頻率和四聯體核苷酸使用頻率對模擬數據集進行分類。本文評估分類性能是根據分類準確率來判斷，公式如下:

2.3 K均值分類效果對比及分析

由于模擬數據集比較多，限于篇幅作者只選取了一些代表性的數據集分類效果進行展示。選取的數據集和分類效果分別列在表1和表2。

表1 選取進行性能對比的模擬數據集

表2 MetaCluster3.0與堿基對關聯性以及三聯體核苷酸關聯性方法的分裝準確率對比

通過分析以上2張表格，我們可以發現在一樣使用K均值分類方法時，MetaCluster3.0軟件分類只有在同綱異目的物種層次時才能對模擬數據集進行較好的分裝，在同科異屬的包含3種微生物的數據集中表現出色，可能是由于數據集選擇的微生物物種基因組中四聯核苷酸使用頻率恰好很相似的原因，因為作者通過對所有數據集分裝后發現這是一個特例。通過數據也可以發現這款分裝軟件對宏基因組數據集的分類層次很敏感。反觀三聯體核苷酸關聯性特征參數，分類效果多數情況下都要優于單純使用四聯核苷酸使用頻率作為特征參數的Meta-Cluster3.0軟件。

但是通過數據分析我們也發現，利用三聯體核苷酸關聯性特征參數進行分類存在兩個問題:在對只包含真核微生物的數據集(數據集6-1-3，11-1-2)和豐度比不均勻的模擬數據集進行分類時效果不如意(數據集3-3-2)。作者認為可能因為真核微生物基因組較為復雜，其中包括大量非編碼序列，而原核微生物基因組基本都是編碼序列。而密碼子就是由三位堿基組成，所以三聯體核苷酸關聯性特征參數對真核微生物的分類效果不好。對于豐度比不均勻導致分類效果不好，作者認為可能是由于使用的單純的使用K均值分類法太過于簡單，對豐度比比較敏感。

3 總結

得益于高通量測序技術的日漸成熟，宏基因組學在微生物群落的研究上進展迅速，但是同時伴隨著一個問題。如何將大量的DNA短片段快速正確的分開，從而可以對宏基因組中的各種微生物進行研究，包括發現新物種、鑒定新基因及其發現微生物群的新功能等。宏基因組分裝技術就是解決這個問題的核心，目前分裝技術的主流是基于堿基使用偏向性特征的無監督算法。所以如何提取一種有效而又簡單的序列特征參數是提升分裝技術的重點。本文提出了一種三聯體核苷酸關聯性特征參數，它并不是單純的使用三聯體核苷酸使用頻率，它體現的是序列中每一個三聯體核苷酸單堿基與堿基對之間的關聯性的量化。本文通過使用K均值分類法和支持向量機對模擬的不同復雜度的宏基因組數據集進行分類，發現利用三聯體核苷酸關聯性特征效果對大多數的數據集分類效果要優于無監督分裝算法MetaCluster3.0，同時也好于單純的利用三聯、四聯體核苷酸使用頻率的分類效果。特別是在對種的分類中保持較高的準確率。但是這種特征值仍然存在問題，在使用無監督算法進行分類時，無法有效分類包含真核微生物和豐度比不均勻的模擬數據集。作者推測是由于三聯體是密碼子的結構，而真核微生物的基因組非編碼序列要遠遠多于原核微生物基因組，導致三聯核苷酸關聯性特征失效。針對這個問題我們下一步準備基于四聯核苷酸的關聯性特征來設計分裝算法。對于無法分類豐度比不均勻數據集的問題，我們準備嘗試不同的機器學習算法或者對K均值分類法的距離進行改進，進一步的提高基于堿基關聯性特征分裝方法的準確性。

［1］Handelsman J，Rondon M R，Brady S F，et al.Molecular Biological Access to the Chemistry of Unknown Soil Microbes:A New Frontier for Natural Products［J］.Chemistry and Biology，1998，5(10):R245-R249.

［2］McHardy A.Rigoutsos I.What’s in the Mix:Phylogenetic Classifcation of Metagenome Sequence Samples［J］.Current Opinion in Microbiology 2007，10:499-503.

［3］Mavromatis K，Ivanova N，Barry K，et al.Use of Simulateddata Sets to Evaluate the Fidelity of Metagenomic Processing Methods［J］.Nature Method，2007，4(6):495-500.

［4］Huson D H，Auch A F，Qi J，et al:MEGAN Analysis of Metagenomic Data［J］.Genome Res，2007，17(3):377-386.

［5］Cole J R，Chai B，Farris R J，et al.The Ribosomal Database Project(RDP-Ⅱ):Sequences and Tools for High-Throughput rRNA A-nalysis［J］.Nucleic Acids Res，2005，33(Database issue):D294-296.

［6］Karlin S，Ladunga I.Comparisons ofEukaryotic Genomic Sequences［J］.Proc Natl Acad Sci USA，1994，91(26):12832-12836.

［7］Karlin S，Burge C.Dinucleotide Relative Abundance Extremes:A Genomic Signature［J］.Trends Genet，1995，11(7):283-290.

［8］Yang B，Peng Y，Leung H C M，et al.Unsupervised Binning of Environmental GenomicFragmentsBased on an ErrorRobust Selection of l-Mers［J］.BMC Bioinformatics，2010，11:1471.

［9］Chan C K，Hsu A L，Tang S L，et al.Using Growing Self-Organising Maps to Improve the Binning Process in Environmental Whole-Genome Shotgun Sequencing［J］.J Biomed Biotechnol，2008，2008:513701.

［10］Chan C K，Hsu A L，Halgamuge S K，et al.Binning Sequences U-sing very Sparse Labels within a Metagenome［J］.BMC Bioinformatics，2008，9:215.

［11］孫嘯，傅靜，焦典，等.生物信息學若干前沿問題的探討:利用序列統計特征分析基因組序列［M］.中國科學技術大學出版社，2004:51-56.

［12］Yang B，Peng Y，Leung H C M，et al.Unsupervised Binning of Environmental GenomicFragmentsBased on an ErrorRobust Selection of l-Mers［C］//BMC Bioinformatics 2010，11:1471.Liu Y Z，Moll J L，Spicer W E.Appl Phys Lett，1970，17(2):60-62.

［13］Markowitz V M，Chen I M A，Palaniappan K，et al.IMG:the Integrated Microbial Genomes Database and Comparative Analysis System［J］.Nucleic Acids Research，2012，40:115-122.