GigaScience | 基因组所刘永鑫组-短读长宏基因组数据的去宿主污染软件评估

背景 ：宏基因组测序技术的快速发展为探索微生物组在宿主健康与疾病中的复杂作用，以及揭示微生物群落的未知结构和功能提供了重要机遇。然而，宏基因组数据的迅速积累也给数据分析带来了巨大挑战。宿主DNA的污染可能显著影响结果的准确性，并因包含非目标序列而增加额外的计算资源消耗。

结果 ：在本研究中，我们评估了计算去除宿主DNA污染对下游分析的影响，强调了其在高效生成准确结果中的重要性。此外，我们还对常用工具（如 KneadData、Bowtie2、BWA、KMCP、Kraken2 和 KrakenUniq）的性能进行了评估，这些工具在不同应用场景下各具优势。此外，我们强调了准确的宿主参考基因组的重要性，并指出其缺失会对所有工具的去污染效果产生负面影响。

结论 ：我们的研究结果强调了在宏基因组分析中慎重选择去污染工具和参考基因组的必要性。这些见解为提高微生物组研究的可靠性和可重复性提供了有价值的指导。

第二代测序技术和数据分析方法的进步极大地促进了微生物组研究，拓宽了我们对微生物组对宿主广泛影响的认知。相比于扩增子测序，宏基因组鸟枪法测序能够以更少的偏倚对细菌群落进行全面评估，并在物种、菌株或功能水平上提供更高的解析度。随着测序成本的迅速下降和测序深度的不断增加，宏基因组数据的规模呈指数级增长。一些研究甚至在单个样本中生成超过 100 Gbps（十亿碱基对）的数据，以解析人类肠道微生物组中的“暗物质”。

然而，在分析复杂的宿主相关微生物群（如唾液、咽喉和阴道拭子）时，宏基因组数据面临重大挑战。这些样本中宿主DNA污染严重，导致超过 90% 的测序读段与人类基因组对齐。这种污染不仅影响微生物群的解析，特别是对低丰度物种的鉴定，还会导致微生物组成的偏倚观察。此外，当宿主为人类时，隐私问题尤为重要，进一步凸显了去除宿主污染的必要性。

尽管实验阶段（特别是在DNA提取过程中）已经采取措施减少宿主污染，但残留的宿主DNA仍然广泛存在于样本中。不同实验方案的去污染效果存在差异，并可能偏向于特定微生物类群的富集。即便宿主DNA经过显著去除，高生物量样本（如黏膜微生物群）在宏基因组测序数据中仍可能保留高达 90% 的宿主污染。这主要是由于动物细胞与微生物在细胞和基因组大小上的差异。同样，在低生物量样本（如内生微生物）中，即使经过微生物富集，宿主DNA读序的比例仍可能高达70%。如此高水平的污染往往需要更深度的测序，以确保足够的微生物读段被捕获。然而，对这些非目标宿主DNA进行测序并随后从大规模 NGS（新一代测序）数据集中计算去除，既浪费资源又耗时，不仅影响下游分析的准确性，还消耗宝贵的研究时间和计算资源。因此，在数据分析阶段开发一种准确高效的宿主污染去除工具至关重要。

我们检索了包含“metagenome”和“microbiome”关键词的 2,853 篇文献，发现 57.94% 的研究涉及宿主污染去除，并且从 2015 年到 2024 年呈现出明显的增长趋势。由于缺乏统一的宿主去污染软件选择标准，目前已有 51 种不同的工具被使用，其中许多工具依赖于比对和 k-mer 策略。在这些工具中，10 种工具表现出较高的流行度，主要采用两种策略：比对（alignment-based）和 k-mer 方法。比对方法的软件，如 Bowtie2 和 BWA，通过将测序读段比对到参考基因组进行宿主污染去除。而基于 k-mer 方法的软件（如 Kraken2 和 KMCP）则通过在参考数据库中寻找读段的小读段（k-mers）的精确匹配来去除宿主序列。此外，一些宿主污染去除流程整合了这些模块，例如 DeconSeq 集成了 BWA 的改进版本，而 KneadData 则整合了 Bowtie2。近年来，一些新工具（如 Hostile 和 HoCoRT）被开发出来，以提高宿主去污染的准确性。尽管已有研究探索不同宿主 DNA 含量对微生物组分析的影响，但去除宿主 DNA 污染对生物信息学下游分析及微生物基因组组装的影响仍不明确。

