宏基因组分析流程 | 全新的 bioBakery 3

本文整理自 Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3：https://www.biorxiv.org/content/10.1101/2020.11.19.388223v1。本人水平有限，很多时候不能准确表达文章的意思，同时对侧重点做了删减，建议感兴趣的同学阅读原文。

为了进一步扩大微生物群落研究的范围，新版本的 bioBakery 平台中引入了一套全新的计算方法。BioBakery 3 流程包括新版本的数据质控软件 KneadData，分类学注释分析软件 MetaPhlAn 3，功能分析软件 HUMAnN 3，菌株分析软件 StrainPhlAn3 和 PanPhlAn 3 以及系统发育分析软件 PhyloPhlAn 3。同时这些工具几乎都使用了新版本的 ChocoPhlAn 3 泛基因组数据库，该数据库包含系统组织和注释的微生物基因组和基因家族簇。基准测试显示，每个单独的工具都比其先前版本和其他类似算法更准确且有效，有的甚至在灵敏度和特异性上提高了两倍以上。所有组件都已开源，具体请访问 http://segatalab.cibio.unitn.it/tools/biobakery 和 http://huttenhower.sph.harvard.edu/biobakery

bioBakery 提供了完整的宏组学工具套件和分析环境，包括用于宏组学的上游分析方法，下游统计分析，可复现的集成工作流程，储存于开源库（GitHub，Conda，PyPI， 和 R/Bioconductor）的标准化软件包和文档，可网格和云部署的镜像（AWS，GCP 和 Docker），在线学习教程和演示数据以及一个公共讨论社区。对于任何样本集，都可以通过该流程生成质量控制，分类分析，功能分析，菌株分析以及所得的输出数据和报告文本，同时保持版本控制和日志记录。新版本中已更新了所有方法，并进行了相应优化。例如，优化 Docker 镜像大小以方便云部署，并将流程分别移植到 AWS（Amazon Web Services）和 Terra/Cromwell（Google Compute Engine）。所有基本镜像和依赖也已更新，包括最新的 Python（v3.7+）和 R（v4.0+）。

高质量的参考数据集可改善宏基因组分析

bioBakery 3 套件中的大多数软件都利用了全新的参考数据集 ChocoPhlAn 3。ChocoPhlAn 使用公开可用的基因组和标准化的基因和基因家族信息来生成 marker。这些标记将用于 MetaPhlAn 3，StrainPhlAn 3 和 PanPhlAn 3 进行宏基因组的分类和菌株水平分析，PhyloPhlAn 3 进行基因组和 MAGs 的系统发育分析以及 HUMAn 3 进行宏基因组的功能分析。

ChocoPhlAn 3 是一个包含了来自 16.8k 个物种的 99.2k 个高质量，全注释的参考微生物基因组以及对应的 87.3M UniRef90 基因家族注释数据库。通过这些数据，ChocoPhlAn 会为每个物种生成带功能注释的泛基因组。这些泛基因组为随后的整个 bioBakery 3 流程提供了统一的参考注释资源，HUMAnN 3 和 PanPhlAn 3 直接基于完整的泛基因组进行整体功能和菌株谱分析，而其他工具则只使用一些附加信息进行注释，比如，PhyloPhlAn 3 着眼于保守的核心基因家族（即该基因家族存在于一个物种的所有菌株中）以推断准确的系统发育，MetaPhlAn 3 和 StrainPhlAn 3 则将核心基因家族细化为物种特有的独特基因家族，以产生用于宏基因组物种鉴定和菌株水平遗传表征的生物标记。

MetaPhlAn 3 提高了定量分类学分析的准确性

MetaPhlAn 使用进化枝特异性标记基因的覆盖度估算宏基因组中微生物类群的相对丰度。MetaPhlAn 3 整合了 13.5k 个 species（比 MetaPhlAn 2 多出两倍多），并整合了全新的 ChocoPhlAn 3 数据库，其包含从 16.8k 物种泛基因组中选择的一组全新的 1.1M 标记基因（每个 species 平均 84±47）。采用 UniRef90 基因家族可进一步扩展核心基因的鉴定程序，然后将潜在核心基因与所有可用的整个微生物基因组进行比对，以确保这些标记的独特性。此外，核心算法也进行了多项改进和升级，比如标记物比对的质量控制过程和可探究由未知微生物组成的宏基因组含量。

