丁 可，閻星羽，廉振穎，成麗娟，刁玉濤，王家林

（1.山東第一醫(yī)科大學(xué)，山東 泰安 271016；2.山東第一醫(yī)科大學(xué)基礎(chǔ)醫(yī)學(xué)院/基礎(chǔ)醫(yī)學(xué)研究所，山東 濟(jì)南 250062；3.山東中醫(yī)藥大學(xué)第二附屬醫(yī)院，山東 濟(jì)南 250001；4.山東第一醫(yī)科大學(xué)附屬腫瘤醫(yī)院，山東 濟(jì)南 250117）

細(xì)菌微生物16S rRNA 基因擴(kuò)增子二代測(cè)序數(shù)據(jù)分析的后續(xù)任務(wù)之一是通過(guò)系統(tǒng)樹(shù)形圖（dendrogram/phylogenetic tree）展示微生物種群之間在進(jìn)化或種系發(fā)生上的相互關(guān)系情況，dendrogram/phylogenetic tree 圖分3 種類型，分別為進(jìn)化分支圖（cladogram），僅有拓?fù)浣Y(jié)構(gòu)，不能從數(shù)量上說(shuō)明各分支間進(jìn)化距離的大??；系統(tǒng)發(fā)生（發(fā)育）圖（phylogram），各分枝長(zhǎng)度表示堿基替換數(shù)，因而能從數(shù)量上說(shuō)明各分支間進(jìn)化距離的大?。粫r(shí)序圖（chronogram）則用各分枝長(zhǎng)度表示進(jìn)化時(shí)間，其中系統(tǒng)發(fā)生圖（phylogram）無(wú)疑是微生物分子遺傳學(xué)研究最常用到的圖形化方法[1]。目前系統(tǒng)發(fā)生圖制作軟件生成的系統(tǒng)發(fā)生樹(shù)文件為Newick 和Nexus 格式，并且大多數(shù)系統(tǒng)發(fā)生樹(shù)查看軟件也主要兼容上述2 種文件格式，但Newick 和Nexus 格式的樹(shù)文件只包含物種間的遺傳距離和進(jìn)化拓?fù)潢P(guān)系，結(jié)構(gòu)比較單一。而R語(yǔ)言中的“ggtree”包對(duì)ape、ggplot2 等繪圖包的功能進(jìn)行了優(yōu)化，“ggtree”除了能支持眾多的樹(shù)文件格式，其繪圖方法更豐富，繪制的圖形更準(zhǔn)確和美觀，最大的優(yōu)點(diǎn)是能通過(guò)各種附加信息對(duì)樹(shù)形圖進(jìn)行多方面的注釋[2]。因此，本研究構(gòu)建了微生物基因組擴(kuò)增子二代測(cè)序數(shù)據(jù)的系統(tǒng)發(fā)生樹(shù)形圖的制作流程，通過(guò)QIIME 命令代碼對(duì)測(cè)序數(shù)據(jù)進(jìn)行預(yù)處理，聚類生成OTUs 表數(shù)據(jù)，對(duì)OTUs 表格數(shù)據(jù)進(jìn)行篩選和轉(zhuǎn)換，并用QIIME 分析管道（pipeline）的相關(guān)命令生成Newick 格式的系統(tǒng)發(fā)生樹(shù)文件，以此文件為輸入，用R 語(yǔ)言“ggtree”包的函數(shù)和代碼繪制包含更多研究信息的不同類型的系統(tǒng)發(fā)生樹(shù)圖，現(xiàn)將結(jié)果報(bào)道如下。

1 材料與方法

