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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

1.2 可變區(qū)OTU聚類(lèi)

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

圖3 可變區(qū)聚類(lèi)后OTU數(shù)目對(duì)比

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

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

2.1 進(jìn)化樹(shù)的構(gòu)建

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

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

圖4 V2可變區(qū)進(jìn)化樹(shù)

圖5 V3可變區(qū)進(jìn)化樹(shù)

圖6 V4可變區(qū)進(jìn)化樹(shù)

圖7 全長(zhǎng)序列進(jìn)化樹(shù)

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

2.2 算法描述

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

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

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

(1)

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

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

(2)

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

(3)

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

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

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

2.4 算法比較

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

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

3 討論

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