本研究将通过去除宿主污染来比较宏基因组测序在微生物组分析中的效率，并系统评估当前最先进计算工具在宿主 DNA 去污染方面的准确性和速度。研究结果将为研究人员提供指导，帮助其合理选择适用于不同宏基因组数据集的工具。

在本研究中，我们使用 CAMISIM 模拟了三组数据（S1、S2、S3），其中包含 90% 的宿主污染（表 S1-2）。这些数据包括来自 30 种微生物的序列，每种以相等丰度表示，并额外加入了来自 Homo sapiens（GRCh38）的人类序列。我们使用近年来广泛应用的宿主去污染软件 KneadData，从原始数据（Raw）中去除宿主污染，得到去宿主数据（Remove），同时以 30 种微生物组数据作为阴性对照（Microbiome）。我们的宏基因组分析涵盖关键步骤，包括物种组成分析、多样性分析、功能分析以及宏基因组组装基因组（MAG）评估（图 1A）。基于这些模拟数据，我们评估了宿主去污染对宏基因组分析的影响，重点关注内存使用情况、处理时间及其对结果准确性的影响（图1A）。

在整个分析过程中，原始数据（Raw）、去宿主数据（Remove）和微生物组数据（Microbiome）在高内存需求步骤（如超过 100 GB 内存的物种分类注释（Kraken2）、去冗余分析（dRep）和 MAG 注释（GTDBtk））中，内存使用量未见显著差异。然而，与原始数据相比，去除宿主序列的数据显著缩短了下游分析的运行时间（图 1B）。具体而言，宿主去污染后的数据，在基因组分箱（MetaWRAP）中的处理时间缩短了 5.98 倍，在功能注释（HUMAnN3）中缩短了 7.63 倍，在组装（MEGAHIT）中缩短了 20.55 倍。去宿主数据的平均处理时间为 139.14 分钟，而原始数据的处理时间为 832.64 分钟；HUMAnN3 处理去宿主数据需要 308.92 分钟，而原始数据则需 2357.95 分钟；MEGAHIT 处理去宿主数据需 106.59 分钟，而原始数据需 2190.27 分钟。此外，阴性对照数据（Microbiome）的内存使用量和时间消耗与去宿主数据相似。

相比于微生物组数据（Microbiome），原始数据（Raw）改变了微生物群落的相对丰度，而去宿主数据（Remove）显示出与 Kraken2 注释的物种组成相似的模式（图 1C）。有趣的是，在去除原始数据中的脊索动物门（Chordata，即 Homo sapiens 所属门）后（de-chordata 组），其余分类群与微生物组数据的分类群组成相似（图 1C）。微生物组数据与去宿主数据在物种丰富度指数上没有显著差异，而原始数据的丰富度指数显著低于两者（图 S2A）。主坐标分析（PCoA）可视化了群落组成的变化，结果显示 PCoA 第一轴解释了 100% 的总体变异，表明样本组之间具有低维度性和明显的分离性。具体而言，原始数据（Raw）与微生物组数据（Microbiome）及去宿主数据（Remove）在 PCo1 轴上明显分开（图 S2B）。这一发现强调了宿主去污染对于揭示潜在微生物群落结构的有效性。

尽管模拟数据包含 30 种微生物，每种微生物均有 10 万条序列，但最终仅获得 14 个 MAG。微生物组数据、原始数据和去宿主数据在 MAG 完整度和污染率方面无显著差异（图 S2C）。然而，相比原始数据，微生物组数据和去宿主数据中获得的 MAG 数量明显更多，唯一的例外是 Bifidobacterium breve，该菌种在所有三组原始数据（S1、S2、S3）中均被检测到，而在微生物组数据和去宿主数据中仅在 S2 组检测到（图 1D）。接下来，我们将去宿主数据的基因本体（GO）术语与微生物组数据进行比较，发现去宿主数据与微生物组数据之间的 GO 术语相关性更高，而原始数据与微生物组数据的相关性较低（图 1E，图 S2D）。这一结果表明，宿主去污染过程能够提高基因功能注释的准确性。

图 1. 宿主污染消耗了额外的计算资源，并影响了宏基因组分析结果的准确性

