劉爽爽， 帖 云， 齊 林， 劉峰輝， 王 磊

(1. 鄭州大學(xué) 信息工程學(xué)院 河南 鄭州 450001； 2. 鄭州大學(xué)第一附屬醫(yī)院 河南 鄭州 450052； 3. 河南省人民醫(yī)院 口腔醫(yī)學(xué)中心 河南 鄭州 450003)

0 引言

16S rRNA基因是微生物生態(tài)學(xué)分析中最常使用的一類分子標志物，存在于所有細菌的基因組中，具有高度的特異性和保守性，序列長度約為1 500 bp[1]。除保守區(qū)外，16S rRNA基因序列還存在V1～V9共9個可變區(qū)[2]，不同可變區(qū)的長度范圍為100～300 bp，新一代測序技術(shù)可以使用短配對堿基輕易覆蓋，使得 16S rRNA 序列可變區(qū)的測量更加便捷。

可變區(qū)的特異性能夠反映出不同微生物的特征核苷酸序列，用于分析復(fù)雜生物環(huán)境中微生物的物種多樣性[3]、相對豐度[4]、物種鑒定及進化距離等[5]。文獻[6]對16S rRNA 基因兩個可變區(qū)進行了比較研究，結(jié)果顯示，在腸道菌群物種多樣性分析及物種鑒定能力上，選擇V1～V3可變區(qū)片段進行測序，得到了與全長序列更為接近的結(jié)果。文獻[7]利用16S rRNA全長序列和部分基因作為熱測序的靶點，在不同水平上分析了16S rRNA基因在基因組內(nèi)因異質(zhì)性引起的高估問題，結(jié)果表明，對于細菌使用針對V4和V5可變區(qū)的引物可以將這種高估最小化。文獻[8]采用Illumina Miseq測序技術(shù)，測定了蘇尼特和阿拉善雙峰駝的自然發(fā)酵駝乳中微生物16S rRNA的V3、V4可變區(qū)序列，并對群落結(jié)構(gòu)和物種多樣性進行了比較分析。文獻[9]對HIV-1包膜蛋白gp120進行分析時，找到了可變區(qū)V1可能作為傳播選擇靶位點的證據(jù)。文獻[10]分析了大西洋鮭魚細菌16s rRNA基因全長序列及不同可變區(qū)對微生物群落結(jié)構(gòu)的影響，發(fā)現(xiàn)不同可變區(qū)對微生物分布和系統(tǒng)發(fā)育有著不同的影響。目前可變區(qū)在物種進化關(guān)系中表現(xiàn)如何的研究較少，但其對物種進化來源分析具有重要的指導(dǎo)意義。本文以核糖體數(shù)據(jù)庫項目(RDP)所提供的細菌16S rRNA基因數(shù)據(jù)為基礎(chǔ)，構(gòu)建不同可變區(qū)及全長序列進化樹，使用層次距離矩陣算法分析了V2、V3、V4可變區(qū)與全長序列所構(gòu)建的進化樹之間的距離差異值，并對可變區(qū)與全長序列進化關(guān)系的相似性進行了分析。

1 數(shù)據(jù)獲取與預(yù)處理

1.1 可變區(qū)截取與篩選

原始數(shù)據(jù)采用RDP中細菌16S rRNA的全部序列。壓縮文件大小為3 GB,解壓后大小為76 GB，共包含約320萬條16S rRNA序列，數(shù)據(jù)格式為fasta。由于V9區(qū)在實際研究中應(yīng)用較少，故只選用V1～V8可變區(qū)進行相關(guān)研究。使用MEGA6軟件在該序列中分別尋找各可變區(qū)兩側(cè)保守序列，保守序列及可變區(qū)位置如表1所示。

表1 可變區(qū)兩側(cè)保守序列及位置

確定了可變區(qū)位置后，使用biopython函數(shù)庫(https:∥biopython.org/)中的序列切片方法進行可變區(qū)片段的截取。初步截取后的序列中仍存在一些含有實際堿基數(shù)目較少且信息量較低的序列，因此需要分析可變區(qū)片段實際堿基長度并篩選出含有一定信息量的序列。使用biopython庫中的seq.parse函數(shù)讀取初步處理后的序列，統(tǒng)計每個序列中的實際堿基數(shù)目，并使用matplotlib庫(https:∥matplotlib.org/)繪制可變區(qū)各序列的實際堿基數(shù)目分析圖，以便對序列進行初步篩選。對序列的初步處理操作由python腳本完成，所使用的核心函數(shù)庫為biopython庫。