1.1 數(shù)據(jù)來(lái)源 分析數(shù)據(jù)來(lái)自文獻(xiàn)所用的原始測(cè)序數(shù)據(jù)[3]，基于人體微生物在不同部位并隨時(shí)間推移的變化，通過(guò)QIIME 分析方法，只選取其中的部分?jǐn)?shù)據(jù)，即每天分別從2 個(gè)人的舌部、左手手掌、右手手掌和腸道共取得34 個(gè)標(biāo)本，在Illumina HiSeq 2000 平臺(tái)上進(jìn)行微生物16s rDNA 擴(kuò)增子測(cè)序。此外，從“https://data.qiime2.org/2017.6/tutorials/moving-pictures/emp -single -end -sequences/barcodes.fastq.gz”網(wǎng)站獲得經(jīng)過(guò)整理后的單向測(cè)序的數(shù)據(jù)。

1.2 軟件系統(tǒng) 通過(guò)QIIME 完成系統(tǒng)發(fā)育樹(shù)構(gòu)建和可視化的工作流程，并使用是R 語(yǔ)言中的“ggtree”包（3.6.0 版）進(jìn)行樹(shù)文件的可視化和注釋。

1.3 分析流程 構(gòu)建可視化系統(tǒng)發(fā)生樹(shù)的工作流程，關(guān)鍵步驟是形成BIOM（Biological Observation Matrix）格式的OTUs 表。QIIME 系統(tǒng)中有多個(gè)Python腳本用于生成OTUs 表，按聚類算法分為使用從頭聚類方法（de novo）和使用參考數(shù)據(jù)庫(kù)聚類方法。在此工作流程中，本研究使用的是基于Greengene 參考數(shù)據(jù)庫(kù)的封閉參考（Closed reference）聚類分析方法生成OTUs 表數(shù)據(jù)進(jìn)行分析的，具體步驟如下：①生成OTUs 表文件（otu_table.biom）：目的是依據(jù)上一步生成的序列文件seqs.fna，通過(guò)與最新的Green－Genes 數(shù)據(jù)庫(kù)比對(duì)進(jìn)行聚類分析，其主要輸出是OTUs 表格(out_table.biom)，表格中的內(nèi)容為記錄每個(gè)OTU 在每個(gè)樣品（微生物群落標(biāo)本）中被觀察到的次數(shù)；②OTUs 的過(guò)濾：通過(guò)filter_otus_from_otu_table.py 命令代碼篩選OTUs，只保留相對(duì)含量＞1‰（或其它比例）的OTUs；③OTUs 表文件格式轉(zhuǎn)換：目的是將上一步生成的biom 格式的OTUs 表文件轉(zhuǎn)換為純文本（csv、tsv 或txt 格式）文件，利于下一步數(shù)據(jù)處理；④選取代表序列：由于每個(gè)OTU 中的序列不完全相同，因此需要通過(guò)pick_rep_set.py 腳本選取一條代表性序列作為該OTU 的序列，用于后續(xù)分析；⑤物種分級(jí)注釋：用assign_taxonomy.py 腳本命令對(duì)每個(gè)OTU 代表性序列進(jìn)行物種分類信息的注釋，可以認(rèn)為每個(gè)OTU 近似為一個(gè)物種，反之一個(gè)分類到“種（species）”或“屬（genus）”水平的物種類別可以對(duì)應(yīng)一個(gè)到多個(gè)OTU；⑥代表序列的比對(duì)：QIIME 系統(tǒng)的align_seqs.py 腳本提供3 種序列比對(duì)方法，即Py－NAST、MUSCLE 和INFERNAL，本研究中使用QIIME 系統(tǒng)默認(rèn)的PyNAST 方法，它基于NAST 算法，將輸入序列與提供的參考序列數(shù)據(jù)比對(duì)，在數(shù)據(jù)庫(kù)中找到最高匹配的序列；且MUSCLE 不需要提供參考序列，可用于真菌轉(zhuǎn)錄間隔區(qū)（internal transcribed spacer，ITS）測(cè)序比對(duì)分析；而INFERNAL 利用RNA結(jié)構(gòu)和序列相似性進(jìn)行比對(duì)，與PyNAST 一樣需要比對(duì)數(shù)據(jù)庫(kù)；⑦篩選比對(duì)序列：由于上述align_seqs.py 腳本通過(guò)將長(zhǎng)度200～400 bp 的目的序列和16S rRNA 基因的全序列比對(duì)，因此，生成的代表性序列包含空缺（gaps），為了保留代表性序列中的有用信息以構(gòu)建系統(tǒng)發(fā)育樹(shù)，需要通過(guò)filter_alignment.py 腳本對(duì)上述代表性序列進(jìn)行篩選，去除堿基空缺等無(wú)用信息；⑧建樹(shù)：運(yùn)用python make_phylogeny.py 腳本構(gòu)建系統(tǒng)發(fā)育樹(shù)；⑨樹(shù)文件的圖形化：通過(guò)FigTree軟件和R 語(yǔ)言“ggtree”包對(duì)Newick 格式的系統(tǒng)發(fā)育樹(shù)文件進(jìn)行圖形化處理。系統(tǒng)發(fā)育樹(shù)文件生成與可視化的數(shù)據(jù)處理工作流程見(jiàn)圖1。