为了进一步比较不同软件工具在宿主去除方面的差异，我们使用 1080 组模拟宏基因组数据集进行分析，其中包括单菌种数据集（SinBac）和合成群落数据集（SynCom），覆盖 10 Gbps、30 Gbps 和 60 Gbps 三种数据规模。这些数据分别针对人类（Homo sapiens）和水稻（Oryza sativa indica）宿主进行了模拟，每种宿主分别包含 90%、50% 和 10% 的宿主污染水平（图 2A，更多细节见方法部分）。

为了方便描述，我们对不同的数据集进行了命名，例如 SinBac10-1 代表一个 10 Gbps 数据集，其中 90% 读长来自宿主基因组，10% 读长来自单菌种基因组。同样，SynCom30-2 代表一个 30 Gbps 数据集，其中 50% 读长来自宿主基因组，50% 读长来自合成群落基因组。而 SynCom60-3 代表一个 60 Gbps 数据集，其中 10% 读长来自宿主基因组，90% 读长来自合成群落基因组（图 2A）。

基于这些模拟数据，我们对六种现有的宿主去除工具的计算资源消耗和宿主污染去除性能进行了比较。这些工具包括基于比对的方法（KneadData、Bowtie2 和 BWA），以及基于 k-mer 策略的方法（KMCP、Kraken2 和 KrakenUniq）。

在进行宿主去除前，宿主参考基因组的索引至关重要。本研究中，我们分别构建了人类（GRCh38）和水稻（GWHBFPX00000000）的参考基因组，大小分别约为 3.1 Gbps 和 373.8 Mbps。Kraken2 在索引这两种基因组时计算资源消耗最低，仅使用 0.3 GB 内存，并在 6.94 分钟内完成了人类基因组的自定义数据库构建（图 2B）。相比之下，其余五种工具平均需要 18.05 GB 内存，并耗时 117.98 分钟。

接下来，我们比较了六种软件在宿主污染去除过程中的资源消耗情况。总体而言，在所有模拟数据集中，Bowtie2（1.95 GB (0.410, 3.42)）和 Kraken2（2.47 GB (0.710, 4.12)）分别是比对方法和 k-mer 方法中最大内存使用量最低的工具（图 2C，图 S3）。相比之下，其他四种工具的内存需求更高，BWA 需要 3.995 GB (1.40, 6.74)，KneadData 需要 15.17 GB (6.47, 30.27)，KMCP 需要 14.45 GB (4.27, 25.110)，KrakenUniq 需要 22.410 GB (11.13, 33.67)。

在数据规模方面，处理 60 Gbps 数据的资源消耗明显高于 10 Gbps 数据，特别是在 KneadData 和 Kraken2 中。此外，宿主类型也显著影响所有软件的资源消耗（表 S3-1），人类数据的资源需求显著高于水稻数据（P < 0.05）。然而，除 KrakenUniq 之外，不同类型的微生物群落对资源消耗没有显著影响（表 S3-1）。

在运行时间方面，k-mer 方法（KMCP 156.10 分钟 (90.47, 231.16)，Kraken2 29.34 分钟 (13.42, 55.45)，KrakenUniq 59.23 分钟 (26.66, 98.14)）明显快于比对方法（BWA 582.26 分钟 (300.15, 1065.64)，Bowtie2 209.00 分钟 (111.66, 512.61)，KneadData 501.38 分钟 (287.41, 1177.06)），其中 Kraken2 的执行时间显著短于其他工具（图 2E，表 S3-2）。微生物群落的多样性对所有六种软件的处理时间没有影响（表 S3-2），但宏基因组数据的规模显著影响了执行时间（图 2D，2E），这凸显了在大规模宏基因组数据处理中使用高效工具的重要性。此外，相比水稻数据，人类数据的处理时间更长，这表明宿主基因组的复杂性会增加计算时间（图 S4，表 S3）。值得注意的是，宿主基因组污染比例较高时，BWA、Bowtie2、KneadData 和 KrakenUniq 的运行速度显著降低（图 S4，表 S3-2）。例如，处理一个 60 Gbps、90% 宿主污染的人类宏基因组数据集时，与 10% 污染的相同数据集相比，KrakenUniq 的运行时间增加了 1.35 倍，BWA 增加了 2.59 倍，KneadData 增加了 5.36 倍，Bowtie2 增加了 6.76 倍。这表明这些工具可能不太适用于宿主污染比例较高的数据集。然而，这种显著的运行时间下降在 KMCP 和 Kraken2 中并未观察到。

