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

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

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

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

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

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

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

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

1 基于序列關聯(lián)性特征的分裝方法

1.1 堿基對關聯(lián)性

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

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

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

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

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

1.2 三聯(lián)體核苷酸關聯(lián)性

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

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

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

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

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

1.3 機器學習算法

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

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

2 實驗結果及其討論

2.1 模擬宏基因組測序數(shù)據(jù)集

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

2.2 分類效果及其分析

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

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

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

表1 選取進行性能對比的模擬數(shù)據(jù)集

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

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

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

3 總結

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

［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，F(xiàn)arris 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］孫嘯，傅靜，焦典，等.生物信息學若干前沿問題的探討:利用序列統(tǒng)計特征分析基因組序列［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.