關(guān)聯(lián)性特征在宏基因組分裝中的應(yīng)用*

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

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

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

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

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

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

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

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

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

1.1 堿基對(duì)關(guān)聯(lián)性

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

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

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

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

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

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

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

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

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

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

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

1.3 機(jī)器學(xué)習(xí)算法

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

本文之所以選擇K均值分類算法是因?yàn)槭紫冗@種算法分別是無監(jiān)督方法，其次因?yàn)榕c作者提取的特征值進(jìn)行對(duì)比的MetaCluster3.0軟件里使用的分類算法就是K均值分類法，這樣使用同一種分類法進(jìn)行測(cè)試效果才有公平性。

2 實(shí)驗(yàn)結(jié)果及其討論

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

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

2.2 分類效果及其分析

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

2.3 K均值分類效果對(duì)比及分析

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

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

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

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

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

3 總結(jié)

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