图 2. 基准测试评估六种宿主去除软件在模拟人类和水稻宏基因组数据上的计算资源消耗

我们使用四个指标（准确性、召回率、精确度和F1得分）来评估六种宿主去污软件在基于图2A规则生成的1080个模拟数据上的宿主去污准确性表现。我们观察到六种软件在准确性、召回率、精确度和F1得分上存在显著差异（P < 0.05）（表S3-4、S3-5、S3-6、S3-7）。在准确性方面，基于比对的软件（BWA，0.9989（0.9966，0.9998）；Bowtie2，0.9997（0.9988，0.9998）；KneadData，0.9997（0.9989，0.9998））表现优于基于k-mer的软件（KMCP，0.8947（0.8133，0.9748）；Kraken2，0.9891（0.9832，0.9974）），但KrakenUniq（0.9998（0.9994，0.9999））始终表现出高效且稳定的性能（图3A）。然而，基于比对的软件表现出较低的精确度（BWA，0.9980（0.9853，0.9996）；Bowtie2，0.9999（0.9999，0.9999）；KneadData，0.9981（0.9971，0.9998）），可能导致宿主基因组相关的假阳性增加。这意味着一些微生物读取可能被错误地映射为宿主基因组的一部分，并随后被作为污染物去除。相反，基于k-mer的软件（KMCP，0.7686（0.7477，0.7925）；Kraken2，0.9787（0.9787，0.9823）；KrakenUniq，0.9999（0.9999，1））表现出较低的召回率，导致宿主基因组相关的假阴性增多。这意味着一些宿主读取可能被错误地未映射，从而在后续分析中保留了一些宿主污染（图3B）。对于F1得分（2 * 精确度 * 召回率 /（精确度 + 召回率）），微生物类型、宿主类型以及宿主基因组的比例都对这些工具的表现产生了影响。值得注意的是，BWA、KneadData和KrakenUniq在人体数据集上的表现明显优于水稻数据集。相反，Bowtie2在水稻数据集上表现出色（图3C）。

我们比较了使用上述六种软件进行宿主去污后，合成群落（SynCom）数据集的组成（图2A）。在BWA、KneadData和KrakenUniq中，类群Magnoliopsida和Mammalia的相对丰度（即对数转化的相对丰度）值较低（图3D）。在这些工具中，BWA和KneadData，作为基于比对的工具，在去除宿主污染方面表现优于KrakenUniq，后者作为基于k-mer的工具，效果较差。除了KMCP之外，所有工具都识别出一些属于放线菌、梭状芽胞杆菌、负膜菌等类群的宿主污染物，并将其从原始数据中去除（图3D）。为了对比资源消耗和宿主去污性能，我们对所有数据进行了最小-最大归一化。这使得我们能够比较不同软件在模拟数据集中的计算效率和宿主污染去除效果（图3E，图S5）。根据归一化数据的总结，Kraken2在高宿主污染水平（90%）下表现出显著的优势（P < 0.05），在人类（7.8093（7.7236，7.8215））和水稻（7.8268（7.7109，7.8330））数据集中的表现均为优异。此外，在比较基于比对的软件在高宿主污染水平（90%）下的宿主去除表现时，重点关注归一化准确性（NA）、归一化精确度（NP）、归一化召回率（NR）和归一化F1得分（NF1），KneadData在人类数据集中的表现显著优越（P < 0.05）（3.9996（3.9996，3.9997））。

图3. 评估各软件在宿主污染去除中的准确性

接下来，我们评估了缺少宿主参考基因组对现有宿主去污染软件效果的影响。选取了三种水稻物种——粳稻（Oryza sativa japonica, GWHBFOO00000000, Osj）、籼稻（Oryza sativa indica, GWHBFPT00000000, Osi）和野生稻（Oryza rufipogon, GWHBFHN00000000, Or）作为宿主宏基因组读取的资源来生成模拟数据集。籼稻（Oryza sativa indica, GWHBFPX00000000, Refer）被选为参考基因组（图4A）。三种物种（Osj、Osi、Or）与参考基因组之间的平均核苷酸相似性（ANI）显示，Osi与参考基因组的相似性最高（99.01%），其次是Osj（97.94%）和Or（97.49%）（图4B）。如前所述，基于单一细菌（SinBac）生成了模拟的宏基因组数据，产生了三个数据集（OsjSinBac、OsiSinBac、OrSinBac），其中包含不同程度的宿主DNA污染（10%、50%和90%），每个数据集大小为10 Gbps。随后，使用参考基因组构建了六种工具的索引数据库，我们比较了在这些条件下宿主去污染工具（BWA、Bowtie2、KneadData、KMCP、Kraken2和KrakenUniq）的性能（图4A）。