2 結(jié)果

本結(jié)果來(lái)自對(duì)34 個(gè)取樣標(biāo)本16S rRNA 基因測(cè)序數(shù)據(jù)的分析，34 個(gè)標(biāo)本分別來(lái)自人體的不同部位和不同取樣時(shí)間，為方便起見(jiàn)，本分析流程沒(méi)有考慮采樣時(shí)間，僅將不同取樣部位作為分組因素進(jìn)行分析。首選通過(guò)QIIME-1.9.1 系統(tǒng)的pick_closed_reference_otus.py 腳本生成OTUs 表和系統(tǒng)發(fā)生樹(shù)文件，然后根據(jù)物種分類信息，將物種名稱標(biāo)記到“樹(shù)葉”上，構(gòu)建優(yōu)化的物種聚類樹(shù)文件，該樹(shù)文件作為“ggtree”包的輸入文件，再結(jié)合研究設(shè)計(jì)的分組信息繪制不同類型的系統(tǒng)發(fā)生樹(shù)圖。

2.1 樣本信息和OTUs 聚類結(jié)果 本示例數(shù)據(jù)來(lái)自人體的腸道（gut）、左側(cè)手掌（left palm）、右側(cè)手掌（right palm）、和舌部（tongue），為了分析方便，省略了取樣時(shí)間、抗生素使用情況等其它分組信息。通過(guò)封閉參考（closed reference）聚類法將有效測(cè)序序列按照≥97%的相似性歸為一個(gè)OTU。采用RDP 算法與Greengene 16S rRNA 數(shù)據(jù)庫(kù)比對(duì)，并將各OTU注釋到所屬的分類單元。使用Greengene 數(shù)據(jù)庫(kù)可以將大多數(shù)OTU 序列注釋到“屬”水平，少部分OTU可以鑒定到“種”水平。通過(guò)上述方法將示例數(shù)據(jù)中的共177882 條序列聚類成4403 個(gè)OTU，在“門(mén)”和“屬”兩個(gè)水平上對(duì)同類的OTU 進(jìn)行分類鑒定，分別歸屬于24 個(gè)門(mén)，615 個(gè)屬，有關(guān)樣本的人體部位分組情況、OTU 聚類和種屬鑒定的具體情況見(jiàn)表1。OTUs 聚類及分類鑒定結(jié)果與同時(shí)生成的系統(tǒng)發(fā)生樹(shù)文件，可以作為后續(xù)統(tǒng)計(jì)分析，構(gòu)建系統(tǒng)發(fā)生樹(shù)圖的基礎(chǔ)。

圖1 系統(tǒng)發(fā)育樹(shù)文件生成與可視化的數(shù)據(jù)處理工作流程

表1 樣本測(cè)序序列與OTU 聚類情況統(tǒng)計(jì)表