8個可變區(qū)使用相同的方法進行實際堿基數(shù)目統(tǒng)計，以V2可變區(qū)為例，實際堿基數(shù)目結(jié)果如圖1所示?？梢钥闯?，可變區(qū)片段中實際堿基數(shù)目出現(xiàn)了明顯的拐點，有約70萬個V2可變區(qū)片段中堿基缺失較為嚴重。將拐點處放大，可以觀察到部分V2可變區(qū)片段實際堿基數(shù)目在80以下，表明在這些序列的測序過程中，V2可變區(qū)的測序出現(xiàn)了遺漏或者并未對V2可變區(qū)進行測序。因此，需要按照拐點處實際堿基數(shù)目對約300萬個可變區(qū)片段進行篩選，以保留含有一定信息量的可變區(qū)片段。

圖1 V2可變區(qū)實際堿基數(shù)目

篩選操作仍由python腳本完成，在完成了8個可變區(qū)片段的初步提取和篩選后，就得到了可使用命令行工具進行處理的數(shù)據(jù)。對提取出的可變區(qū)片段進行去冗余與去除嵌合體操作，以V2可變區(qū)為例，兩端對齊可以發(fā)現(xiàn)，序列數(shù)目為2 487 500，片段長度為1 264 bp，實際堿基數(shù)目為80～120，分別使用unique.seqs與chimera.uchime函數(shù)去除冗余部分序列和包含嵌合體較多的序列，再次對序列進行總結(jié)分析，此時序列數(shù)目縮減為533 136，片段長度縮減為771 bp。在完成了所有可變區(qū)的篩選、過濾操作后，繪制8個可變區(qū)預(yù)處理前后序列堿基數(shù)目對比圖，結(jié)果如圖2所示。通過可變區(qū)截取、篩選等數(shù)據(jù)預(yù)處理后，從序列堿基數(shù)目對比圖可以看出，V2、V3、V4可變區(qū)預(yù)處理后序列堿基數(shù)目較其他可變區(qū)多，且V2、V4可變區(qū)序列數(shù)目比V3可變區(qū)多，V2、V4兩個可變區(qū)較其他可變區(qū)包含更多的序列信息。

圖2 可變區(qū)預(yù)處理前后序列堿基數(shù)目對比

1.2 可變區(qū)OTU聚類

為了方便區(qū)分序列，不同的16S rRNA基因序列若相似性高于97%，就可以把它定義為一個操作分類單元(operational taxonomic unit，OTU)，每個OTU對應(yīng)于一個不同的微生物種。通過OTU聚類分析可以簡化數(shù)據(jù)結(jié)構(gòu)，得到樣品中微生物多樣性水平以及不同微生物的豐度。對各可變區(qū)分別按相似度97%、98%及100%進行OTU聚類，取均值作為各可變區(qū)OTU數(shù)目，對比結(jié)果如圖3所示。通過可變區(qū)聚類后的OTU數(shù)目可以看出，V2、V4兩個可變區(qū)在多樣性水平上較其他可變區(qū)更接近全長序列。

圖3 可變區(qū)聚類后OTU數(shù)目對比

通過以上分析可以看出，V2、V4可變區(qū)在含有獨特序列數(shù)目和實際堿基長度上均優(yōu)于其他6個可變區(qū)，使用包含V2、V4可變區(qū)的基因測序片段對序列進行分類，將得到更加接近全長序列的分類結(jié)果。

2 進化關(guān)系相似性分析

2.1 進化樹的構(gòu)建

在各個可變區(qū)分析結(jié)果的基礎(chǔ)上，為了更好地反映可變區(qū)對物種進化關(guān)系的相似性且盡可能減少計算的復(fù)雜度，選用的數(shù)據(jù)必須在門、綱、目、科、屬、種層次上具有良好的區(qū)分度。使用RDP網(wǎng)站的Browser在數(shù)據(jù)庫中選取了56條序列，這些序列來自11種不同門、16種不同綱的56個不同屬種的細菌。為了將V2、V4可變區(qū)與其他可變區(qū)進行對比分析，基于OTU聚類的結(jié)果，選取V3可變區(qū)作為對照組，分別對V2、V3、V4可變區(qū)構(gòu)建進化樹，并與全長序列進化樹進行比較。為了準確截取出V2、V3、V4可變區(qū)，將一條細菌序列添加到多重比對中，使用MEGA軟件搜索特定保守位點序列，重新比對后3個可變區(qū)位點為V2(203～397)、V3(583～657)、V4(735～855)。截取出可變區(qū)后進行實際堿基數(shù)目分析，發(fā)現(xiàn)在V2、V3、V4可變區(qū)序列中S000583665、S000830684、S000346245、S000120585這4條序列中存在實際堿基數(shù)目較少的可變區(qū)片段。為了保證使用完全一致的物種類別構(gòu)建進化樹，只保留3個可變區(qū)中實際堿基數(shù)目均較高的序列，最終得到可用于構(gòu)建進化樹的V2、V3、V4可變區(qū)以及全長序列數(shù)據(jù)。