在时间和内存消耗方面，KneadData和KrakenUniq使用了更多的内存（图S6A），而比对软件所需的时间比k-mer软件更长（图S6B），如前所述。与参考宏基因组数据相比，处理OsjSinBac、OsiSinBac、OrSinBac数据时，所有软件都需要显著更多的时间（图S6B）。值得注意的是，缺乏紧密比对且准确的宿主参考基因组对所有工具的去污染性能产生了负面影响（图4C）。具体而言，代表OsjSinBac、OsiSinBac和OrSinBac的数据集的准确性、召回率和F1得分显著低于与参考基因组对齐的参考数据。精确度没有显著差异，表明去污染的数据集仍然包含一些残留的宿主读取。

然后，我们比较了缺乏准确宿主参考基因组时所有软件的综合性能。基于运行过程中的资源消耗和宿主去污染指标（准确性、精确度、召回率和F1得分）的归一化数据总结。所有软件在refer数据上表现更好（5.41（5.15，5.91）），相比之下，在OsjSinBac（4.93（4.32，5.50））、OsiSinBac（4.54（3.91，5.35））和OrSinBac（4.64（4.08，5.19））数据集上的表现较差，但Kraken2和KrakenUniq在这些数据集之间的差异较小（图S6C）。此外，虽然不同工具在参考基因组的高宿主污染样本上表现更好，但在没有参考基因组时这种优势并不明显。在90%宿主污染的情况下，所有工具的准确性都有显著下降（图4D，表S4-1），这突显了宿主参考基因组在高污染宏基因组数据中的重要性。

图4. 缺乏宿主参考基因组对宿主去污染工具性能的影响

宿主污染在宏基因组数据分析中仍然是一个重要的挑战，特别是对于来源于复杂环境或超高深度测序的数据集。在我们的研究中，我们观察到57.94%的文献在分析数据之前消除了宿主污染，而其余的研究则选择直接分析数据，避免在去除宿主污染的过程中可能牺牲宝贵的微生物组信息。一些软件被开发用于去除宿主基因组污染，但很少有研究对它们之间的区别进行了比较。为了加速宏基因组数据的可重复性、可比性和标准化，理解这些不同工具的影响，广泛地被称为下游分析，是非常重要的。而且，功能注释、物种注释、序列拼接和分箱在计算资源消耗方面有显著影响，因此去除宿主污染，特别是在超高深度测序数据中的必要性更为突出。

我们生成了一个模拟数据集，其中包含30个微生物物种的10万条序列和27百万条人类（Homo sapiens）序列。在微生物物种中，六个属于放线菌门（Actinomycetota），15个属于芽孢菌门（Bacillota），六个属于拟杆菌门（Bacteroidota），各一个属于伪单胞菌门（Pseudomonadota）、热硫还原细菌门（Thermodesulfobacteriota）和变形菌门（Verrucomicrobiota）。然而，相对丰度的比例并不相等，并且有两个未包括在模拟中的门（Streptophyta和Campylobacterota）被注释出来。这种差异可能归因于Kraken2注释过程的限制，尽管我们使用了最全面的数据库（PlusPFP），以及模拟数据的异质性。此外，尽管每个模拟数据包含30个微生物群体，但我们只获得了14个MAGs（宏基因组单元）。值得注意的是，Bifidobacterium breve（B. breve）仅在使用严格的分箱细化参数（-c 50，-x 10）时从原始样本中恢复出来。通过调整MetaWRAP细化参数为-c 0 -x 100，在微生物组和去除数据集中恢复了更多的B. breve MAGs，并改善了完整性和污染度指标。这个发现突出了在严格优先考虑MAG质量和在更宽松的设置下增强微生物多样性恢复之间的权衡。它强调了需要进一步完善宏基因组分析软件，以在这些目标之间取得更好的平衡，尤其是在复杂的数据集中。尽管如此，宿主去污染的重要性仍然在基因功能注释方面得到了体现。同时，微生物组（负对照）数据与原始数据或去除数据之间的显著α多样性指数差异，突出了建立宿主去除“金标准”的重要性。