通过 OPAL 基准框架，使用来自 2nd CAMI Challenge 的 113 个合成样本的 118 个合成宏基因组评估了 MetaPhlAn 3 的分类学性能。这些代表了来自五个人体相关身体部位和鼠类肠道的典型微生物群，并为它们补充了另外 五 个新生成的合成的非人类相关的宏基因组。除了 MetaPhlAn 3 外，还比较了 MetaPhlAn 2.7，mOTU 2.51 和 Bracken 2.5。多个评价指标均显示 MetaPhlAn 3 胜过其他算法。

对于 F1 score，MetaPhlAn 3 的表现无论在哪种样本中都优于其他方法。相对于其他工具，MetaPhlAn 3 检出的假阳性菌种数量极少。对于 recall，Bracken 和 mOTU 在某些情况下要优于 MetaPhlAn 3，但有很高的假阳性。MetaPhlAn 通过设置更高比例的阳性物种（ --stat_q 参数）可降低假阳率，但总体而言，默认设置下的 F1 值仍然高于其他算法。

除了更准确的物种检测外，在量化物种丰度上 MetaPhlAn 3 也比 MetaPhlAn 2，mOTUs2 和 Bracken 更准确。

在计算效率方面，MetaPhlAn 3 比 MetaPhlAn 2 快 3 倍以上，并且几乎与 Bracken 的速度（每秒11k reads）差不多。MetaPhlAn 3 的内存使用量比 MetaPhlAn 2（2.1Gb）略高（2.6Gb），但优于其他方法（mOTU 为 4.3 Gb，Bracken 为 32.5 Gb）。

HUMAnN 3 准确量化物种对群落功能的贡献

HUMAnN 3 使用来自 ChocoPhlAn 泛基因组的 UniRef90 功能注释数据以进行物种功能分析。同样用上文介绍的 30 个 CAMI 和 5 个其他合成宏基因组比较了 HUMAnN 2 和最近发表的 Carnelian 算法的性能。之所以选择 Carnelian 是因为它是在 HUMAnN 2 之后发表的，同时其算法策略也类似 HUMAnN ，即直接从宏基因组测序数据而不是从组装的重叠群估计分子功能的相对丰度。虽然 HUMAnN 2 和 3 都可以通过宏基因组估计多种功能注释的相对丰度，但选择了4级酶促（EC）注释作为与 Carnelian 比较的基础，因为 Carnelian 为 EC 量化提供了基准指标。

通过 HUMAnN 3 流程，产生了 30 个 CAMI 宏基因组中社区水平 EC 丰度的高准确率估算值（Bray-Curtis 相似性平均值 = 0.93±0.03）。HUMAnN 2 的注释准确率紧随其后为 0.70±0.04，而 Carnelian 的准确率仅 0.49±0.04。虽然 HUMAnN 3 受益于更新的序列数据库，但 HUMAnN 2 的数据库（2014年）也比 Carnelian 早了数年。

与其他功能分析流程相比（比如 Carnelian），HUMAnN 3（和2）的主要优势之一是它们能够根据贡献物种对群落功能概况进行分层。此功能在 HUMAnN 3 中更为准确和有用，因为他有更广泛的泛基因组注释。在 CAMI 数据中，对至少具有 1x 平均深度的物种，在 HUMAnN 3 中 EC 的准确率为 0.81±0.16，而 HUMAnN 2 为 0.51 ± 0.15 。与 HUMAnN 2 相比，HUMAnN 3 倾向于在此覆盖度内检测到更多预期物种。HUMAnN 的物种间功能灵敏度对于样本中 1x 深度以下的物种自然较低，但对每个物种的准确率仍然很高，并且新版本中，v3 与 v2 相比略有提高 (0.95 ± 0.08 vs. 0.91 ± 0.07)。