將用于進化樹構(gòu)建的數(shù)據(jù)在ClustalX中重新比對，比對完成后將ClustalX生成的.dnd文件使用TreeView軟件打開，構(gòu)建V2、V3、V4可變區(qū)及全長序列進化樹，結(jié)果如圖4～圖7所示。

圖4 V2可變區(qū)進化樹

圖5 V3可變區(qū)進化樹

圖6 V4可變區(qū)進化樹

圖7 全長序列進化樹

結(jié)果顯示，V3可變區(qū)片段生成的進化樹與全長序列生成的進化樹有著較大的偏差，而V2、V4可變區(qū)在進化樹結(jié)構(gòu)上與全長序列很相似。對于進化距離較小的序列對，使用V2、V4可變區(qū)構(gòu)建進化樹仍能得到比較接近全長序列的進化關(guān)系，例如S000543677與S000587182、S000691981與S000946165、S000007759與S000649409。以S000543677與S000587182為例，在4棵進化樹中這兩條序列之間的進化距離都非常小，將這兩個序列視為一個結(jié)點，可以看出，在全長序列中與該結(jié)點進化距離最小的序列為S000345627，這三者的距離關(guān)系與在V2、V4可變區(qū)進化樹中的距離關(guān)系是吻合的，而與V3可變區(qū)構(gòu)建進化樹的結(jié)果相距甚遠。這說明3個可變區(qū)在進化關(guān)系上與全長序列之間都存在一定的差異，但V2、V4兩個可變區(qū)所構(gòu)建的進化樹在可信度上要略優(yōu)于V3可變區(qū)。

2.2 算法描述

為了定量分析可變區(qū)進化樹與全長序列進化樹之間的相似程度，需要建立一個能夠衡量兩棵進化樹之間相似度的方法。因此，本文提出一種能評價相同序列不同片段構(gòu)成的不同進化樹之間進化關(guān)系相似度的方法，即對任意一棵進化樹，都可以建立一個層次距離矩陣，通過對兩棵樹層次距離矩陣的比較，可以得到它們之間的距離差異值，進而可以分析進化關(guān)系的相似程度。層次距離矩陣算法示意圖如圖8所示。對于樹中任意兩個葉子結(jié)點，兩兩計算層次距離，結(jié)點間的層次距離定義為：若兩個結(jié)點有同一父結(jié)點，則兩個結(jié)點間層次距離為0；否則，層次距離為兩個結(jié)點向上到達第一個共同祖先結(jié)點的距離權(quán)值。

圖8 層次距離矩陣算法示意圖

若M、N為樹中兩個葉子結(jié)點，則其層次距離可以表示為

(1)

其中：LM和LN分別表示結(jié)點M和N下方路徑所在層數(shù)；LMN表示結(jié)點M、N最近公共祖先結(jié)點下方路徑所在層數(shù)；W表示路徑上方結(jié)點的距離權(quán)值。

A、D兩個葉子結(jié)點之間的最近公共祖先結(jié)點為根結(jié)點，若各層的距離權(quán)值如圖8所示，則A、D兩個結(jié)點間的層次距離為根結(jié)點的距離權(quán)值3。在完成了兩棵樹的層次距離計算后，得到任意兩個序列在兩棵樹中層次距離的差異值，可以表示為

(2)

將兩個層次距離矩陣中所有對應(yīng)位置的差異值相加，則兩棵樹之間的進化關(guān)系相似度指標可以表示為

(3)

其中：S為所有葉子結(jié)點的集合。將可變區(qū)的層次距離矩陣與全長序列層次距離矩陣按位置相減并取絕對值，可以得到差異值矩陣，該矩陣可以較準確地反映任意兩個序列在不同片段構(gòu)建的進化樹的差異程度。相似度指標D能一定程度地反映兩棵進化樹在整體上的層次距離差異。D值越大，兩棵樹中對應(yīng)的序列對之間層次距離的差異越大，表明兩棵樹之間的相似度越小。

2.3 算法實現(xiàn)與結(jié)論