值得注意的是，我们观察到KneadData在高（90%）宿主污染数据集中的准确性、精确度和F1得分显著高于低（10%）宿主污染数据集。其他对齐软件也观察到类似的趋势。这一现象可能归因于两个因素。首先，我们的模拟数据源自宿主基因组，作为索引参考数据库，因此在高宿主污染数据（90%）中，宿主读取的比例较高，这些宿主读取会对齐到宿主参考基因组中（在对齐软件如BWA、Bowtie2、KneadData中）。其次，来自微生物组的序列可能很难与参考基因组区分，导致在低宿主污染（10%）数据集中，更多的微生物读取被丢弃。

我们同意，比较不同软件中参数设置的影响是优化宿主污染去除的重要考虑因素。不同工具有不同的默认参数，这些参数可以显著影响去污染的准确性和计算效率。例如，一些研究比较了Bowtie2和Kraken2在各种参数设置下的性能，突出了计算资源和宿主去除效果之间的权衡。未来的工作应集中在系统地比较多种工具的参数设置，以便提供有关如何平衡这些因素的有价值见解，特别是在处理大数据集或高宿主污染时。这种比较方法将有助于标准化工作流程，提高可重复性，并增强宿主去污染过程在宏基因组研究中的效率。

此外，当宿主基因组不可用时，预测宿主基因组的挑战也需要考虑，可能需要到属、科或目水平。有趣的是，尽管Osi与参考基因组的ANI比Osj更接近，但OsiSinBac的去污染性能并不如OsjSinBac那样有效。而且，高宿主污染导致所有软件的准确性显著下降，强调了宿主参考基因组在高污染宏基因组数据中的重要性。缺乏准确的宿主参考基因组可能导致残留宿主序列，从而降低随后的功能注释的精度。一个可能的解决方案是结合基于对齐和k-mer的方法，以更准确地去除宿主污染。然而，一个独特的挑战仍然存在，即如何区分由水平基因转移（HGT）而非宿主污染引起的微生物序列。微生物组与宿主基因组之间的HGT，通常涉及可移动遗传元件（MGEs），在微生物适应多样化环境中起着至关重要的作用。宏基因组测序，特别是短读技术，面临着准确识别这些水平转移基因区域的重大困难。未来，人工智能算法和长读测序的进展可能不仅有助于克服宿主污染问题，还可能解决细菌的水平基因转移问题。

总之，宿主去污染不仅加速了下游分析，还提高了基因功能注释的准确性，特别是在超高深度测序数据中。而且，每种工具（BWA、Bowtie2、KneadData、Kraken2、KMCP和KrakenUniq）都具有独特的优势，可以根据研究的具体需求加以利用。简而言之，Bowtie2和KneadData提供了更准确的去除能力，尽管需要更高的计算资源。Kraken2和KrakenUniq则提供了快速且用户友好的解决方案，而KMCP能够保留更多低丰度的分类群。在缺乏参考基因组的情况下，BWA和KneadData在对齐软件中受到的影响较小，而Kraken2和KrakenUniq在k-mer软件中受到的影响较小。

理解速度、准确性和计算资源之间的权衡，对于选择最合适的宿主DNA去除工具至关重要。随着研究日益关注宿主污染对微生物组注释的影响，尤其是在低丰度分类群方面[27]，本研究提供了一个全面的评估，为工具和方法的改进奠定了基础。最终，这些进展将使研究人员能够从复杂的宏基因组数据集中提取有意义的生物学见解。

2024年4月18日，在Web of Science核心合集数据库中进行了文献搜索，使用的搜索词为“metagenome”和“microbiome”。仅使用研究文章收集关于软件使用的信息。我们排除了仅关注扩增子测序数据、长读长宏基因组数据、环境样本或食品样本的文献。随后，我们汇总了提及宿主污染去除的文献比例和每个软件的文献数量。

