MPB：青岛大学苏晓泉组分享基于分类学和系统发育的宏基因组比较DMS算法
为进一步提高《微生物组实验手册》稿件质量，本项目新增大众评审环节。文章在通过同行评审后，采用公众号推送方式分享全文，任何人均可在线提交修改意见。公众号格式显示略有问题，建议电脑端点击文末阅读原文下载PDF审稿。在线文档(https://kdocs.cn/l/cL8RRqHIL)大众评审页面登记姓名、单位和行号索引的修改建议。修改意见的征集截止时间为推文发布后的72小时，文章将会结合有建设性的修改意见进一步修改后获得DOI在线发表，同时根据贡献程度列为审稿人或致谢。感谢广大同行提出宝贵意见。Dynamic Meta-Storms算法：基于物种水平的生物分类学和系统发育信息对宏基因组进行全面比较张玉凤1, 2, #，荆功超2, #，陈俞竹1, 徐健2，苏晓泉1, 2, *, $1青岛大学，青岛市，山东省；2单细胞中心，中国科学院青岛生物能源与过程研究所，青岛市，山东省； $现工作单位：计算机科学技术学院，青岛大学，青岛市，山东省*通讯作者邮箱: suxq@qdu.edu.cn#共同第一作者/同等贡献摘要：精确、全面地计算宏基因组样本之间的距离对于理解微生物组的β多样性起到至关重要的作用。目前针对宏基因组的距离计算方法往往忽略了物种之间的进化关系，抑或是舍弃掉了无法定位到系统发育树叶子结点的未知物种，从而导致对β多样性进行错误地解读。为解决以上问题，我们提出了Dynamic Meta-Storms（以下简称DMS）算法，可基于生物分类和系统发育信息，在物种水平上对宏基因组样本进行全面地比较。在计算两个微生物组的距离时，对于其中能够精确注释到物种（species）水平的成分，DMS算法将这些物种映射到系统发育树的叶子结点上进行比较；而对于只能注释到属或更高水平分类信息的成分，DMS算法将其动态、合理地放置到系统发育树的虚拟中间节点上，从而最大限度地利用宏基因组的信息来计算样本之间的距离。同时，得益于并行计算技术的优化，DMS通量大且消耗资源少，在单个计算节点上计算10万个宏基因组样本的距离矩阵，仅用6.4个小时即可完成。最新版本的DMS软件可在https://github.com/qibebt-bioinfo/dynamic-meta-storms下载。输入多个宏基因组样本的物种组成后，DMS便可计算出样本之间距离矩阵(distance matrix)，用于后续的β多样性分析。目前DMS软件已经内置了与MetaPhlAn2兼容的生物分类信息和系统发育树，同时也支持用户自定义的生物分类信息和系统发育树。关键词: 宏基因组分析工具，Dynamic Meta-Storms (DMS)，生物信息算法，鸟枪宏基因组，距离矩阵，β多样性仪器设备1.DMS仅需要具有约2 GB RAM（内存）的标准计算机即可支持用户的操作。为了获得最佳性能，我们建议您使用以下规格的计算机：内存：4 GB +CPU：4核+软件1.DMS软件(Jing等， 2019)可运行于Linux、Mac OS、Windows 10内置Linux子系统等类Unix操作系统中。Dynamic Meta-Storms依赖于OpenMP库来实现并行计算。常见版本的Linux系统已经安装了OpenMP，无需额外配置。在Mac OS中，需要安装支持OpenMP的编译器，建议使用homebrew软件包管理器，通过以下命令来配置OpenMP库brew install gccDMS软件已经内置了与MetaPhlAn2(Segata等， 2012; Truong等， 2015)兼容的生物分类信息和系统发育树，其中包含了MetaPhlAn2的默认数据库中所有的细菌物种。我们建议用MetaPhlAn2对宏基因组序列进行预处理，将得到的物种名称和丰度结果作为DMS的输入来计算距离矩阵（详见实验步骤2.1）。DMS内置的分类信息和系统发育树同样也适用于由mOTUs、Karken等其他软件得出来的宏基因组物种信息。同时，DMS也支持用户自定义的生物分类信息和系统发育树（详见实验步骤3.2）。实验步骤1.安装DMS我们建议选择步骤1a中自动安装的方式来配置DMS软件。但如果自动安装程序失败，可以按照步骤1b中的步骤手动安装DMS软件。a.自动安装（首选方案）1)下载安装包git clone https://github.com/qibebt-bioinfo/dynamic-meta-storms.git2)安装运行以下安装命令：cd dynamic-meta-stormssource install.sh按照上述步骤操作，该软件包可以在1分钟内安装到计算机上。示例数据集在安装包内“example”文件夹下，可以查看“example/Readme”中的内容来获取演示运行的详细信息，或直接运行：sh Readme来演示示例数据集的距离矩阵的计算。b.手动安装（备用方案）1)下载安装包git clone https://github.com/qibebt-bioinfo/dynamic-meta-storms.git2)配置环境变量将以下内容写入环境变量配置文件（一般默认的文件是“~/.bashrc”）export DynamicMetaStorms=“Specify the full path to DMS package here”export PATH=$PATH:$DynamicMetaStorms/bin/并启用环境变量source ~/.bashrc3)编译源代码cd dynamic-meta-stormsmake2.从宏基因组序列获得物种丰度信息a.用MetaPhlAn2获得宏基因组的物种组成及相对丰度信息（profiling）以单个宏基因组序列文件“sample_1.fasta”为例：metaphlan2.py sample_1.fasta --input_type fasta --tax_lev s --ingore_viruses --ignore_eukaryotes --ignore_archaea > profiled_sample_1.sp.txt得到的输出文件“profiled_sample_1.sp.txt”格式如下（表1）表1. 单个宏基因组样本的物种组成及相对丰度信息文件#SpeciesAbundances__Rothia_aeria22.78s__Actinomyces_naeslundii13.9s__Lautropia_mirabilis12.49s__Corynebacterium_matruchotii11.27s__Corynebacterium_durum10.36s__Streptococcus_sanguinis8.13s__Actinomyces_oris6.24s__Actinomyces_massiliensis5.67s__Cardiobacterium_hominis4.83s__Porphyromonas_sp_oral_taxon_2794.33其中，第一列为物种名称，第二列为相对丰度数据。如果已经由其他软件获得了该格式的物种信息或者步骤2b中的相对丰度表（如软件自带的示例数据集“example/dataset1.sp.abd”）那么可以忽略这个步骤。但我们强烈建议所有的样本都用相同的软件和参数来处理宏基因组序列。b.将多个样本的物种信息文件合并生成物种相对丰度表将多个步骤2a中生成的样本的物种信息文件路径，汇总成一个列表文件（如samples.list.txt），格式如下（表2）表2. 样本的物种信息文件路径列表文件Sample_1profiled_sample_1.sp.txtSample_2profiled_sample_2.sp.txtSample_3profiled_sample_3.sp.txtSample_4profiled_sample_4.sp.txtSample_5profiled_sample_5.sp.txt其中第一列是样本ID，第二列是物种信息文件的路径。然后运行DMS的以下命令：MS-single-to-table -l samples.list.txt -o samples.sp.abd得到的输出文件“sample.sp.abd”格式如下（表3）表3. 物种相对丰度表SampleIDSample_1Sample_2Sample_3s__Rothia_aeria22.7816.327.65s__Actinomyces_naeslundii13.91.749.32s__Lautropia_mirabilis12.4921.185.83s__Corynebacterium_matruchotii11.271.220s__Corynebacterium_durum10.367.4111.38s__Streptococcus_sanguinis8.1314.255.8s__Actinomyces_oris6.2417.1518.46s__Actinomyces_massiliensis5.6718.3217.07s__Cardiobacterium_hominis4.832.4110.95s__Porphyromonas_sp_oral_taxon_2794.33013.54其中第一行是样本ID，第一列是物种名称，表中的数值是某物种在某样本中的相对丰度。如果已经获得了以上格式的物种相对丰度表（如软件自带的示例数据集“example/dataset1.sp.abd”），那么可以忽略这个步骤。3.使用DMS计算距离矩阵a.基于DMS内置系统发育树和物种分类信息计算距离矩阵：MS-comp-taxa-dynamic -T samples.sp.abd -o samples.sp.dist其中“-T”参数指定的输入文件“samples.sp.abd”为表3格式的物种相对丰度表。由“-o”参数指定的输出文件“samples.sp.dist”即为输入的所有样本之间的DMS距离矩阵，格式如下（表4）表4. 样本之间的DMS距离矩阵Sample_1Sample_2Sample_3Sample_100.0710880.070619Sample_20.07108800.088648Sample_30.0706190.0886480其中第一行和第一列是样本ID，表中的数值是样本之间的DMS距离，在0-1之间。数值越大，表示距离越远。b.基于用户自定义的系统发育树和物种分类信息计算距离矩阵1)自定义系统发育树和物种分类信息并生成新的DMS库自定义DMS库需要newick格式的二叉系统发育树（如“tree.newick”），树的叶子结点为物种名称，并提供所有叶子结点物种从并提供从Kingdom到Species水平的完整的生物分类信息（如“tree.taxonomy”），格式如下（表5）表5. 完整的生物分类信息文件KingdomPhylumClassOrderFamilyGenusSpeciesk__Archaeap__Euryarchaeotac__Methanopyrio__Methanopyralesf__Methanopyraceaeg__Methanopyruss__Methanopyrus_kandlerik__Archaeap__Euryarchaeotac__Methanobacteriao__Methanobacterialesf__Methanothermaceaeg__Methanothermuss__Methanothermus_fervidusk__Archaeap__Euryarchaeotac__Methanobacteriao__Methanobacterialesf__Methanobacteriaceaeg__Methanothermobacters__Methanothermobacter_thermautotrophicusk__Archaeap__Euryarchaeotac__Methanobacteriao__Methanobacterialesf__Methanobacteriaceaeg__Methanothermobacters__Methanothermobacter_marburgensisk__Archaeap__Euryarchaeotac__Methanobacteriao__Methanobacterialesf__Methanobacteriaceaeg__Methanobacteriums__Methanobacterium_formicicum利用以上文件，生成自定义的DMS库：MS-make-ref -i tree.newick -r tree.taxonomy -o tree.dms输出的“tree.dms”即为用户自定义的DMS库。2)根据自定义的系统发育树和物种分类信息计算样本之间的DMS距离矩阵：MS-comp-taxa-dynamic -D tree.dms -T samples.sp.abd -o samples.sp.dist4.DMS软件包中工具的汇总介绍a.基本工具1)MS-comp-taxa-dynamic： 计算宏基因组间的DMS距离矩阵，示例用法见3.1。可以运行MS-comp-taxa-dynamic -h了解详细的参数信息。2)MS-comp-taxa：计算宏基因组间普通的meta-storms距离矩阵(Su等， 2012; Su等， 2014) ，该方法忽略了宏基因组中未注释到物种水平的成分。用法与3.1中MS-comp-taxa-dynamic类似，可以运行MS-comp-taxa –h了解详细的参数信息。b.高级工具1)MS-single-to-table：将输出的多个单样本文件（表1）合并到同一张相对丰度表中（表3），示例用法见步骤2b。可以运行MS-single-to-table –h了解详细的参数信息。2)MS-table-to-single：将相对丰度表拆分为多个单样本输出文件，是MS-single-to-table相反的操作。可以运行MS-table-to-single –h了解详细的参数信息。3)MS-make-ref：根据自定义系统发育树和生物分类信息生成自定义DMS库，示例用法见步骤3b。可以运行MS-make-ref -h了解详细的参数信息。计算方法1.采用Meta-Storms算法计算可以识别的物种间的距离在比较宏基因组样本时，对于可以直接识别的物种，DMS使用Meta-Storms算法(Su等， 2012; Su等， 2014)基于物种水平的系统发育树计算两者之间的距离。具体来说，Meta-Storms算法计算出两个样本在叶子结点（物种）上共享的丰度（公式1），并根据系统发育树将叶子节点上剩余的丰度加权后作为其公共父节点的丰度（公式2）。最终，通过遍历树的所有节点得到样本间整体相似度（图1）。
(1)
(2)如图1所示，有两个宏基因组样本S1和S2，它们的系统发育树是典型的二叉树，其中有三个物种X，Y，Z及每个物种的比例。第一步是根据S1和S2在X和Y上共享的丰度获得相似性。然后，将X和Y的剩余丰度乘以1-Dist作为它们的共同祖先的丰度（图1B），并继续在X和Y的祖先结点与叶子结点Z上进行比较（图1C）。最后，通过在根结点的比较，得到这两个样本的总体相似度为92.8％（图1D）。
图1. Meta-Storms基于系统发育树来计算微生物组之间的相似度和距离2.采用虚拟结点映射计算未识别的群落成分对于不能映射到物种水平（species level）但有更高层级（比如属、科、目等）分类信息的群落成分，它们不能直接映射到任何特定的叶结点或内部父节点，例如，在图2右侧的分类信息（taxonomy）中，属A的成分（g__A，即物种水平的s__A_unclassified）包含的四个物种在图2左侧系统发育树（phylogeny）四个不同的分支上（s__A_sp1，s__A_sp2，s__A_sp3和s__A_sp4），但该系统发育树中并没有明确的中间结点来代表属A（g__A）。例如以上四个物种在系统发育树中的共同父节点为结点b，但该结点b下也同时包含了属B的物种s__B_sp5和s__B_sp6。为解决此问题，DMS将未能识别到物种水平群落成分，根据其更高层级分类信息，动态地插入到系统发育树中适当位置的虚拟结点来计算。该虚拟结点仅包含同一分类单元下的所有子分支。例如，在图2中，将属A映射到仅包含s__A_sp1，s__A_sp2，s__A_sp3和s__A_sp4的虚拟结点b’（不包括其他属的物种）。在计算以上四个叶子结点上的相似度时不将其添加到总相似度中，而是当计算虚拟结点时，将其平均值作为虚拟节点的相似度加到总结果中。
图2. 对于缺少物种水平分类信息的群落成分，DMS根据其更高层级分类信息，动态地插入到系统发育树中适当位置的虚拟结点来计算其对距离的贡献程度结果与分析为了验证DMS算法的准确性、可靠性与计算性能，本工作采用三个宏基因组数据集（表6）对DMS算法进行测试，并与Weighted UniFrac和Bray-Curtis距离进行比较。表6. 测试数据集数据集样本量来源Synthetic Dataset 140基于48个细菌物种的模拟合成样本Real Dataset 12,355Human Microbiome Project Phase 1 (Peterson等， 2009)Synthetic Dataset 2100,000基于3,688个细菌物种的模拟合成样本以上测试中所有的数据集均可在DMS软件下载页面的“Supplementary”部分中下载。结果1：基于模拟数据的算法验证根据图3a的样本组成结构模拟合成的样本组成结构，模拟4组宏基因组样本（每组10个，共40个；Synthetic Dataset 1）。预期的样本距离如图1b所示，即m1与m2之间的距离大于m1与m3之间的距离（Case I）且m1与m2之间的距离大于m1与m4之间的距离（Case II）。 分别利用DMS，Weighted UniFrac和Bray-Curtis算法计算其距离矩阵，并进行主坐标分析（PCoA）来获得其β多样性分布。
图3. 模拟样本的组成结构（a）与预期的样本距离（b）
图4. 基于DMS，Weighted UniFrac和Bray-Curtis距离的PCoA结果图4的结果显示，只有DMS距离的结果与预期相符，Bray-Curtis距离由于没有考虑物种间的进化距离，在Case I和II均出现错误；Weighted UniFrac距离由于忽略了未注释到系统发育树叶子结点（即species水平）的成分，在Case II出现错误。结果2：基于真实数据的算法测试将人类微生物组计划（HMP）第1阶段(Peterson等， 2009)中的2,355个真实的人体宏基因组样本作为测试数据集（来自肠道、口腔、皮肤和生殖道四个身体部位；Real Dataset 1）分别利用DMS，Weighted UniFrac和Bray-Curtis算法计算其距离矩阵，Anosim检验（1,000次重复）结果显示DMS距离能够更明显地区分样本来源的身体部位（R = 0.965, P = 0.01；图5）。
图5. 根据DMS，Weighted UniFrac和Bray-Curtis距离进行PCoA分析和Anosim检验结果3：算法运行效率测试随机选取不同数目（从10,000到100,000）的宏基因组模拟样本（Synthetic Dataset 2），利用DMS算法计算其距离矩阵，并将Striped UniFrac算法(Mcdonald等， 2018)作为参照，比较两者的总体运行时间和RAM（内存）的最大使用值。所有测试均在同一个具有80个线程的非共享计算节点上完成。图6中的结果显示，DMS计算10万个宏基因组样本的距离矩阵，仅用6.4个小时即可完成，比参照方法快20％，且节省了40％以上的内存消耗。随着宏基因组样本的数量迅速增长，此功能将发挥更大的作用。
图6. DMS与Striped UniFrac计算不同数量样本的距离矩阵所消耗的运行时间和内存使用量的对比失败经验问题1.安装提示：“make: g++: command not found”问题原因：没有安装DMS所需要的g++编译器。解决方法：根据不同的操作系统，利用相应的命令安装g++，常见的操作系统：Ubuntu系统：sudo apt-get install g++CentOS系统：sudo yum install g++MacOS系统：brew install gcc问题2.运行提示：“MS-comp-taxa-dynamic: command not found”问题原因：环境变量设置失败。解决方法：请参考实验步骤1.2.2中手动配置环境变量的方法将DMS所需要的环境变量添加到配置文件中。问题3.运行提示：“Error: Please set the environment variable "DynamicMetaStorms" to the directory”问题原因：DMS的库文件没有配置到环境变量中。解决方法：请参考实验步骤1.2.2中手动配置环境变量的方法将DMS的库文件配置到环境变量中。问题4.运行提示：“Error: Cannot open file: XXX”。问题原因：输入了错误的输入/输出文件路径。解决方案：请检查正确的输入文件路径（可在输入时用Tab键自动补全），并确保用户在输出路径下有足够的写权限。问题5.运行提示：“Argument #X Error : Arguments must start with -”。问题原因：运行命令中所有参数选项名称必须以“-”开头。解决方法：请检查第X个参数并更正。问题6.运行提示：“Error: Features: s__XXXX does not have XX Samples”问题原因：输入的物种丰度表中，物种“s__XXXX”所在行的丰度信息列数与样本数不统一。解决方法：请按照表3的格式，检查样本数量（第一行）与“s__XXXX”所在行的丰度信息个数是否一致。如果某样本中不包含该物种，则丰度标为0。如果输入的是该文件的转置格式（即每一行表示一个样本，每一列表示一个物种），请在运行MS-comp-taxa和MS-comp-taxa-dynamic时增加参数“-R F”。致谢本项工作得到了国家自然科学基金委员会31771463和32070086，中国科学院KFZD-SW-219-5，山东省自然科学基金会ZR2017ZB0421和ZR201807060158，中国博士后科学基金会2018M630807的资助。参考文献1. Jing, G., Zhang, Y., Yang, M., Liu, L., Xu, J. and Su, X. (2019). Dynamic Meta-Storms enables comprehensive taxonomic and phylogenetic comparison of shotgun metagenomes at the species level. Bioinformatics(7): 7. https://doi.org/10.1093/bioinformatics/btz9102. Mcdonald, D., Vázquez-Baeza, Y., Koslicki, D., Mcclelland, J., Reeve, N., Xu, Z., Gonzalez, A. and Knight, R. (2018). Striped UniFrac: enabling microbiome analysis at unprecedented scale. Nature Methods 15. https://doi.org/10.1038/s41592-018-0187-83. Peterson, J., Garges, S., Giovanni, M., McInnes, P., Wang, L., Schloss, J. A., Bonazzi, V., McEwen, J. E., Wetterstrand, K. A., Deal, C., Baker, C. C., Di Francesco, V., Howcroft, T. K., Karp, R. W., Lunsford, R. D., Wellington, C. R., Belachew, T., Wright, M., Giblin, C., David, H., Mills, M., Salomon, R., Mullins, C., Akolkar, B., Begg, L., Davis, C., Grandison, L., Humble, M., Khalsa, J., Little, A. R., Peavy, H., Pontzer, C., Portnoy, M., Sayre, M. H., Starke-Reed, P., Zakhari, S., Read, J., Watson, B., Guyer, M. and Grp, N. H. W. (2009). The NIH Human Microbiome Project. Genome Research 19(12): 2317-2323. https://doi.org/10.1101/gr.096651.1094. Segata, N., Waldron, L., Ballarini, A., Narasimhan, V., Jousson, O. and Huttenhower, C. (2012). Metagenomic microbial community profiling using unique clade-specific marker genes. Nature Methods 9(8): 811-+. https://doi.org/10.1038/Nmeth.20665. Su, X., Wang, X., Jing, G. and Ning, K. (2014). GPU-Meta-Storms: computing the structure similarities among massive amount of microbial community samples using GPU. Bioinformatics(7): 1031-1033. https://doi.org/10.1093/bioinformatics/btt7366. Su, X., Xu, J. and Ning, K. (2012). Meta-Storms: efficient search for similar microbial communities based on a novel indexing scheme and similarity score for metagenomic data. Bioinformatics 28(19): 2493. https://doi.org/10.1093/bioinformatics/bts4707. Su, X. Q., Wang, X. T., Jing, G. C. and Ning, K. (2014). GPU-Meta-Storms: computing the structure similarities among massive amount of microbial community samples using GPU. Bioinformatics 30(7): 1031-10338. Truong, D. T., Franzosa, E. A., Tickle, T. L., Scholz, M., Weingart, G., Pasolli, E., Tett, A., Huttenhower, C. and Segata, N. (2015). MetaPhlAn2 for enhanced metagenomic taxonomic profiling. Nature Methods 12(10): 902-903. https://doi.org/10.1038/nmeth.3589