2.2 各OTU 的物種歸類 人體不同部位的微生物構(gòu)成存在顯著差異，通過(guò)QIIME 系統(tǒng)的種群多樣性分析命令集core_diversity_analyses.py 可全面分析不同取樣部位之間物種的兩類多樣性指數(shù)Alpha 多樣性指數(shù)（α-diversity）與Beta 多樣性指數(shù)（β-diversity）的差異，更重要的是可以發(fā)現(xiàn)人體特定部位的菌群相對(duì)含量（豐度）并發(fā)現(xiàn)其優(yōu)勢(shì)菌群。在“門(mén)”和“屬”水平對(duì)菌群分類結(jié)果顯示，具體菌群分類與豐度信息可以在多樣性分析輸出目錄（“taxa_summary_plots”）中查看，該命令集還會(huì)給出不同物種的豐度在樣本組（人體部位）之間差異的統(tǒng)計(jì)學(xué)假設(shè)檢驗(yàn)結(jié)果（“group_significance_SampleType.txt”），其中擬桿菌屬（bacteroides） 在腸道中的相對(duì)含量為60.48%，遠(yuǎn)高于其它身體部位（P＜0.05）而鏈球菌屬（streptococcus）只存在與手掌部和舌部，在腸道未檢出（P＜0.05），見(jiàn)圖2。

2.3 QIIME 系統(tǒng)生成的原始系統(tǒng)發(fā)生樹(shù) FigTree 軟件打開(kāi)原始系統(tǒng)生成樹(shù)文件，F(xiàn)igTree 的運(yùn)行需要JAVA 語(yǔ)言支持，需根據(jù)軟件提示安裝相應(yīng)的JAVA運(yùn)行環(huán)境，依次點(diǎn)擊file--open--rep_set3.tre 輸入原始樹(shù)文件rep_set3.tre，運(yùn)行后生成rectangular、polar 和radial 3 種類型的系統(tǒng)發(fā)生樹(shù)圖，結(jié)果顯示其只有物種間的進(jìn)化距離關(guān)系，但“樹(shù)葉”只能用OTU的編碼表示，rectangular 形狀的樹(shù)形圖見(jiàn)圖3。

2.4 通過(guò)QIIME 和R 語(yǔ)言“ggtree”包優(yōu)化系統(tǒng)發(fā)生樹(shù)

2.4.1 優(yōu)化樹(shù)形 為了使樹(shù)文件中包含更少的樹(shù)枝以增加圖形的可分辨性，將篩選OTUs 的豐度閾值定為5%。在R 語(yǔ)言運(yùn)行環(huán)境下調(diào)用ggtree 包讀取篩選后的Newick 格式的樹(shù)文件rep_set4.tre 繪圖，定義樹(shù)形圖的顏色、線條的形狀和樹(shù)形圖布局顯示物種間進(jìn)化距離及比例尺、標(biāo)注內(nèi)部節(jié)點(diǎn)和樹(shù)枝末端、顯示OTUs 編號(hào)。根據(jù)以上條件繪制的樹(shù)形圖見(jiàn)圖4。

2.4.2 向樹(shù)圖中添加物種分類信息 通過(guò)QIIME 的assign_ assign_taxonomy.py 腳本將篩選后的OTU 代表序列進(jìn)行物種注釋，即每個(gè)OTU 對(duì)應(yīng)一個(gè)物種名稱，將物種注釋文件轉(zhuǎn)換為csv 格式并命名為sample_rep_set4_tax_assignments.csv，需要給文件中的每一列命名，其中第一列“taxa”即為系統(tǒng)發(fā)生樹(shù)文件中的OTU 編號(hào)，“taxonomy”列為精確到“種”水平的物種分類信息，最后通過(guò)“%＜+%”操作符將“taxonomy”列中與OTU 對(duì)應(yīng)的物種分類信息添加到系統(tǒng)發(fā)生樹(shù)中。運(yùn)行R 代碼生成的系統(tǒng)發(fā)生樹(shù)圖見(jiàn)圖5。