我们模拟了三组数据（S1、S2、S3），宿主污染率为90%，包括来自30个微生物物种和人类基因组（GRCh38）的基因组。30个微生物物种是基于先前发布的人类相关微生物群落随机选择的。我们使用CAMISIM生成了30个微生物基因组的宏基因组数据，模拟了成对端读数（PE150），每个微生物物种有三次重复（number_of_samples = 3）。每个数据集包含300万条微生物读数，每个微生物物种0.1百万条读数，外加2700万条来自人类基因组的读数。数据在混合前已按相应标签标记。在微生物物种中，六个属于放线菌门，15个属于芽孢菌门，六个属于拟杆菌门，另外一个分别属于假单胞菌门、热硫还原菌门和厚壁菌门。为了去除原始数据中的宿主污染（Raw），我们使用了专为此目的设计的生物信息学工具KneadData。KneadData处理的输出构成了宿主去除数据（Remove）。30个微生物群体作为负对照（Microbiome）来评估宿主污染去除过程的准确性和有效性。

为了比较直接数据（Raw数据）与去除宿主污染数据（Remove数据）的分析差异，我们选择了10 Gbps的合成群落数据集，包含不同宿主污染水平的人类数据。资源消耗按以下描述进行了测试，宏基因组分析直接参考了EasyMetagenome 1.10流程的步骤。简而言之，分类学分析使用Kraken2和PlusPFP数据库进行，使用Bracken计算相对丰度。功能分析通过HUMAnN3使用Uniref90基因家族进行。使用Megahit 1.0进行宏基因组数据的组装，宏基因组分箱和分箱优化使用MetaWRAP进行。MetaWRAP优化旨在提高MAG分箱的质量，利用参数-c 50 -x 10，仅保留完整性大于50%且污染小于10%的分箱。使用dRep v2.6.2去除宏基因组组装基因组（MAGs）的冗余。使用GTDBtk v2.3.2对MAGs进行注释，并使用CheckM2评估其质量。基因预测使用Prodigal v2.6.3，基因聚类使用CD-HIT v4.8.1，基因定量使用salmon v1.8.0，基因注释使用emapper v2.1.6。然后，使用R 4.2.3进行alpha多样性和beta多样性分析，分析方法参考EasyAmplicon。宏基因组组装基因组（MAGs）的完整性和污染率是正常测量数据，因此以中位数（P25，P75）表示。我们还使用eggnog-mapper对GO术语进行了注释，并计算了它们与微生物组数据的相关性，仅保留来自原始数据的微生物组，使用Spearman进行计算。

我们选择了在人类学和医学上具有重要经济和医学意义且基因组已知的稻米和人类作为研究对象。选择了六个工具进行分析，其中三个是基于比对的软件（BWA、Bowtie2和KneadData），其余是基于k-mer的方法（KMCP、Kraken2、KrakenUniq）。使用CAMISIM生成模拟数据集，采用默认或作者推荐的参数进行分析。为了确保可比性和可靠性，每个数据集包含五个重复。这些数据集涵盖了不同的数据大小（10 Gbps、30 Gbps、60 Gbps）、不同宿主DNA污染水平（90%、50%、10%）和多样的微生物复杂性（单细菌SinBac或合成群落SynCom）来自人类和稻米样本。每个数据集每种条件下有五个重复。对于稻米SynCom，我们随机选择了14个常见的稻米相关物种。人类SynCom由35个物种组成，基于先前发布的人类相关微生物群落。对于SinBac模拟，我们在CAMISIM中使用默认参数生成成对端读数（PE150），对于宿主和单一细菌基因组，进行五次重复模拟（number_of_samples = 5）。然后将宿主和微生物基因组的读数按不同的宿主DNA污染比例混合。在SynCom模拟中，我们采用与SinBac相同的方法生成宿主的宏基因组数据，但对于微生物宏基因组数据，我们使用CAMISIM中的差异模式。随后，根据不同宿主DNA污染水平混合这些数据。物种的分类信息及其基因组ID提供在表中。为了确保全面评估，我们使用CAMISIM生成了1080个模拟数据集。这些数据集包括三种不同的大小（10 Gbps、30 Gbps和60 Gbps），代表简单（SinBac）和复杂（SynCom）微生物群落。模拟分别针对人类（Homo sapiens, GRCh38）和籼稻（Oryza sativa indica, GWHBFPX00000000）宿主进行，每个宿主具有三种宿主污染水平（10%、50%、90%），以便在不同条件下深入探讨宿主基因组污染去除的影响。对于含有多个染色体的物种，我们下载了所有的染色体并模拟信息，所有参考和fasta信息附在表中。在这里，为了便于描述，我们为各种数据集分配了缩写。例如，SinBac10-1表示一个10 Gbps的数据集，90%的读数来自宿主基因组，10%的读数来自单一细菌基因组。同样，SynCom30-2表示一个30 Gbps的数据集，50%的读数来自宿主基因组，50%的读数来自合成群落基因组。基于这些模拟数据，我们评估了去除宿主污染对微生物组分析的影响。