Carnelian 是这三种方法中计算效率最高的方法，分析 CAMI 的时间为 26.4±2.7 CPU 小时，而 HUMAnN 2 为 38.1±12.8 CPU小时，HUMAnN 3 为52.5±19.2 CPU小时。峰值内存使用量（MaxRSS）的趋势与此类似，其中 Carnelian 要求 11.9±0.0 GB，而 HUMAnN 2 为 17.0±0.3 GB，HUMAnN 3 为 21.5±1.9 GB。这些差异主要归因于算法检索的数据库大小：尽管 Carnelian 仅关注 EC 可注释的序列子集，但 HUMAnN 的目标是首先量化数千万个唯一的 UniRef90，其中最终只有 12.5％ 是由 EC 注释。与 HUMAnN 2 相比，HUMAnN 3 运行时间的增加同样归因于其更大的翻译检索数据库（UniRef90 87.3M v.s. 23.9M），因为即使在前面的核苷酸水平检索层中已解释了大多数数据，翻译后的检索层也是 HUMAnN 算法的速率限制步骤。这种现象也解释了 HUMAnN 的运行时间差异，因为运行时间与翻译后的检索层数据成反比。值得注意的是，当翻译检索步骤，HUMAnN 3 可以在 5.8±0.8 CPU小时中解释大多数 CAMI 宏基因组读数（每个样本70.9±9.6％），当然这通常仅适用于已知被相关参考序列充分覆盖的菌群。

使用一组富含非人类相关物种的合成数据进行评估，可以看到这三种方法之间的相对准确率和运行效率趋势相似。因此，HUMAnN 3 的强大性能不仅限于由宿主相关物种组装而成的微生物群落。此外，相对于 HUMAnN 2，MetaPhlAn 3 对非宿主相关物种的灵敏度提高了，从而提高了 HUMAnN 3 的准确性和性能（通过更快更准确的泛基因组检索步骤注释了绝大多数 reads）。

算法原理

bioBakery workflow

GitHub : https://github.com/biobakery/biobakery

bioBakery3 是一套从宏组学数据中分析微生物群落的计算方法，可产生分类，功能，系统发育和菌株水平概况，以直接解释或包括在下游统计分析中。

用 KneadData 进行 reads 质量控制后，MetaPhlAn3 将估计样品中存在的微生物种类及其相对丰度。进一步，StrainPhlAn3 通过完善由 MetaPhlAn3 鉴定的菌种的菌株水平基因型来加深遗传学表征。而 HUMAnN3 则更侧重于鉴定和定量在宏基因组中编码或在宏转录组中表达的分子功能，PanPhlAn3 可以解析其菌株的基因型。PhyloPhlAn3 提供了一种全面的方法来解释基于组装算法的基因组草图。这些 bioBakery3 模块都基于 ChocoPhlAn3 数据集，目前共包括 99,227 个基因组和 87.3M 个基因家族，这几乎是 bioBakery 第一个发行版所包含数据的 100 倍。

用 AnADAMA 2 串联流程

GitHub : https://github.com/biobakery/anadama2

bioBakery3 分析流程使用 AnADAMA v2 进行串联。简单来说，AnADAMA 打包了 Python 的依赖管理器 doit （ http://pydoit.org ），该模块为分析任务定义，版本和输出跟踪，变更管理，文档，网格，云部署和自动报告提供了一种简单但可扩展的方式。如果修改了工作流程或更改了输入文件，则只会重新运行那些受影响的任务。使用默认的记录器和命令行报告器记录所有任务的基本信息，确保流程的可重复性。记录的信息包括流程的各项命令和参数，可执行文件的版本以及每个任务的输出文件。同时也可用于将其他 bioBakery3 模块连接在一起。

用 KneadData 进行 reads 质控

GitHub : https://github.com/biobakery/kneaddata

bioBakery 3 包括一个原始数据质控模块 KneadData，该模块可自动化清洗原始宏基因组和宏转录组数据。主要内容包括：

• 使用 Trimmomatic 去除 • 1）低质量的碱基（默认值：4-mer 窗口，平均 Phred 质量 <20） • 2）截短的 reads（默认值：<50％ 的预修剪长度） • 3）adapter 和 barcode • 使用 FastQC 去除过表达序列（默认为 > 0.1％ 频率），使用 TRF 去除低复杂度的序列。但其实软件默认情况下并没有加上这个参数 --run-trim-repetitive ，因为扩增子序列通常具有大量重复 reads，会干扰结果。但强烈建议在宏基因组流程中加上该参数。 • 使用 bowtie2 比对并去除参考基因组序列。软件自带一个扩展的人类参考基因组（包含 hg37，已知的污染序列 和 decoy 序列）。

