Dynamic Meta-Storms算法：基于物种水平的生物分类学和系统发育信息对宏基因组进行全面比较   

仪器设备

1. DMS仅需要具有约2 GB RAM (内存) 的标准计算机即可支持用户的操作。为了获得最佳性能，我们建议您使用以下规格的计算机：
   内存：4 GB +
   CPU：4核+
   

软件

1. DMS软件 (Jing et al., 2020) 可运行于Linux、Mac OS、Windows 10内置Linux子系统等类Unix操作系统中。Dynamic Meta-Storms依赖于OpenMP库来实现并行计算。常见版本的Linux系统已经安装了OpenMP，无需额外配置。在Mac OS中，需要安装支持OpenMP的编译器，建议使用homebrew软件包管理器，通过以下命令来配置OpenMP库
   
   brew install gcc
   DMS软件已经内置了与MetaPhlAn2 (Segata et al., 2012; Truong et al., 2015) 兼容的生物分类信息和系统发育树，其中包含了MetaPhlAn2的默认数据库中所有的细菌物种。我们建议用MetaPhlAn2对宏基因组序列进行预处理，将得到的物种名称和丰度结果作为DMS的输入来计算距离矩阵 (详见实验步骤2.1) 。DMS内置的分类信息和系统发育树同样也适用于由mOTUs、Kraken等其他软件得出来的宏基因组物种信息。同时，DMS也支持用户自定义的生物分类信息和系统发育树 (详见实验步骤3.2) 。
   

实验步骤

1. 安装DMS
   我们建议选择步骤1.1中自动安装的方式来配置DMS软件。但如果自动安装程序失败，可以按照步骤1.2中的步骤手动安装DMS软件。
   1.1

   自动安装 (首选方案)

   1)

   下载安装包
   

   git clone https://github.com/qibebt-bioinfo/dynamic-meta-storms.git

   2)

   安装 运行以下安装命令：
   运行以下安装命令：
   

   cd dynamic-meta-storms

   source install.sh

   按照上述步骤操作，该软件包可以在1分钟内安装到计算机上。
   示例数据集在安装包内“example”文件夹下，可以查看“example/Readme”中的内容来获取演示运行的详细信息，或直接运行：
   

   sh Readme

   来演示示例数据集的距离矩阵的计算。
   1.2

   手动安装 (备用方案)

   1)

   下载安装包
   

   git clone https://github.com/qibebt-bioinfo/dynamic-meta-storms.git

   2)

   配置环境变量
   将以下内容写入环境变量配置文件 (一般默认的文件是“~/.bashrc”)
   

   export DynamicMetaStorms=“Specify the full path to DMS package here”

   export PATH=$PATH:$DynamicMetaStorms/bin/

   并启用环境变量
   

   source ~/.bashrc

   3)

   编译源代码
   

   cd dynamic-meta-storms

   make
2. 从宏基因组序列获得物种丰度信息
   2.1

   用MetaPhlAn2获得宏基因组的物种组成及相对丰度信息 (profiling)
   以单个宏基因组序列文件“sample_1.fasta”为例：
   metaphlan2.py sample_1.fasta --input_type fasta --tax_lev s --ignore_viruses --ignore_eukaryotes --ignore_archaea > profiled_sample_1.sp.txt
   得到的输出文件“profiled_sample_1.sp.txt”格式如下 (表1)
   
   表1. 单个宏基因组样本的物种组成及相对丰度信息文件
   
   
   其中，第一列为物种名称，第二列为相对丰度数据。如果已经由其他软件获得了该格式的物种信息或者步骤2.2中的相对丰度表 (如软件自带的示例数据集“example/dataset1.sp.abd”) 那么可以忽略这个步骤。但我们强烈建议所有的样本都用相同的软件和参数来处理宏基因组序列。
   2.2

   将多个样本的物种信息文件合并生成物种相对丰度表 将多个步骤2.1中生成的样本的物种信息文件路径，汇总成一个列表文件 (如samples.list.txt) ，格式如下 (表2)
   
   表2. 样本的物种信息文件路径列表文件
   
   
   其中第一列是样本ID，第二列是物种信息文件的路径。然后运行DMS的以下命令：
   MS-single-to-table -l samples.list.txt -o samples.sp.abd
   得到的输出文件“samples.sp.abd”格式如下 (表3)
   

   
   表3. 物种相对丰度表
   
   
   其中第一行是样本ID，第一列是物种名称，表中的数值是某物种在某样本中的相对丰度。如果已经获得了以上格式的物种相对丰度表 (如软件自带的示例数据集“example/dataset1.sp.abd”) ，那么可以忽略这个步骤。
   
3. 使用DMS计算距离矩阵
   3.1

   基于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距离矩阵
   
   
   其中第一行和第一列是样本ID，表中的数值是样本之间的DMS距离，在0-1之间。数值越大，表示距离越远。
   3.2

   基于用户自定义的系统发育树和物种分类信息计算距离矩阵

   1)

   自定义系统发育树和物种分类信息并生成新的DMS库
   自定义DMS库需要newick格式的二叉系统发育树 (如“tree.newick”) ，树的叶子结点为物种名称，并提供所有叶子结点物种从Kingdom到Species水平的完整的生物分类信息 (如“tree.taxonomy”) ，格式如下 (表5)。
   
   表5. 完整的生物分类信息文件
   
   
   利用以上文件，生成自定义的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.dist
