宏基因组：从海量基因组到物种目录的构建、注释与亚型解析

天意生信云

发布于 2025-11-20 16:29:17

170

文章被收录于专栏：天意生信俱乐部天意生信俱乐部

经过前面的处理，我们手中握有一个包含近20万个经过严格质控的基因组（MAGs + 参考基因组）的庞大集合。然而，这片数据的海洋充满了冗余，且个体的身份仍然未知。

本教程将指导你完成最后三项、也是最具价值的任务：

物种聚类与去冗余: 采用预分桶结合 dRep 的高效策略，将数十万基因组精确聚类为非冗余的物种级别基因组集 (SGBs)。
权威分类学注释: 使用 GTDB-Tk 为每个 SGB 赋予当前最前沿的分类学名称。
群内精细结构解析: 利用 Mash 和层次聚类深入 SGB 内部，判定亚种/谱系的存在。

最佳实践流程概览

直接对近20万个基因组进行 all-vs-all 的 ANI 比较是不现实的。我们将采用一种更智能的多级策略：

快速预分桶 (Pre-binning): 使用超快速的 Mash 将整个库粗略地划分为数千个“桶”，确保同一物种大概率落入同一桶。
桶内精确去冗余: 对每个小规模的桶，独立运行 dRep 进行高精度的 95% ANI 物种聚类。
全局代表合并: 合并所有桶的代表基因组，再进行一次最终的 dRep 去冗余，消除分桶边界的重复。
分类学注释与亚群判定: 对最终的 SGB 代表集进行 GTDB-Tk 注释和群内距离分析。

准备工作：环境部署

conda activate metagenome

# 一次性安装所有需要的工具
# dRep, Mash, GTDB-Tk, 以及 R 环境 (用于层次聚类)
conda install -c bioconda -c conda-forge drep mash gtdb-tk r-fastcluster -y

# GTDB-Tk 需要下载其庞大的数据库，首次使用前必须运行
# 这将需要很长时间和超过 100GB 的磁盘空间
download-gtdb-tk.sh

步骤一：Mash 快速预分桶

软件作用

Mash 基于 MinHash 算法，可以极快地估算基因组之间的距离（Mash distance），这个距离与 ANI 高度相关。我们用它来执行第一轮的粗略聚类，避免 dRep 直接面对海量数据。

软件用法

这是一个两步过程：计算距离矩阵，然后基于矩阵进行聚类。

# 1. 为所有基因组创建 sketch (签名文件)
mash sketch -p 48 -o all_genomes_sketch FINAL_GENOME_DATABASE/*.fa

# 2. 计算两两距离 (输出一个制表符分隔的文件)
mash dist -p 48 all_genomes_sketch.msh all_genomes_sketch.msh > all_genomes.mash.dist.tsv

# 3. 编写脚本进行分桶
# 这一步需要一个简单的脚本 (Python/Perl/R)
# 逻辑：读取 all_genomes.mash.dist.tsv，使用一个距离阈值 (如 0.1，约对应 90% ANI)
# 进行单链接聚类 (single-linkage clustering)，将聚在一起的基因组分到一个“桶”
# 输出：多个文件，每个文件是一个桶的基因组列表 (e.g., bin_1.txt, bin_2.txt, ...)

步骤二：dRep 桶内精确去冗余与全局合并

软件作用

dRep 基于更精确的 fastANI 算法，在每个桶内部执行 95% ANI 的物种聚类，并根据我们预设的质量评分，选出每个物种簇（SGB）中唯一的、最佳的代表。

软件用法

桶内并行 dRep

需要编写一个循环或使用 GNU Parallel 来为上一步生成的每个 bin_*.txt 运行 dRep。

# 示例：处理一个桶 bin_1.txt
dRep dereplicate drep_bin_1 \
    -g bin_1.txt \
    -p 24 \
    --S_algorithm fastANI \
    -pa 0.9 \
    -sa 0.95 \
    -nc 0.30 \
    --genomeInfo all_genome_quality.csv # 提供包含所有基因组质量评分的文件

全局合并

# a. 将所有桶的 dRep 输出目录中的 dereplicated_genomes/*.fa 收集到一个新目录 ALL_REPS/

# b. 对这个已经显著缩小的集合再运行一次 dRep
dRep dereplicate drep_FINAL \
    -g ALL_REPS/*.fa \
    -p 48 \
    --S_algorithm fastANI \
    -sa 0.95 \
    -nc 0.30 \
    --genomeInfo all_genome_quality.csv

最终，drep_FINAL/dereplicated_genomes/ 目录中的基因组就是我们最终的、非冗余的 SGB 代表集。

步骤三：GTDB-Tk 分类学注释

软件作用

GTDB-Tk 是目前微生物分类学的黄金标准。它基于全基因组系统发育分析，将基因组放置在 GTDB (Genome Taxonomy Database) 的参考树上，从而提供一个稳定、动态更新的分类学注释。

软件用法

# 对最终的 SGB 代表集进行注释
gtdbtk classify_wf \
    --genome_dir drep_FINAL/dereplicated_genomes/ \
    --out_dir gtdbtk_output \
    --cpus 48

# 结果解读:
# 核心输出是 gtdbtk_output/gtdbtk.bac120.summary.tsv (细菌) 和
# gtdbtk.ar53.summary.tsv (古菌)，其中包含了每个 SGB 的完整七级分类 (域-门-纲-目-科-属-种)。

# 映射到 NCBI 分类 (可选但常用)
# 使用社区提供的脚本 (如 gtdb_to_ncbi_majority_vote.py)
# 可以将 GTDB 名称尽可能地映射回大家更熟悉的 NCBI 分类体系。

步骤四：Mash 亚群结构判定

软件作用

在定义了 SGB 之后，我们可能还想知道一个 SGB 内部是否包含多个亚种或谱系。我们可以回到原始的、属于同一个 SGB 的所有基因组，用 Mash 重新计算它们之间更精细的距离。

软件用法

以一个 SGB（例如 SGB_1）为例，该 SGB 包含10个基因组（从 dRep 的 Cdb.csv 文件中可知）。

# 1. 提取属于 SGB_1 的所有基因组文件
# 2. 对这10个基因组计算 Mash 距离
mash sketch -o SGB_1.msh SGB_1_genomes/*.fa
mash dist SGB_1.msh SGB_1.msh > SGB_1.dist.tsv

# 3. 在 R 中进行层次聚类
# R 脚本示例:
# library(fastcluster)
# dist_matrix <- as.dist(read.table("SGB_1.dist.tsv"))
# hc <- hclust(dist_matrix, method = "complete")
# plot(hc) # 可视化聚类树
#
# # 根据阈值判定亚群
# # Mash 距离 0.03 ≈ 97% ANI (亚种?)
# # Mash 距离 0.01 ≈ 99% ANI (谱系?)
# subclusters_97_ani <- cutree(hc, h = 0.03)
# subclusters_99_ani <- cutree(hc, h = 0.01)
# print(paste("在 97% ANI 水平下有", max(subclusters_97_ani), "个亚群"))

总结

至此，已经完成了宏基因组分析中从数据到知识的终极飞跃。通过这套高效、严谨的流程，不仅将一个庞杂的基因组集合提炼成了一个结构清晰、身份明确的 SGB 目录，还具备了深入探索物种内部微观多样性的能力。这个最终的 SGB 集合，是开启所有下游比较基因组学、系统发育学和生态学分析的坚实基石。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2025-10-29，如有侵权请联系 cloudcommunity@tencent.com 删除

数据