有研究（Human contamination in bacterial genomes has created thousands of spurious proteins）对 NCBI RefSeq 数据库中的细菌和古细菌基因组进行了大规模扫描，结果发现有 2250 个基因组被人类序列污染。污染物序列主要来自高拷贝的重复区域。在某些情况下，污染的重叠群被错误地注释为包含蛋白质编码序列，该序列随着时间的流逝已传播到多个原核和真核基因组中，从而形成了伪蛋白“家族”。

• 在转录组数据中使用 SILVA 数据库去除微生物核糖体和结构 RNA 序列。

建议在进行后续分析前使用 KneadData 进行质控，bioBakery workflows 工作流程会默认对所有类型的数据执行此操作。

构建 ChocoPhlAn 3 泛基因组数据库

构建 ChocoPhlAn 数据库的目的就是根据物种分类水平来组织微生物参考基因组，并整合相关序列和注释数据。在检索 UniProt 基因组和基因注释后，使用所有经过初始质控的微生物参考基因组来生成特定物种的泛基因组（比如，物种的基因家族需存在于其至少一个基因组中）。然后从整个泛基因组库中鉴定出核心基因组（即，一个物种的所有基因组中都存在的基因家族），并将其用作 PhyloPhlAn 3 中的标记。进一步处理核心基因组，以提取独特的标记基因（即，与一个物种唯一相关的核心基因家族 ），构成 MetaPhlAn 3 和 StrainPhlAn 3 的标记数据库。最后，对带有功能注释的泛基因组进行整理，构成 PanPhlAn 3 和 HUMAnN 3 的参考数据库。

泛基因组：某一物种全部基因的总称，包括核心基因组（core genome），在所有菌株中都存在的基因；以及非必须基因组（dispensable genome）， 在部分菌株中存在的基因。在实际研究中，有时候，泛基因组也可以分成核心基因组（在所有菌株中都存在的基因）、非必须基因组（在2个以及2个以上的菌株中存在的基因），以及菌株特有基因（strains-specific gene，仅在某一个菌株中存在的基因）。

数据检索

ChocoPhlAn 依赖于 UniProt 核心数据资源（release January 2019）和 NCBI 分类学和基因组存储库（release January 2019）。ChocoPhlAn 共纳入两种基本序列数据类型，即微生物的原始基因组以及这些基因组上鉴定的微生物蛋白质/基因数据。ChocoPhlAn 用微生物分类整合基因组数据，用蛋白家族整合微生物蛋白质数据。

ChocoPhlAn 采用 NCBI 分类数据库进行分类注释。完整数据可从 NCBI FTP 服务器 ftp.ncbi.nlm.nih.gov/pub/taxonomy/ 下载。ChocoPhlAn 使用 unidentified ， sp. ， Candidatus ， bacterium 来标记 species。也有一些关键词作为低质量的菌种的标志，具体来说，用于过滤低质量分类注释的正则表达式如下所示：

“ (C|c)andidat(e|us) | _sp(_.*|$) | (.*_|^)(b|B)acterium(_.*|) | .*(eury|)archaeo(n_|te|n$).*   |.*(endo|)symbiont.* | .*genomosp_.* | .*unidentified.* | .*_bacteria_.* | .*_taxon_.* | .*_et_al_.* |.*_and_.* 




    
| .*(cyano|proteo|actino)bacterium_.*) ”

根据 UniProt 蛋白质组资源中的 GCA 编号从 ftp.ncbi.nlm.nih.gov/genomes/genbank/assembly_summary_genbank.txt 文件中获取对应 URL，再从 NCBI GenBank FTP 服务器 ftp.ncbi.nlm.nih.gov/genomes/all/GCA 下载每个蛋白质组的参考基因组（ fasta 格式，后缀 .fna ）和相关的基因组注释（ .gff ）。丢弃了 12,598 个缺少 GenBank 登录号的蛋白质组，最终得到 99,227 个基因组（997 古细菌，97,941 细菌，339 真核生物）。