4. DMS软件包中工具的汇总介绍
   4.1

   基本工具

   1)

   MS-comp-taxa-dynamic： 计算宏基因组间的DMS距离矩阵，示例用法见3.1。可以运行
   

   MS-comp-taxa-dynamic -h

   了解详细的参数信息。

   2)

   MS-comp-taxa：计算宏基因组间普通的meta-storms距离矩阵 (Su et al., 2012; Su et al., 2014)，该方法忽略了宏基因组中未注释到物种水平的成分。用法与3.1中MS-comp-taxa-dynamic类似，可以运行
   

   MS-comp-taxa –h

   了解详细的参数信息。
   4.2

   高级工具

   1)

   MS-single-to-table：将输出的多个单样本文件 (表1) 合并到同一张相对丰度表中 (表3) ，示例用法见步骤2.2。可以运行
   

   MS-single-to-table –h 

   了解详细的参数信息。

   2)

   MS-table-to-single：将相对丰度表拆分为多个单样本输出文件，是MS-single-to-table相反的操作。可以运行
   

   MS-table-to-single –h 

   了解详细的参数信息。

   3)

   MS-make-ref：根据自定义系统发育树和生物分类信息生成自定义DMS库，示例用法见步骤3.2。可以运行
   

   MS-make-ref -h 

   了解详细的参数信息。
5. 计算方法
   5.1

   采用Meta-Storms算法计算可以识别的物种间的距离
   在比较宏基因组样本时，对于可以直接识别的物种，DMS使用Meta-Storms算法 (Su et al., 2012; Su et al., 2014) 基于物种水平的系统发育树计算两者之间的距离。具体来说，Meta-Storms算法计算出两个样本在叶子结点 (物种) 上共享的丰度 (公式1) ，并根据系统发育树将叶子结点上剩余的丰度加权后作为其公共父结点的丰度 (公式2) 。最终，通过遍历树的所有结点得到样本间整体相似度 (图1) 。
   
   
   
   
   如图1所示，有两个宏基因组样本S1和S2，它们的系统发育树是典型的二叉树，其中有三个物种X，Y，Z及每个物种的比例。第一步是根据S1和S2在X和Y上共享的丰度获得相似性。然后，将X和Y的剩余丰度乘以1-Dist作为它们的共同祖先的丰度 (图1. B) ，并继续在X和Y的祖先结点与叶子结点Z上进行比较 (图1. C) 。最后，通过在根结点的比较，得到这两个样本的总体相似度为92.8％ (图1. D) 。
   
   
   图1. Meta-Storms基于系统发育树来计算微生物组之间的相似度和距离
   
   
   5.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. 测试数据集


以上测试中所有的数据集均可在DMS软件下载页面的“Supplementary”部分中下载。

1. 基于模拟数据的算法验证
   根据图3. A的样本组成结构模拟合成的样本组成结构，模拟4组宏基因组样本 (每组10个，共40个；Synthetic Dataset 1) 。预期的样本距离如图3. B所示，即m1与m2之间的距离大于m1与m3之间的距离 (Case I) 且m1与m2之间的距离大于m1与m4之间的距离 (Case II) 。分别利用DMS, Weighted UniFrac和Bray-Curtis算法计算其距离矩阵，并进行主坐标分析 (PCoA) 来获得其β多样性分布。
   
   
   图3. 模拟样本的组成结构 (A) 与预期的样本距离 (B)
   
   图4的结果显示，只有DMS距离的结果与预期相符，Bray-Curtis距离由于没有考虑物种间的进化距离，在Case I和II均出现错误；Weighted UniFrac距离由于忽略了未注释到系统发育树叶子结点 (即species水平) 的成分，在Case II出现错误。
   
   
   图4. 基于DMS, Weighted UniFrac和Bray-Curtis距离的PCoA结果
   
   
2. 基于真实数据的算法测试
   将人类微生物组计划 (HMP) 第1阶段 (Peterson et al., 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 et al., 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的资助。

参考文献

1. Jing, G., Zhang, Y., Yang, M., Liu, L., Xu, J. and Su, X. (2020). Dynamic Meta-Storms enables comprehensive taxonomic and phylogenetic comparison of shotgun metagenomes at the species level. Bioinformatics 36(7): 2308-2310. 
2. McDonald, D., Vazquez-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(11): 847-848. 
3. Peterson, J., Garges, S., Giovanni, M., McInnes, P., Wang, L., Schloss, J. A., Bonazzi, V., McEwen, J. E., Wetterstrand, K. A., Deal, C. et al. (2009). The NIH Human Microbiome Project. Genome Research 19(12): 2317-2323.
4. 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-814. 
5. 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 30(7): 1031-1033. 
6. 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-2501. 
7. 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-90