第34卷第1期2021年1月
医学信息
Journal of Medical Information
Vol.34No.1
Jan.2021论著基于细菌16S rRNA基因扩增子测序数据的
丁可打阎星羽2，廉振颖2，成丽娟3，刁玉涛2，王家林4
(1.山东第一医科大学，山东泰安271016；
2.山东第一医科大学基础医学院/基础医学研究所，山东济南250062；
3.山东中医药大学第二附属医院，山东济南250001；
4.山东第一医科大学附属肿瘤医院，山东济南250117)
摘要：目的通过QIIME系统和R语言的野1ggtree”包绘制基于细菌16S rRNA基因扩增子测序数据的系统发生树图，并构建系统发生树图操作流程遥方法通过QIIME将16S rRNA基因测序数据转化为OTUs表，对OTUs表数据进行筛选和转换，用相关命令生成Newick格式的系统发生树文件，以此文件为输入，用R语言野ggtree冶包的相关函数和代码绘制不同类型的系统发生树图遥结果成功开发了绘制系统发生树图的工作流程，使用该流程对公开发表的细菌16S rRNA基因扩增子测序数据进行处理，生成了包含物种丰度和分类信息的OTUs表，以该表为输入数据，绘制了物种相对丰度累积条图和多种形式的系统发生树图遥结论使用本流程绘制的图形更准确和美观，最大的优点是能通过各种附加信息对树形图进行多方面的注释，还可以针对系统发生与微生物丰度联合作图遥 关键词：16S rRNA基因；QIIME系统;R语言“ggtree”包；系统发生树；操作流程
中图分类号：Q752文献标识码:A DOI：10.3969/j.issn.1006-1959.2021.01.019文章编号：1006—1959(2021)01—0067—08
A Method for the Mapping of Phylogenetic Tree Based on Sequencing Data
of16S rRNA Gene Amplicon
DING Ke',YAN Xing-yu2,LIAN Zhen-ying2,CHENG Li-juan3,DIAO Yu-tao2,WANG Jia-lin4
(l.Shandong First Medical University,Taian271016,Shandong,China;
2.School of Basic Medicine/Institute of Basic Medicine,Shandong First Medical University,Jinan250062,Shandong,China;
3.The Second Affiliated Hospital of Shandong University of Traditional Chinese Medicine,Jinan250001,Shandong,China;
4.Affiliated Cancer Hospital,Shandong First Medical University,Jinan250117,Shandong,China)
Abstract:Objective To draw a phylogenetic tree diagram based on the sequencing data of bacterial16S rRNA gene amplicons through the QIIME system and R language"ggtree"package,and construct the phylogenetic tree diagram operation process.Methods The16S rRNA gene sequencing data is converted into OTUs table through QIIME,the OTUs table data is filtered and converted,and the phylogenetic tree file in Newick format is generated with related commands.This file is used as input,and the R language"ggtree"package is used as input.Related functions and codes draw different types of phylogenetic tree diagrams.Results A workflow for drawing phylogenetic tree diagrams was successfully developed.This process was used to process the publicly published bacterial16S rRNA gene amplicon sequencing data,and an OTUs table containing species abundanc
e and classification information was generated,using this table as input data,the relative abundance accumulation bar chart of species and various forms of phylogenetic tree chart are drawn.Conclusion The graphics drawn by this process are more accurate and beautiful.The biggest advantage is that the tree diagram can be annotated in many ways throughseveralof additionalinformation.It can also be used to map phylogeny and microbial abundance. Key words：rRNA gene;QIIME system;R language"ggtree"package;Phylogenetic tree;Operation process
细菌微生物16S rRNA基因扩增子二代测序数据分析的后续任务之一是通过系统树形图(dendro-gram/phylogenetic tree)展示微生物种之间在进化或种系发生上的相互关系情况,dendrogram/phyloge-netic tree图分3种类型，分别为进化分支图(cladogram)，仅有拓扑结构，不能从数量上说明各分支间进化距离的大小；系统发生(发育)图(phylogram)，各分枝长度表示碱基替换数，因而能从数量上说明各分支间进化距离的大小；时序图(chronogram)则用各分枝长度表示进化时间，其中系统发生图
基金项目：山东第一医科大学学术提升计划(编号：2019QL007)
作者简介：丁可(1994.10-)，女，山东荷泽人，硕士研究生，主要从事生物统计学研究
通讯作者：王家林(1966.8-)，男，山东济南人，本科，研究员，主要从事流行病与卫生统计学研究(ph
ylogram)无疑是微生物分子遗传学研究最常用到的图形化方法⑴。目前系统发生图制作软件生成的系统发生树文件为Newick和Nexus格式，并且大多数系统发生树查看软件也主要兼容上述2种文件格式，但Newick和Nexus格式的树文件只包含物种间的遗传距离和进化拓扑关系，结构比较单一。而R 语言中的“ggtree”包对ape、ggplot2等绘图包的功能进行了优化“ggtree”除了能支持众多的树文件格式，其绘图方法更丰富，绘制的图形更准确和美观，最大的优点是能通过各种附加信息对树形图进行多方面的注释［2］。因此，本研究构建了微生物基因组扩增子二代测序数据的系统发生树形图的制作流程，通过QIIME命令代码对测序数据进行预处理，聚类生成OTUs表数据，对OTUs表格数据进行筛选和转换，并用QIIME分析管道(pipeline)的相关命令生
67
第34卷第1期2021年1月
医学信息
Journal of Medical Information
Vol.34No.1
Jan.2021
论著
成Newick格式的系统发生树文件，以此文件为输入，用R语言“ggtree”包的函数和代码绘制包含更多研究信息的不同类型的系统发生树图，现将结果报道如下。
1材料与方法
1.1数据来源分析数据来自文献所用的原始测序数据叫基于人体微生物在不同部位并随时间推移的变化，通过QIIME分析方法，只选取其中的部分数据，即每天分别从2个人的舌部、左手手掌、右手手掌和肠道共取得34个标本，在Illumina HiSeq 2000平台上进行微生物16s rDNA扩增子测序。此夕卜，从"data.qiime
1.3分析流程构建可视化系统发生树的工作流程，关键步骤是形成BIOM(Biological Observation Matrix)格式的OTUs表。QIIME系统中有多个Python 脚本用于生成OTUs表，按聚类算法分为使用从头聚类方法(de novo)和使用参考数据库聚类方法。在此工作流程中，本研究使用的是基于Greengene参考
数据库的封闭参考(Closed reference)聚类分析方法生成OTUs表数据进行分析的，具体步骤如下:①生成OTUs表文件(otu_table.biom)：目的是依据上一步生成的序列文件seqs.fna，通过与最新的Green-Genes数据库比对进行聚类分析，其主要输出是O-TUs表格(out_table.biom)，表格中的内容为记录每个OTU在每个样品(微生物落标本)中被观察到的次数；②OTUs的过滤：通过filter_otus_from_otu_table. py命令代码筛选OTUs，只保留相对含量＞1%。(或其它比例)的OTUs;③OTUs表文件格式转换：目的是将上一步生成的biom格式的OTUs表文件转换为纯文本(csv、tsv或txt格式)文件，利于下一步数据处理;④选取代表序列：由于每个OTU中的序列不完全相同，因此需要通过pick_rep_set.py脚本选取一条代表性序列作为该OTU的序列，用于后续分析;⑤物种分级注释：用assign_taxonomy.py脚本命令对每个OTU代表性序列进行物种分类信息的注释，可以认为每个OTU近似为一个物种，反之一个分类到“种(species)”或“属(genus)”水平的物种类别可以对应一个到多个OTU；⑥代表序列的比对:QIIME系统的align_seqs.py脚本提供3种序列比对方法，即Py-NAST、MUSCLE和INFERNAL，本研究中使用QIIME系统默认的PyNAST方法，它基于NAST算法，将输入序列与提供的参考序列数据比对，在数据库中到最高匹配的序列；且MUSCLE不需要提供参考序列，可用于真菌转录间隔区(internal transcribed spacer,ITS)测序比对分析；而INFERNAL利用RNA 结构和序列相似性进行比对，与PyNAST一样需要比对数据库;⑦筛选比对序列：由于上述align_seqs. py脚本通过将长度200~400bp的目的序列和16S rRNA基因的全序列比对，因此，生成的代表性序列包含空缺(gaps)，为了保留代表性序列中的有用信息以构建系统发育树，需要通过filter_alignment.py脚本对上
述代表性序列进行筛选，去除碱基空缺等无用信息;⑧建树：运用python make_phylogeny.py脚本构建系统发育树；⑨树文件的图形化：通过FigTree 软件和R语言"ggtree”包对Newick格式的系统发育树文件进行图形化处理。系统发育树文件生成与可视化的数据处理工作流程见图1o
2结果
本结果来自对34个取样标本16S rRNA基因测序数据的分析,34个标本分别来自人体的不同部位和不同取样时间，为方便起见，本分析流程没有考虑采样时间，仅将不同取样部位作为分组因素进行分析。首选通过QIIME-1.9.1系统的pick_closed_reference_otus.py脚本生成OTUs表和系统发生树文件，然后根据物种分类信息，将物种名称标记到“树叶”上，构建优化的物种聚类树文件，该树文件作为“ggtree”包的输入文件，再结合研究设计的分组信息绘制不同类型的系统发生树图。
2.1样本信息和OTUs聚类结果本示例数据来自人体的肠道(gut)、左侧手掌(left palm)、右侧手掌(right palm)、和舌部(tongue)，为了分析方便，省略了取样时间、抗生素使用情况等其它分组信息。通过封闭参考(closed reference)聚类法将有效测序序列按照逸97%的相似性归为一个OTU。采用RDP算法与Greengene16S rRNA数据库比对，并将各OTU 注释到所属的分类单元。使用Greengene数据库可以将大多数OTU序列注释到“属”水平，少部分OTU 可以鉴定到“种”水平。通过上述方法将示例数据中的共177882条序列聚类成4403个OTU，在“门”和“属”两个水平上对同类的OTU进行分类鉴定，分别归属
于24个门,615个属，有关样本的人体部位分组情况、OTU聚类和种属鉴定的具体情况见表1o OTUs聚类及分类鉴定结果与同时生成的系统发生树文件，可以作为后续统计分析，构建系统发生树图的基础。
68
第34卷第1期2021 年 1 月
医学信息
Journal of Medical Information Vol. 34 No.1Jan. 2021
论著
图1系统发育树文件生成与可视化的数据处理工作流程
表1样本测序序列与OTU 聚类情况统计表
标本取样部位有效reads 数OTU 数不同分类等级的物种数
标本取样部位有效 reads 数OTU 数不同分类等 级的物种数
门
属门
属
L1S105
gut 8795428756L3S294right palm 14752659111L1S140gut 8261403
4
44L3S313
right palm 1313
289
12
122
L1S208
gut 9223481755
L3S341right palm 124235114147L1S257gut 6649
394750L3S360right palm 125445216194L1S281
gut 7058
423
554L3S378
right palm 1503
180
536
L1S57gut 95183905
41
L4S112
right palm 974999719297L1S76gut 8615412644L4S137right palm 1131266612
165L1S8
gut
8164286
637
L4S63
right palm
10820
1341
21378
L2S155left palm 4796
62016222L5S104tongue 2483147748L2S175
left palm 542765316204L5S155
tongue 1945
120
633L2S204left palm 4047
66215235
L5S174tongue 2117117735
L2S222
left palm 4072937
20
327L5S203
tongue 2295154741
L2S240left palm 6836
275679
L5S222tongue 2712
127740L2S309
left palm 183241410160L5S240tongue 2103125643L2S357left palm 302045910
169
L6S20
tongue 7707
212654L2S382left palm
5022
48513117L6S68tongue 7145
225
651
L3S242
right palm 1097107
5
35
L6S93
tongue
827527110
57
69
第34卷第1期2021年1月
Vol.34No.1
Jan.2021医学信息
Journal of Medical Information
2.2各OTU的物种归类人体不同部位的微生物构成存在显著差异，通过QIIME系统的种多样性分析命令集core_diversity_analyses.py可全面分析不同取样部位之间物种的两类多样性指数Alpha多样性指数(琢-diversity)与Beta多样性指数([3-diversi-ty)的差异，更重要的是可以发现人体特定部位的菌相对含量(丰度)并发现其优势菌。在“门”和“属”水平对菌分类结果显示，具体菌分类与丰度信息可以在多样性分析输出目录("taxa_summa-ry_plots")中查看，该命令集还会给出不同物种的丰度在样本组(人体部位)之间差异的统计学假设检验结果("group_''”，其中拟杆菌属(bacteroides)在肠道中的相对含量为60.48%，远高于其它身体部位(P<0.05)而链球菌属(streptococcus)只存在与手掌部和舌部，在肠道未检出(P<0.05)，见图2。
2.3 QIIME系统生成的原始系统发生树FigTree软件打开原始系统生成树文件,FigTree的运行需要JAVA语言支持，需根据软件提示安装相应的JAVA 运行环境，依次点击file--open--rep_set
<输入原始树文件，运行后生成rectangular、po-lar和radial3种类型的系统发生树图，结果显示其只有物种间的进化距离关系，但“树叶”只能用OTU 的编码表示,rectangular形状的树形图见图3。
2.4通过QIIME和R语言“ggtree”包优化系统发生树2.4.1优化树形为了使树文件中包含更少的树枝以增加图形的可分辨性，将筛选OTUs的丰度阈值定为5%o在R语言运行环境下调用ggtree包读取筛选后的
Newick格式的树文件绘图，定义树形图的颜、线条的形状和树形图布局显示物种间进化距离及比例尺、标注内部节点和树枝末端、显示O-TUs编号。根据以上条件绘制的树形图见图4o
2.4.2向树图中添加物种分类信息通过QIIME的assign_assign_taxonomy.py脚本将筛选后的OTU代表序列进行物种注释，即每个OTU对应一个物种名称，将物种注释文件转换为csv格式并命名为sam-ple_rep_set4_tax_assignments.csv，需要给文件中的每一列命名，其中第一列"taxa"即为系统发生树文件中的OTU编号，“taxonomyy"列为精确到"种"水平的物种分类信息，最后通过"％<+%"操作符将"taxonomy"列中与OTU对应的物种分类信息添加到系统发生树中。运行R代码生成的系统发生树图见图5o
A
Actino b acteria
BacteroiceLes
Flirmicutes
Fustibateriii
■Dyatw b acteria
PrrM.pooncLer»B
Cs lilrimiina
Bracnyt-acte n um
二DemiQbDctET
二□ermacDCCLis
D4?]|D lin US
Pijirjic；
left palm right palm tongue
D
C
注:A：门分类水平下标本菌分类与相对丰度;B：门分类水平下人体不同部位菌分类与相对丰度;C：属分类水平下标本菌分类与相对丰度；D：属分类水平下人体不同部位菌分类与相对丰度
图2示例样品菌分类与相对丰度累积条图
70
第34卷第1期2021年1月
医学信息
Journal of Medical Information
Vol.34No.1
Jan.2021论著£
■I.__I忙叮H
9G!?731
IMC™
嗣IM'
——-——WW
W5
醫此ILb-XL'l
1~~—甲
轴悴b'用
图3FigTree软件查看rectangular布局的原始系统发生树图
注:A：斜形(oblique)树;B长方形(rectangular)树;C无根(unrooted)树;D圆形(circular)树
图4野ggtree冶包生成的没有物种注释的系统发生树
71