對于圖4～圖7中所示的進化樹，將各葉子結(jié)點序列名稱以數(shù)字0～52進行編號，并以大小寫英文字母A～Z、a～z為中間結(jié)點編號，可得到各序列在不同樹中從根結(jié)點出發(fā)的路徑。在兩兩路徑之間進行層次距離矩陣計算時，首先倒序?qū)ふ覂蓚€序列中第一個相同的中間結(jié)點，接著提取該中間結(jié)點的層次權(quán)值，得到兩路徑之間的層次距離。利用式(1)～(3)進行差異值矩陣計算，可得V2、V3、V4可變區(qū)與全長序列之間的距離差異值分別為59 052、87 154和45 848，相似度分別為41%、13%和55%?？梢钥闯觯琕4可變區(qū)與全長序列之間兩兩序列的距離差異值要小于V2、V3可變區(qū)，且V4可變區(qū)與全長序列之間進化關(guān)系的相似度為55%，大于V2、V3可變區(qū)與全長序列之間進化關(guān)系的相似度。因此，V4可變區(qū)所構(gòu)建的進化樹更接近全長序列所構(gòu)建的進化樹，在可信度上要優(yōu)于V2、V3可變區(qū)，在細菌物種進化關(guān)系上V4可變區(qū)較V2、V3可變區(qū)片段更為接近全長序列。

2.4 算法比較

與傳統(tǒng)的距離與相似度算法如歐氏距離算法、馬氏距離算法、漢明距離算法相比，層次距離矩陣算法的優(yōu)點在于若兩結(jié)點都是根結(jié)點的子結(jié)點，在不使用距離權(quán)值的情況下，兩結(jié)點的層次距離是很小的，這與實際情況不符，而引入距離權(quán)值之后則可以解決這一問題，即在樹中越接近頂端的結(jié)點與相同層次的結(jié)點之間的序列距離越大，兩結(jié)點的層次距離差距也越大。分別使用層次距離矩陣算法、歐氏距離算法、馬氏距離算法、漢明距離算法計算V2、V3、V4可變區(qū)與全長序列之間的距離差異值，結(jié)果如表2所示。這四種算法的均方值分別為21 096.02、29 887.54、28 318.72和34 960.29，精度分別為79.0%、70.2%、71.7%和65.1%?？梢钥闯?，層次距離矩陣算法與傳統(tǒng)的距離與相似度算法相比，在各可變區(qū)與全長序列距離上的均方值更小，計算距離的精度更高，使用這一算法可以計算樹中任意兩個葉子結(jié)點之間的層次距離，并得到不同樹之間所有序列對在層次距離上的差異值，而引入距離權(quán)值主要是基于對距離偏移大小的考慮。以上分析對物種進化來源分析研究具有重要的指導(dǎo)意義，若要研究某一菌種的突變來源，可利用可變區(qū)分別構(gòu)建進化樹并比較相似度的方法進行分析。

表2 不同算法計算的可變區(qū)與全長序列之間的距離差異值

3 討論

本文以RDP所提供的細菌16S rRNA數(shù)據(jù)為基礎(chǔ)，對不同可變區(qū)物種進化關(guān)系的相似性進行了研究，分別對這些可變區(qū)片段進行重新比對和構(gòu)建進化樹，并且與全長序列構(gòu)建的進化樹進行比較，發(fā)現(xiàn)V4可變區(qū)在進化關(guān)系上與全長序列更為貼近，V4可變區(qū)構(gòu)建進化樹的可信度要優(yōu)于V2、V3可變區(qū)。本文使用的層次距離矩陣算法對物種進化來源分析研究有一定的指導(dǎo)意義，較傳統(tǒng)的距離與相似度算法具有更好的性能。但是，這并不能否認單獨利用可變區(qū)進行物種進化關(guān)系分析存在一定的局限性。另外，本文所提出的相似度計算方法雖然能在一定程度上反映兩棵進化樹之間的相似程度，但也存在一些不足之處。在提出層次距離矩陣計算方法之前，曾嘗試使用樹的層次遍歷進行轉(zhuǎn)換代價計算，在葉子結(jié)點相同的情況下，這種方法可以將一棵樹轉(zhuǎn)換成另一棵樹的形式，通過計算這一轉(zhuǎn)換過程中的代價，可以評價兩棵樹的相似程度，但是該方法的問題在于兩棵樹的非葉結(jié)點數(shù)目可能不同，如果能夠解決這一問題，那么使用轉(zhuǎn)移代價來評價樹的相似程度是更為可行的一種方法。