从 UniProt Knowledgebase（UniProtKB）和 UniProt Archive（UniParc）数据库中检索与至少一个 UniProt 蛋白质组相关并包含在 ChocoPhlAn 数据库中的微生物蛋白质（和基因）。UniProtKB 中包含的蛋白质已从 UniProt 蛋白质组中所有可用参考基因组的 CDS 的翻译衍生而来。ChocoPhlAn 3 还检索并包含了 UniProtKB 条目中存在的相关数据（从 ftp.uniprot.org/pub/databases/uniprot/ 下载 uniprot_sprot.xml.gz ， uniprot_trembl.xml.gz 和 uniparc_all.xml.gz ），比如，功能，系统发育和蛋白质结构域注释（KEGG，KO，EggNOG，GO，EC，Pfam），与外部数据库（GenBank，ENA和BioCyc）进行交叉引用条目的登录名和该基因的名称。

共处理了 UniProtKB 和 UniParc 中 203.9M 蛋白质，其中 126.9M 的蛋白质与 UniProt 蛋白质组条目相关。Bacteria domain 的蛋白质数量最多（194.8M），而古细菌和真核生物分别占 5.0M 和 4.0M 蛋白质。

为了减少数据库冗余，使用 UniProt 提供的 UniProtKB UniRef90。UniRef90 是通过将 UniRef100（将相同的 UniProtKB 蛋白结合在一个簇中）用 CD-HIT 聚类得到的（2019 年 8月后换用 MMseqs2 进行聚类）。UniRef90 中至少有 90％ 的序列同一性。UniRef50 是通过将 UniRef90 种子序列进行聚类而生成的，每个簇都包含具有至少 50％ 同一性的蛋白质。UniRef90 和 UniRef50 都要求每种蛋白质与簇的最长序列重叠至少80％。在 ChocoPhlAn 3 中考虑的 UniRef 条目包含一个代表性蛋白质的序列，该簇中包括的所有条目的登录 ID，UniProtKB 和 UniParc 记录的登录以及其他关联的 UniRef 簇的登录号都包含在 UniProt 条目中。

共处理了 292.1M 个 UniRef cluster（UniRef100，UniRef90和UniRef50分别为172.3M，87.3M和32.5M），并与 ChocoPhlAn 3 中的每种蛋白质和每个基因组相关联。

生成泛蛋白质组

接着为至少由一个 UniProt 蛋白质组表示的物种生成他们的泛蛋白质组，即该物种具有蛋白编码潜力的非冗余序列。通过整合基因组中独特的 UniRef90 和 UniRef50 蛋白家族，可以得到每个物种的蛋白数据。

对于每种泛蛋白，将计算几个分数：

• 核心评分 coreness ：该物种的泛蛋白质组中包含的基因组数目 • 唯一性评分 uniqueness ：其他具有相同泛蛋白质的物种的泛蛋白质组数目 • uniqueness_sp 得分：唯一性评分的一种变换形式，即不包括那些以前被标记为低质量物种的物种 • external_genomes ：拥有相同泛蛋白的其他物种的“泛蛋白质组”的基因组数目（而不是物种或物种的泛蛋白质组）。

这些分数是针对 UniRef50 和 UniRef90 蛋白家族计算的。

在 ChocoPhlAn 3 中，共纳入了 22,096 种泛蛋白质组和总共 87.3M UniRef90 核心蛋白。

生成 MetaPhlAn 3 标记

MetaPhlAn 依赖于一套独特的且特定于物种的核苷酸标记。首先使用 species-level genome bin（SGB）过滤掉了先前被标记为低质量分类的物种。归入同一 SGB 的低质量物种被合并，仅考虑代表 SGB，这里共合并了 1,328 个物种（6％），因为它们可能导致假阳性的分类学注释。对于 NCBI 分类中包含在“species-group”中的多个物种显示出大量具有较高唯一性评分（>30）的标记，主要包括下面这些菌种群：Streptococcus anginosus group，Lactobacillus casei group，Bacillus subtilis group，Enterobacter cloacaecomplex，Pseudomonassyringae group，Pseudomonas stutzeri group，Pseudomonas putida group，Pseudomonasfluorescens group，Pseudomonas aeruginosa group，Streptococcus dysgalactiae group 以及 Bacillus cereus group。在这些情况下，泛基因组都是通过合并所有物种级别的泛基因组并将它们视为一个物种来构建的。