圖2 示例樣品菌群分類與相對(duì)豐度累積條圖

圖3 FigTree 軟件查看rectangular 布局的原始系統(tǒng)發(fā)生樹(shù)圖

圖4 “ggtree”包生成的沒(méi)有物種注釋的系統(tǒng)發(fā)生樹(shù)

圖5 ggtree 包繪制的包含物種分類信息注釋的系統(tǒng)發(fā)生矩形樹(shù)圖

2.4.3 系統(tǒng)發(fā)生和物種豐度聯(lián)合圖 通過(guò)QIIME 的convert_biom.py 腳本將篩選后的biom 格式的OTUs表文件轉(zhuǎn)換為純文本格式并命名為table7.from_biom_w_consensuslineage.txt，文件中第一列taxa 即為≥5%豐度的OTU 編號(hào)，其它列為每個(gè)OTU在不同標(biāo)本中的豐度（相對(duì)含量）值。運(yùn)行R 語(yǔ)言代碼生成的系統(tǒng)發(fā)生樹(shù)與OTU 豐度熱圖見(jiàn)圖6，可以比較方便的查看不同OTU 間進(jìn)化親緣關(guān)系的遠(yuǎn)近以及OTU 數(shù)量在每個(gè)標(biāo)本中的分布情況。同時(shí)，還可以按研究設(shè)計(jì)的標(biāo)本分組情況將標(biāo)本進(jìn)行合并后作圖，對(duì)于本示例數(shù)據(jù)，可以按取樣部位將標(biāo)本合并為腸道（gut）、手掌（palm）和舌部（tongue）來(lái)源共3 個(gè)組，或者將所有標(biāo)本合并為1 個(gè)組后，重新繪制系統(tǒng)發(fā)生和OTU 在不同組間豐度分布的組合圖，運(yùn)行代碼生成的系統(tǒng)發(fā)生樹(shù)與OTU 豐度分布熱圖，結(jié)果見(jiàn)圖7，其中圖7B 所示的聯(lián)合圖系將所有標(biāo)本合并后繪制系統(tǒng)發(fā)生樹(shù)、點(diǎn)圖和條形圖的聯(lián)合圖，點(diǎn)的位置和條圖的長(zhǎng)短均表示相應(yīng)的OTU 的豐度值。

圖6 系統(tǒng)發(fā)生樹(shù)和OTU 在不同標(biāo)本中的豐度分布聯(lián)合圖

圖7 系統(tǒng)發(fā)生樹(shù)和OTU 在不同標(biāo)本組中的豐度分布聯(lián)合圖

3 討論

細(xì)菌16S rRNA 基因是編碼原核生物核糖體小亞基的基因，長(zhǎng)度約為1542 bp，在細(xì)菌進(jìn)化過(guò)程中的突變率相對(duì)較小，并且其分子大小適中，是細(xì)菌系統(tǒng)分類學(xué)研究中最常用和最有用的分子標(biāo)志[4]。通過(guò)對(duì)該基因的V3～V4 單（或雙）可變區(qū)域進(jìn)行PCR 擴(kuò)增和二代測(cè)序（NGS），根據(jù)測(cè)序數(shù)據(jù)預(yù)測(cè)標(biāo)本中微生物群落的種屬信息和種屬的相對(duì)構(gòu)成比，并進(jìn)一步研究標(biāo)本中微生物群落的種屬多樣性、微生物與環(huán)境因子的相互關(guān)系，以及微生物功能基因組與宿主和環(huán)境因子之間的相互依存關(guān)系，這些均是目前微生物學(xué)、環(huán)境科學(xué)和醫(yī)學(xué)研究中的重要課題[5-8]。