为了评估宿主去污染工具在属水平的表现，我们分别使用了三种稻米物种和简单微生物群落：粳稻（Oryza sativa japonica，GWHBFOO00000000，Osj），籼稻（Oryza sativa indica，GWHBFPT00000000，Osi）和野生稻（Oryza rufipogon，GWHBFHN00000000，Or）。对于每个物种，我们通过将宿主基因组与单细菌类群基因组结合生成数据集，得到三个独立的数据集：Or与单细菌类群（OrSinBac）、Osj与单细菌类群（OsjSinBac）和Osi与单细菌类群（OsiSinBac）。每个数据集使用CAMISIM生成，具有不同的宿主DNA污染水平（10%、50%、90%），并且为10 Gbps数据集设置了五次重复。

Oryza sativa indica作为参考基因组（GWHBFPX00000000，Refer），为六个工具（BWA、Bowtie2、KneadData、KMCP、Kraken2、KrakenUniq）创建了索引，具体如上所述。详细的fasta信息见表。这些模拟数据集被用于比较六个工具在属水平去除宿主污染的表现。比较涉及评估计算资源和宿主去污染的有效性。基因组相似度使用fastANI v1.34计算，并通过以下性能指标评估去除宿主污染对微生物组的影响。

所有软件在构建数据库和运行过程时均使用了八个线程，针对模拟数据集进行分析。我们对资源消耗进行了比较分析，重点关注最大内存使用量和处理时间。此外，我们还评估了它们的性能指标，包括真阳性（TP）、真阴性（TN）、假阳性（FP）和假阴性（FN）。计算了精度、召回率和F1分数，公式如下：

准确率 = (TP + TN) / (TP + FN + FP + TN)

精确度 = TP / (TP + FP)

召回率 = TP / (TP + FN)

F1分数 = 2 * 精确度 * 召回率 / (精确度 + 召回率)

然后，所有生物信息学分析都在R 4.2.3中进行。我们对所有数据进行了正态性和均匀性检验。对于正态分布数据，采用均值±标准误（SE）表示；对于非正态分布数据，采用中位数（P25, P75）表示。对于方差齐性且正态分布的两组数据比较，使用配对t检验；非正态分布或方差不齐的数据显示，使用非参数的Wilcoxon检验。在多组比较中，正态分布且方差齐的数据采用ANOVA分析，而非正态分布或方差不齐的数据则采用非参数的Kruskal-Wallis检验。对每组数据进行了Bonferroni事后检验，以分析不同数据集之间的差异。对于平均物种组成图，在处理包含零的数据时，我们采用将一个正常常数（0.001）加到所有叶片性状值上的简便方法，然后取结果的对数平均相对丰度。对于每个软件中宿主去污染资源消耗和性能的综合比较，我们对所有数据进行了最小-最大标准化，这是一种对原始数据进行线性转换的技术。数据可视化使用了ggplot2包，P ≤ 0.05被认为具有统计学意义。

1.数据分析流程

项目名称：HostPurge-DownstreamAnalysis

项目主页：https://github.com/YunyunGao374/HostPurge/blob/main/0HostDecontaminationImpactiononDownstreamAnalysis.sh

操作系统：Linux

编程语言：Bash, R

其他需求：Environment Modules, Conda

2.模拟数据获取流程

项目名称：HostPurge-SoftwareComparison

项目主页：https://github.com/YunyunGao374/HostPurge/blob/main/1HostDecontaminationSoftwareComparison.sh

操作系统：Linux

编程语言：Bash, R

其他需求：Environment Modules, Conda

3.文章图片代码来源

项目名称：HostPurge-paper-figures

项目主页：https://github.com/YunyunGao374/HostPurge

操作系统：Windows, MacOS, Linux

编程语言：R

其他需求：N/A

许可证：GNU GPL 3

北京林业大学讲师、中国农业科学院深圳农业基因组研究所刘永鑫组博士后高云云和中国农业科学院深圳农业基因组研究所在读博士研究生罗豪为 本文的第一作者，中国农科院基因组所食品中心研究员，微生物组与营养健康团队首席，iMeta执行主编，宏基因组公众号创始人刘永鑫为本文的通讯作者。