丁 可,閻星羽,廉振穎,成麗娟,刁玉濤,王家林
(1.山東第一醫(yī)科大學(xué),山東 泰安 271016;2.山東第一醫(yī)科大學(xué)基礎(chǔ)醫(yī)學(xué)院/基礎(chǔ)醫(yī)學(xué)研究所,山東 濟南 250062;3.山東中醫(yī)藥大學(xué)第二附屬醫(yī)院,山東 濟南 250001;4.山東第一醫(yī)科大學(xué)附屬腫瘤醫(yī)院,山東 濟南 250117)
細(xì)菌微生物16S rRNA 基因擴增子二代測序數(shù)據(jù)分析的后續(xù)任務(wù)之一是通過系統(tǒng)樹形圖(dendrogram/phylogenetic tree)展示微生物種群之間在進(jìn)化或種系發(fā)生上的相互關(guān)系情況,dendrogram/phylogenetic tree 圖分3 種類型,分別為進(jìn)化分支圖(cladogram),僅有拓?fù)浣Y(jié)構(gòu),不能從數(shù)量上說明各分支間進(jìn)化距離的大?。幌到y(tǒng)發(fā)生(發(fā)育)圖(phylogram),各分枝長度表示堿基替換數(shù),因而能從數(shù)量上說明各分支間進(jìn)化距離的大?。粫r序圖(chronogram)則用各分枝長度表示進(jìn)化時間,其中系統(tǒng)發(fā)生圖(phylogram)無疑是微生物分子遺傳學(xué)研究最常用到的圖形化方法[1]。目前系統(tǒng)發(fā)生圖制作軟件生成的系統(tǒng)發(fā)生樹文件為Newick 和Nexus 格式,并且大多數(shù)系統(tǒng)發(fā)生樹查看軟件也主要兼容上述2 種文件格式,但Newick 和Nexus 格式的樹文件只包含物種間的遺傳距離和進(jìn)化拓?fù)潢P(guān)系,結(jié)構(gòu)比較單一。而R語言中的“ggtree”包對ape、ggplot2 等繪圖包的功能進(jìn)行了優(yōu)化,“ggtree”除了能支持眾多的樹文件格式,其繪圖方法更豐富,繪制的圖形更準(zhǔn)確和美觀,最大的優(yōu)點是能通過各種附加信息對樹形圖進(jìn)行多方面的注釋[2]。因此,本研究構(gòu)建了微生物基因組擴增子二代測序數(shù)據(jù)的系統(tǒng)發(fā)生樹形圖的制作流程,通過QIIME 命令代碼對測序數(shù)據(jù)進(jìn)行預(yù)處理,聚類生成OTUs 表數(shù)據(jù),對OTUs 表格數(shù)據(jù)進(jìn)行篩選和轉(zhuǎn)換,并用QIIME 分析管道(pipeline)的相關(guān)命令生成Newick 格式的系統(tǒng)發(fā)生樹文件,以此文件為輸入,用R 語言“ggtree”包的函數(shù)和代碼繪制包含更多研究信息的不同類型的系統(tǒng)發(fā)生樹圖,現(xiàn)將結(jié)果報道如下。
1.1 數(shù)據(jù)來源 分析數(shù)據(jù)來自文獻(xiàn)所用的原始測序數(shù)據(jù)[3],基于人體微生物在不同部位并隨時間推移的變化,通過QIIME 分析方法,只選取其中的部分?jǐn)?shù)據(jù),即每天分別從2 個人的舌部、左手手掌、右手手掌和腸道共取得34 個標(biāo)本,在Illumina HiSeq 2000 平臺上進(jìn)行微生物16s rDNA 擴增子測序。此外,從“https://data.qiime2.org/2017.6/tutorials/moving-pictures/emp -single -end -sequences/barcodes.fastq.gz”網(wǎng)站獲得經(jīng)過整理后的單向測序的數(shù)據(jù)。
1.2 軟件系統(tǒng) 通過QIIME 完成系統(tǒng)發(fā)育樹構(gòu)建和可視化的工作流程,并使用是R 語言中的“ggtree”包(3.6.0 版)進(jìn)行樹文件的可視化和注釋。
1.3 分析流程 構(gòu)建可視化系統(tǒng)發(fā)生樹的工作流程,關(guān)鍵步驟是形成BIOM(Biological Observation Matrix)格式的OTUs 表。QIIME 系統(tǒng)中有多個Python腳本用于生成OTUs 表,按聚類算法分為使用從頭聚類方法(de novo)和使用參考數(shù)據(jù)庫聚類方法。在此工作流程中,本研究使用的是基于Greengene 參考數(shù)據(jù)庫的封閉參考(Closed reference)聚類分析方法生成OTUs 表數(shù)據(jù)進(jìn)行分析的,具體步驟如下:①生成OTUs 表文件(otu_table.biom):目的是依據(jù)上一步生成的序列文件seqs.fna,通過與最新的Green-Genes 數(shù)據(jù)庫比對進(jìn)行聚類分析,其主要輸出是OTUs 表格(out_table.biom),表格中的內(nèi)容為記錄每個OTU 在每個樣品(微生物群落標(biāo)本)中被觀察到的次數(shù);②OTUs 的過濾:通過filter_otus_from_otu_table.py 命令代碼篩選OTUs,只保留相對含量>1‰(或其它比例)的OTUs;③OTUs 表文件格式轉(zhuǎn)換:目的是將上一步生成的biom 格式的OTUs 表文件轉(zhuǎn)換為純文本(csv、tsv 或txt 格式)文件,利于下一步數(shù)據(jù)處理;④選取代表序列:由于每個OTU 中的序列不完全相同,因此需要通過pick_rep_set.py 腳本選取一條代表性序列作為該OTU 的序列,用于后續(xù)分析;⑤物種分級注釋:用assign_taxonomy.py 腳本命令對每個OTU 代表性序列進(jìn)行物種分類信息的注釋,可以認(rèn)為每個OTU 近似為一個物種,反之一個分類到“種(species)”或“屬(genus)”水平的物種類別可以對應(yīng)一個到多個OTU;⑥代表序列的比對:QIIME 系統(tǒng)的align_seqs.py 腳本提供3 種序列比對方法,即Py-NAST、MUSCLE 和INFERNAL,本研究中使用QIIME 系統(tǒng)默認(rèn)的PyNAST 方法,它基于NAST 算法,將輸入序列與提供的參考序列數(shù)據(jù)比對,在數(shù)據(jù)庫中找到最高匹配的序列;且MUSCLE 不需要提供參考序列,可用于真菌轉(zhuǎn)錄間隔區(qū)(internal transcribed spacer,ITS)測序比對分析;而INFERNAL 利用RNA結(jié)構(gòu)和序列相似性進(jìn)行比對,與PyNAST 一樣需要比對數(shù)據(jù)庫;⑦篩選比對序列:由于上述align_seqs.py 腳本通過將長度200~400 bp 的目的序列和16S rRNA 基因的全序列比對,因此,生成的代表性序列包含空缺(gaps),為了保留代表性序列中的有用信息以構(gòu)建系統(tǒng)發(fā)育樹,需要通過filter_alignment.py 腳本對上述代表性序列進(jìn)行篩選,去除堿基空缺等無用信息;⑧建樹:運用python make_phylogeny.py 腳本構(gòu)建系統(tǒng)發(fā)育樹;⑨樹文件的圖形化:通過FigTree軟件和R 語言“ggtree”包對Newick 格式的系統(tǒng)發(fā)育樹文件進(jìn)行圖形化處理。系統(tǒng)發(fā)育樹文件生成與可視化的數(shù)據(jù)處理工作流程見圖1。
本結(jié)果來自對34 個取樣標(biāo)本16S rRNA 基因測序數(shù)據(jù)的分析,34 個標(biāo)本分別來自人體的不同部位和不同取樣時間,為方便起見,本分析流程沒有考慮采樣時間,僅將不同取樣部位作為分組因素進(jìn)行分析。首選通過QIIME-1.9.1 系統(tǒng)的pick_closed_reference_otus.py 腳本生成OTUs 表和系統(tǒng)發(fā)生樹文件,然后根據(jù)物種分類信息,將物種名稱標(biāo)記到“樹葉”上,構(gòu)建優(yōu)化的物種聚類樹文件,該樹文件作為“ggtree”包的輸入文件,再結(jié)合研究設(shè)計的分組信息繪制不同類型的系統(tǒng)發(fā)生樹圖。
2.1 樣本信息和OTUs 聚類結(jié)果 本示例數(shù)據(jù)來自人體的腸道(gut)、左側(cè)手掌(left palm)、右側(cè)手掌(right palm)、和舌部(tongue),為了分析方便,省略了取樣時間、抗生素使用情況等其它分組信息。通過封閉參考(closed reference)聚類法將有效測序序列按照≥97%的相似性歸為一個OTU。采用RDP 算法與Greengene 16S rRNA 數(shù)據(jù)庫比對,并將各OTU注釋到所屬的分類單元。使用Greengene 數(shù)據(jù)庫可以將大多數(shù)OTU 序列注釋到“屬”水平,少部分OTU可以鑒定到“種”水平。通過上述方法將示例數(shù)據(jù)中的共177882 條序列聚類成4403 個OTU,在“門”和“屬”兩個水平上對同類的OTU 進(jìn)行分類鑒定,分別歸屬于24 個門,615 個屬,有關(guān)樣本的人體部位分組情況、OTU 聚類和種屬鑒定的具體情況見表1。OTUs 聚類及分類鑒定結(jié)果與同時生成的系統(tǒng)發(fā)生樹文件,可以作為后續(xù)統(tǒng)計分析,構(gòu)建系統(tǒng)發(fā)生樹圖的基礎(chǔ)。
圖1 系統(tǒng)發(fā)育樹文件生成與可視化的數(shù)據(jù)處理工作流程
表1 樣本測序序列與OTU 聚類情況統(tǒng)計表
2.2 各OTU 的物種歸類 人體不同部位的微生物構(gòu)成存在顯著差異,通過QIIME 系統(tǒng)的種群多樣性分析命令集core_diversity_analyses.py 可全面分析不同取樣部位之間物種的兩類多樣性指數(shù)Alpha 多樣性指數(shù)(α-diversity)與Beta 多樣性指數(shù)(β-diversity)的差異,更重要的是可以發(fā)現(xiàn)人體特定部位的菌群相對含量(豐度)并發(fā)現(xiàn)其優(yōu)勢菌群。在“門”和“屬”水平對菌群分類結(jié)果顯示,具體菌群分類與豐度信息可以在多樣性分析輸出目錄(“taxa_summary_plots”)中查看,該命令集還會給出不同物種的豐度在樣本組(人體部位)之間差異的統(tǒng)計學(xué)假設(shè)檢驗結(jié)果(“group_significance_SampleType.txt”),其中擬桿菌屬(bacteroides) 在腸道中的相對含量為60.48%,遠(yuǎn)高于其它身體部位(P<0.05)而鏈球菌屬(streptococcus)只存在與手掌部和舌部,在腸道未檢出(P<0.05),見圖2。
2.3 QIIME 系統(tǒng)生成的原始系統(tǒng)發(fā)生樹 FigTree 軟件打開原始系統(tǒng)生成樹文件,F(xiàn)igTree 的運行需要JAVA 語言支持,需根據(jù)軟件提示安裝相應(yīng)的JAVA運行環(huán)境,依次點擊file--open--rep_set3.tre 輸入原始樹文件rep_set3.tre,運行后生成rectangular、polar 和radial 3 種類型的系統(tǒng)發(fā)生樹圖,結(jié)果顯示其只有物種間的進(jìn)化距離關(guān)系,但“樹葉”只能用OTU的編碼表示,rectangular 形狀的樹形圖見圖3。
2.4 通過QIIME 和R 語言“ggtree”包優(yōu)化系統(tǒng)發(fā)生樹
2.4.1 優(yōu)化樹形 為了使樹文件中包含更少的樹枝以增加圖形的可分辨性,將篩選OTUs 的豐度閾值定為5%。在R 語言運行環(huán)境下調(diào)用ggtree 包讀取篩選后的Newick 格式的樹文件rep_set4.tre 繪圖,定義樹形圖的顏色、線條的形狀和樹形圖布局顯示物種間進(jìn)化距離及比例尺、標(biāo)注內(nèi)部節(jié)點和樹枝末端、顯示OTUs 編號。根據(jù)以上條件繪制的樹形圖見圖4。
2.4.2 向樹圖中添加物種分類信息 通過QIIME 的assign_ assign_taxonomy.py 腳本將篩選后的OTU 代表序列進(jìn)行物種注釋,即每個OTU 對應(yīng)一個物種名稱,將物種注釋文件轉(zhuǎn)換為csv 格式并命名為sample_rep_set4_tax_assignments.csv,需要給文件中的每一列命名,其中第一列“taxa”即為系統(tǒng)發(fā)生樹文件中的OTU 編號,“taxonomy”列為精確到“種”水平的物種分類信息,最后通過“%<+%”操作符將“taxonomy”列中與OTU 對應(yīng)的物種分類信息添加到系統(tǒng)發(fā)生樹中。運行R 代碼生成的系統(tǒng)發(fā)生樹圖見圖5。
圖2 示例樣品菌群分類與相對豐度累積條圖
圖3 FigTree 軟件查看rectangular 布局的原始系統(tǒng)發(fā)生樹圖
圖4 “ggtree”包生成的沒有物種注釋的系統(tǒng)發(fā)生樹
圖5 ggtree 包繪制的包含物種分類信息注釋的系統(tǒng)發(fā)生矩形樹圖
2.4.3 系統(tǒng)發(fā)生和物種豐度聯(lián)合圖 通過QIIME 的convert_biom.py 腳本將篩選后的biom 格式的OTUs表文件轉(zhuǎn)換為純文本格式并命名為table7.from_biom_w_consensuslineage.txt,文件中第一列taxa 即為≥5%豐度的OTU 編號,其它列為每個OTU在不同標(biāo)本中的豐度(相對含量)值。運行R 語言代碼生成的系統(tǒng)發(fā)生樹與OTU 豐度熱圖見圖6,可以比較方便的查看不同OTU 間進(jìn)化親緣關(guān)系的遠(yuǎn)近以及OTU 數(shù)量在每個標(biāo)本中的分布情況。同時,還可以按研究設(shè)計的標(biāo)本分組情況將標(biāo)本進(jìn)行合并后作圖,對于本示例數(shù)據(jù),可以按取樣部位將標(biāo)本合并為腸道(gut)、手掌(palm)和舌部(tongue)來源共3 個組,或者將所有標(biāo)本合并為1 個組后,重新繪制系統(tǒng)發(fā)生和OTU 在不同組間豐度分布的組合圖,運行代碼生成的系統(tǒng)發(fā)生樹與OTU 豐度分布熱圖,結(jié)果見圖7,其中圖7B 所示的聯(lián)合圖系將所有標(biāo)本合并后繪制系統(tǒng)發(fā)生樹、點圖和條形圖的聯(lián)合圖,點的位置和條圖的長短均表示相應(yīng)的OTU 的豐度值。
圖6 系統(tǒng)發(fā)生樹和OTU 在不同標(biāo)本中的豐度分布聯(lián)合圖
圖7 系統(tǒng)發(fā)生樹和OTU 在不同標(biāo)本組中的豐度分布聯(lián)合圖
細(xì)菌16S rRNA 基因是編碼原核生物核糖體小亞基的基因,長度約為1542 bp,在細(xì)菌進(jìn)化過程中的突變率相對較小,并且其分子大小適中,是細(xì)菌系統(tǒng)分類學(xué)研究中最常用和最有用的分子標(biāo)志[4]。通過對該基因的V3~V4 單(或雙)可變區(qū)域進(jìn)行PCR 擴增和二代測序(NGS),根據(jù)測序數(shù)據(jù)預(yù)測標(biāo)本中微生物群落的種屬信息和種屬的相對構(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 基因二代測序數(shù)據(jù)分析的常用軟件主要有 Mothur、QIIME、Uparse、RDP、VAMPS、MEGAN 等,其中Mothur 和QIIME 是最常用的分析工具。QIIME 和R 語言相結(jié)合可以解決幾乎全部的微生物二代測序數(shù)據(jù)的整理、統(tǒng)計分析和分析結(jié)果的可視化處理。QIIME 的主要用途是運用經(jīng)過整理后的測序數(shù)據(jù),通過聚類分析生成可操作的分類單元(OTUs)數(shù)據(jù)表明,OTUs 表中的數(shù)據(jù)為微生物群落的種群分類和數(shù)量組成信息。以O(shè)TUs表數(shù)據(jù)作為輸入,通過R 語言中的相關(guān)程序包可對OTUs 表中的數(shù)據(jù)進(jìn)行圖形化展示。R 語言作為開源的,面向?qū)ο蟮慕换ナ秸Z言,除了能進(jìn)行常規(guī)的數(shù)據(jù)處理和統(tǒng)計學(xué)分析,近年來有眾多學(xué)者針對分子生物學(xué)實驗數(shù)據(jù)處理的要求開發(fā)了大量的數(shù)據(jù)處理工具,并將這些工具連同部分實驗數(shù)據(jù)以包的形式放置到R 語言環(huán)境當(dāng)中,使用者通過相應(yīng)代碼(命令)調(diào)用這些包完成相應(yīng)的分析任務(wù)。在上述分析任務(wù)中,通過系統(tǒng)發(fā)生樹圖展示微生物種群之間在種系發(fā)生上的相互關(guān)系是微生物組數(shù)據(jù)分析的主要內(nèi)容。目前多數(shù)生物信息學(xué)軟件生成的微生物組系統(tǒng)樹文件只包含最基本的物種間的拓?fù)潢P(guān)系,結(jié)構(gòu)比較單一。因此本研究設(shè)計了基于QIIME 系統(tǒng)和R 語言“ggtree”包的構(gòu)建系統(tǒng)發(fā)生樹圖操作流程,通過該操作流程,研究者可以向系統(tǒng)發(fā)生樹圖中添加實驗設(shè)計信息,繪制內(nèi)容更加豐富的系統(tǒng)發(fā)生樹圖,滿足論文發(fā)表和研究報告的生成等科研需求。
本工作流程整合了專業(yè)分析細(xì)菌16S rRNA 基因測序數(shù)據(jù)的QIIME 系統(tǒng)和用于數(shù)據(jù)分析和圖形化展示的R 語言系統(tǒng)。該方法的特點在于分析過程全部在Windows 系統(tǒng)下進(jìn)行,其中QIIME 的工作界面雖然是基于linux 系統(tǒng)的,但可以通過VirtualBox在windows 系統(tǒng)中虛擬出一個Linux 系統(tǒng)。并且本流程是基于代碼的操作,與其他商業(yè)化程序相比具有很高的應(yīng)用靈活性。更重要的是,工作流程是從原始的測序數(shù)據(jù)出發(fā)的,而不是從由Windows 和菜單用戶界面的軟件生成的樹文件開始,因此可以繪制多種類型的系統(tǒng)發(fā)生樹形圖。使用本流程繪制的圖形更準(zhǔn)確和美觀,最大的優(yōu)點是能通過各種附加信息對樹形圖進(jìn)行多方面的注釋,還可以針對系統(tǒng)發(fā)生與微生物豐度聯(lián)合作圖。