用于細(xì)菌16S rRNA 基因二代測(cè)序數(shù)據(jù)分析的常用軟件主要有 Mothur、QIIME、Uparse、RDP、VAMPS、MEGAN 等，其中Mothur 和QIIME 是最常用的分析工具。QIIME 和R 語(yǔ)言相結(jié)合可以解決幾乎全部的微生物二代測(cè)序數(shù)據(jù)的整理、統(tǒng)計(jì)分析和分析結(jié)果的可視化處理。QIIME 的主要用途是運(yùn)用經(jīng)過(guò)整理后的測(cè)序數(shù)據(jù)，通過(guò)聚類分析生成可操作的分類單元（OTUs）數(shù)據(jù)表明，OTUs 表中的數(shù)據(jù)為微生物群落的種群分類和數(shù)量組成信息。以O(shè)TUs表數(shù)據(jù)作為輸入，通過(guò)R 語(yǔ)言中的相關(guān)程序包可對(duì)OTUs 表中的數(shù)據(jù)進(jìn)行圖形化展示。R 語(yǔ)言作為開(kāi)源的，面向?qū)ο蟮慕换ナ秸Z(yǔ)言，除了能進(jìn)行常規(guī)的數(shù)據(jù)處理和統(tǒng)計(jì)學(xué)分析，近年來(lái)有眾多學(xué)者針對(duì)分子生物學(xué)實(shí)驗(yàn)數(shù)據(jù)處理的要求開(kāi)發(fā)了大量的數(shù)據(jù)處理工具，并將這些工具連同部分實(shí)驗(yàn)數(shù)據(jù)以包的形式放置到R 語(yǔ)言環(huán)境當(dāng)中，使用者通過(guò)相應(yīng)代碼（命令）調(diào)用這些包完成相應(yīng)的分析任務(wù)。在上述分析任務(wù)中，通過(guò)系統(tǒng)發(fā)生樹(shù)圖展示微生物種群之間在種系發(fā)生上的相互關(guān)系是微生物組數(shù)據(jù)分析的主要內(nèi)容。目前多數(shù)生物信息學(xué)軟件生成的微生物組系統(tǒng)樹(shù)文件只包含最基本的物種間的拓?fù)潢P(guān)系，結(jié)構(gòu)比較單一。因此本研究設(shè)計(jì)了基于QIIME 系統(tǒng)和R 語(yǔ)言“ggtree”包的構(gòu)建系統(tǒng)發(fā)生樹(shù)圖操作流程，通過(guò)該操作流程，研究者可以向系統(tǒng)發(fā)生樹(shù)圖中添加實(shí)驗(yàn)設(shè)計(jì)信息，繪制內(nèi)容更加豐富的系統(tǒng)發(fā)生樹(shù)圖，滿足論文發(fā)表和研究報(bào)告的生成等科研需求。

本工作流程整合了專業(yè)分析細(xì)菌16S rRNA 基因測(cè)序數(shù)據(jù)的QIIME 系統(tǒng)和用于數(shù)據(jù)分析和圖形化展示的R 語(yǔ)言系統(tǒng)。該方法的特點(diǎn)在于分析過(guò)程全部在Windows 系統(tǒng)下進(jìn)行，其中QIIME 的工作界面雖然是基于linux 系統(tǒng)的，但可以通過(guò)VirtualBox在windows 系統(tǒng)中虛擬出一個(gè)Linux 系統(tǒng)。并且本流程是基于代碼的操作，與其他商業(yè)化程序相比具有很高的應(yīng)用靈活性。更重要的是，工作流程是從原始的測(cè)序數(shù)據(jù)出發(fā)的，而不是從由Windows 和菜單用戶界面的軟件生成的樹(shù)文件開(kāi)始，因此可以繪制多種類型的系統(tǒng)發(fā)生樹(shù)形圖。使用本流程繪制的圖形更準(zhǔn)確和美觀，最大的優(yōu)點(diǎn)是能通過(guò)各種附加信息對(duì)樹(shù)形圖進(jìn)行多方面的注釋，還可以針對(duì)系統(tǒng)發(fā)生與微生物豐度聯(lián)合作圖。