第一步，使用 UniRef90 中长度在 150 至 1500 个氨基酸之间的蛋白质构建泛蛋白组数据库。从核心评分和唯一性评分开始，应用迭代的方法来查找标记，最多保留 150 个标记，并仅保留至少有 10 种标记的物种。根据“唯一性”评分将候选标记分为唯一标记和准标记：“唯一性”为零的标记被报告为“唯一标记”。如果无法识别出唯一的标记，则将使用更宽松阈值得到所谓的“准标记”，即“唯一性”评分大于零。

迭代的方法主要是结合“核心”，“唯一性”和“ external_genomes”评分来定义四个标记质量等级。

• 等级 A：核心评分超过 80％，不超过与 2 个其他泛蛋白质组共享（ Uniqueness_NR90 和 Uniqueness_NR50 ），并且出现在不超过 10 个基因组中（UniRef90，在 UniRef50 中为 5）（即 External_genomes_NR90 和 External_genomes_NR50 ） • 等级 B：核心评分介于 70％ 和 80％ 之间， Uniqueness_NR90 和 Uniqueness_NR50 值为 5 ， External_genomes_NR90 和 External_genomes_NR50 值分别低于 15 和 10 个基因组 • 等级 C：不符合先前条件的标记将包含在 C 层中，其中包括核心评分在 50％ 和 70％ 之间， Uniqueness_NR90 小于10， Uniqueness_NR50 小于15， External_genomes_NR90 小于25， External_genomes_NR50 小于20。 • 等级 U：对于泛蛋白质组中仅包含一个基因组的物种，其核心定义微不足道的标记被分类为“ U”级，它们的唯一性评分为零。

对迭代过程使用如下定义的评分对满足阈值的候选标记进行排序：

该评分将优先选择进化枝中高度保守（核心评分高）但与其他物种的数量尽可能少（唯一性评分低）的候选标记。为每个候选标记注释其等级，如果识别出 50 个以上的候选标记，那么将从排名列表中选择至多 150 个标记。如果没有找到足够的标记（少于50个），则使用下一层的阈值重复该过程。如果使用 C 级阈值仍没有找到符合的标记，则将物种丢弃。

然后，将通过此过程选择的每个标记的核苷酸序列作为 MetaPhlAn 数据库的条目。为了完善通过“唯一性”参数估算的物种数量，将标记序列拆分为150bp 的非重叠块，并使用 bowtie2（ version 2.3.4.3，参数 -a --very-sensitive --no-unal --no-hq --no-sq ）进行比对。如果在识别出的目标参考基因组中发现了至少 150 个标记序列的连续核苷酸，将基于“唯一性”参数注释一个新识别的物种。

对由 Co-Abundance gene group（CAG）获得的基因组物种的标记物进行了额外的质控步骤。为了减少假阳性的数量，如果与提供 CAG 基因组分类的物种共享超过 50％ 的标记，就删除该 CAG 物种。

每个标记都在 MetaPhlAn 数据库中关联了一个条目，该数据库包括标记的物种，共享该标记的物种列表，序列长度和物种分类。病毒标记继承自 v20_m200 MetaPhlAn2 数据库。

共确定了 13,475 个物种的 1.1M 标记。

用 MetaPhlAn 3 进行物种注释

GitHub : https://github.com/biobakery/MetaPhlAn

质控后的宏基因组 reads 由 MetaPhlAn 3 调用 bowtie2 比对到 1.1M 标记基因数据库中。默认的 bowtie2 参数是 very-sensitive ，但可通过 MetaPhlAn 3 进行自定义。在 MetaPhlAn 3 中，输入文件可以是单个 FASTQ 文件，也可以是压缩文件亦或事先比对好的数据。程序内部，MetaPhlAn 3 会估算每个标记的覆盖度，并将进化枝的覆盖度计算为同一进化枝的各个标记覆盖度的平均值。然后对进化枝的覆盖度进行标准化，以获得每个分类单元的相对丰度。