宏基因组 Binning 与 MAGs 提取¶
彭文强 | 2026届 | 生信分析工程师求职
1. 一句话说明¶
Binning 就是把宏基因组组装出的成千上万条 contig 按"谁和谁来自同一个菌"分堆，最终拼出单个细菌的近完整基因组（MAG），让你从"知道有哪些菌"进化到"知道这些菌能干什么"。
2. 为什么要学¶
你目前已经能做的（Kraken2 物种分类）¶
你的 T2D 项目已经用 Kraken2 做了 reads 级别的物种注释，能回答： - 这个样本里有哪些菌？ - 每种菌的相对丰度是多少？
Binning + MAGs 能额外回答的问题¶
| 问题 | 纯物种分类 | MAGs |
|---|---|---|
| 这个菌有哪些代谢通路？ | 只能推测（参考基因组的） | 直接从 MAG 注释得到 |
| T2D 患者的某个菌是否携带耐药基因？ | 无法确认 | MAG 里搜就知道 |
| 数据库里没有的新菌怎么办？ | Kraken2 分不出来，归为 unclassified | binning 照样能把它聚成一个 MAG |
| 同一种菌在不同患者里功能有差异吗？ | 看不出来 | 比较不同样本的 MAG，发现菌株级差异 |
| 菌与菌之间有没有代谢互作？ | 无法分析 | 通过 MAG 的功能注释推断代谢网络 |
白话总结：Kraken2 告诉你"菜市场里有什么菜"，MAGs 告诉你"每道菜的配方是什么"。你的 T2D 项目如果加上 MAGs 分析，就能从"T2D 患者肠道菌群组成变化"深入到"T2D 相关菌株的功能特征"，档次直接升一级。
3. 核心概念白话版¶
3.1 什么是 Binning（分箱/归箱）¶
想象场景：
你家搬家，所有书、衣服、厨具混在一起装了 100 个快递箱。
现在你要把它们分回"书是书、衣服是衣服、厨具是厨具"。
你会怎么分？
1. 看大小和颜色（书大多是长方形、厨具是金属色）→ 类似 GC 含量 + 四核苷酸频率
2. 看包裹单号（同一批寄的东西通常是一起的）→ 类似覆盖度
3. 综合判断，把"大概率来自同一个来源"的东西放到一个箱子里
Binning 做的就是这件事，只不过"东西"是 DNA 片段，"箱子"是不同的细菌基因组。
一句话定义：Binning 是利用序列特征（GC含量、四核苷酸频率）和丰度信息（覆盖度），把组装后的 contig 聚类到不同的"bin"中，每个 bin 理想情况下代表一个物种的基因组。
3.2 Contig vs Scaffold vs MAG¶
| 术语 | 白话解释 | 类比 |
|---|---|---|
| Read | 测序仪直接读出的短片段（150bp） | 一个拼图碎片 |
| Contig | 多个 read 拼接成的连续序列（几百到几万 bp） | 拼好的一小块拼图 |
| Scaffold | 多个 contig 用 paired-end 信息串起来（中间可能有 gap，用 N 填充） | 几小块拼图用线串起来，中间有缺口 |
| MAG | 一个 bin 里所有 contig 的集合，代表一个菌的基因组草图 | 一整个拼图（可能缺几块、可能混了别的拼图的碎片） |
3.3 Binning 用的三大特征¶
特征1：覆盖度（Coverage/Depth）¶
原理：同一个菌的所有 DNA 片段在测序时被测到的次数应该差不多。
白话：如果肠道里某个菌占 10%，那它的所有 contig 被 reads 覆盖的深度都应该在同一水平。
如果 contig_A 覆盖度 50x，contig_B 覆盖度 50x → 大概率来自同一个菌
如果 contig_A 覆盖度 50x，contig_C 覆盖度 5x → 大概率来自不同的菌
多样本的优势：
样本1 里 contig_A 和 contig_B 覆盖度都高 → 只能说可能同源
样本1 和 样本2 里 contig_A 和 contig_B 覆盖度变化趋势一致 → 更确定同源
（所以多样本联合 binning 效果更好！）
特征2：GC 含量¶
原理：同一个菌的 DNA 里 G+C 碱基的比例是比较稳定的。
白话：就像每个人说话有自己的口音——有的菌 GC 含量高（60%），有的低（30%）。
同一个菌产出的 contig，GC 含量应该差不多。
特征3：四核苷酸频率（Tetranucleotide Frequency, TNF）¶
原理：DNA 里相邻 4 个碱基的组合频率（AAAA, AAAT, AATG, ...共 256 种），
同一个物种的 TNF 模式非常稳定，不同物种之间差异明显。
白话：就像每个人的"写字笔迹"——同一个人写的不同段落笔迹相似，
不同人的笔迹区别很大。TNF 就是细菌 DNA 的"笔迹"。
3.4 完整性和污染度（CheckM 评估指标）¶
完整性（Completeness）：这个 bin 包含了多少"应该有的基因"？
- 通过检查单拷贝标记基因（每个菌都应该只有一份的基因，比如核糖体蛋白基因）
- 如果 104 个标记基因里找到了 95 个 → 完整性 ≈ 91.3%
- 白话：这个拼图你拼了多完整？100% = 全拼好了，50% = 缺了一半
污染度（Contamination）：这个 bin 里混入了多少"不该有的基因"？
- 如果单拷贝标记基因出现了重复 → 说明混入了别的菌的序列
- 标记基因有 5 个出现了 2 次 → 污染度 ≈ 5/104 ≈ 4.8%
- 白话：你的拼图里混进了多少块别人的碎片？0% = 完全纯净
3.5 物种级 MAG vs 菌株级 MAG¶
| 维度 | 物种级 MAG | 菌株级 MAG |
|---|---|---|
| ANI（平均核苷酸一致性） | 同物种：ANI > 95% | 同菌株：ANI > 99.9% |
| 数据需求 | 常规宏基因组（10-20 Gb/样本） | 深度测序或长读长测序 |
| 应用场景 | 物种功能注释 | 传播链追踪、菌株进化 |
| 你的项目 | 适合（T2D 肠道菌群功能分析） | 暂时不需要（数据量不够） |
4. 完整 MAGs 提取流程¶
流程总览¶
质控后的 Clean Reads
│
▼
① 组装（MEGAHIT）
│ 产出：contigs.fa
▼
② 比对回帖（BWA + samtools）
│ 产出：sorted.bam → 覆盖度信息
▼
③ Binning（MetaBAT2 / MaxBin2 / CONCOCT）
│ 产出：多个 bin（.fa 文件）
▼
④ Bin 精炼（DAS Tool）
│ 产出：精选的高质量 bin
▼
⑤ 质量评估（CheckM2）
│ 产出：完整性 + 污染度报告
▼
⑥ 物种注释（GTDB-Tk）
│ 产出：每个 MAG 的物种分类
▼
⑦ 功能注释（Prokka / DRAM）
│ 产出：基因预测 + 代谢通路
▼
最终成果：一组高质量 MAGs + 功能/物种注释
环境准备¶
# ============================================================
# 创建专用 conda 环境
# 白话：建一个独立的工具箱，把 binning 要用的工具都装进去
# ============================================================
# 创建环境（用 mamba 更快，mamba 是 conda 的加速版）
mamba create -n binning python=3.10 -y # 创建名为 binning 的环境
conda activate binning # 激活环境
# 安装组装工具
mamba install -c bioconda megahit -y # MEGAHIT：宏基因组组装器
# 安装比对工具
mamba install -c bioconda bwa-mem2 samtools -y # bwa-mem2：比对器（bwa升级版，更快）
# samtools：BAM 文件处理瑞士军刀
# 安装 binning 工具（三个主流工具都装上，后面用 DAS Tool 整合）
mamba install -c bioconda metabat2 maxbin2 concoct -y
# 安装 bin 精炼工具
mamba install -c bioconda das_tool -y # DAS Tool：整合多个 binner 的结果
# 安装质控工具
mamba install -c bioconda checkm2 -y # CheckM2：评估 bin 质量（比 CheckM1 更快）
# 安装注释工具
mamba install -c bioconda gtdbtk prokka -y # GTDB-Tk：物种分类
# Prokka：基因预测与注释
# GTDB-Tk 需要下载参考数据库（约 110 GB 解压后，提前下好）
# 白话：GTDB-Tk 要靠一个巨大的"菌种图鉴"来对照，不下载就没法分类
# 注意：运行 GTDB-Tk 时 pplacer 步骤需要约 210 GB 内存
download-db.sh # GTDB-Tk 自带的数据库下载脚本
# 或手动设置路径：
export GTDBTK_DATA_PATH=/path/to/gtdbtk_data/release220
4.1 组装（MEGAHIT）¶
# ============================================================
# 第一步：宏基因组组装
# 目的：把短 reads 拼成长 contigs（相当于把碎纸片拼成段落）
# 工具：MEGAHIT（内存友好，适合服务器资源有限的情况）
# ============================================================
# 设置变量（方便后续引用，不用每次都打完整路径）
SAMPLE="T2D_sample01" # 样本名
CLEAN_R1="clean_data/${SAMPLE}_R1.fq.gz" # 质控后的正向 reads
CLEAN_R2="clean_data/${SAMPLE}_R2.fq.gz" # 质控后的反向 reads
OUTDIR="assembly/${SAMPLE}" # 输出目录
# 运行 MEGAHIT
megahit \
-1 ${CLEAN_R1} \ # 正向 reads（R1）
-2 ${CLEAN_R2} \ # 反向 reads（R2）
-o ${OUTDIR} \ # 输出目录
--min-contig-len 1000 \ # 只保留 ≥1000bp 的 contig（太短的不好 binning）
--k-list 21,29,39,59,79,99,119,141 \ # 多个 k-mer 大小（MEGAHIT 会合并最优结果）
-t 16 \ # 用 16 个 CPU 线程
--memory 0.5 # 最多用 50% 内存（防止服务器爆内存被 kill）
# 白话解释参数：
# --min-contig-len 1000：太短的 contig（<1000bp）信息太少，binning 工具分不准，直接扔掉
# --k-list：k-mer 就是"用多长的碎片去拼图"，多试几种长度，取最好的结果
# --memory 0.5：保险起见只用一半内存，避免被系统杀进程
# 检查组装结果
# 白话：看看拼出了多少条 contig、最长的多长、N50 是多少
cat ${OUTDIR}/log # 查看 MEGAHIT 运行日志
grep ">" ${OUTDIR}/final.contigs.fa | wc -l # 统计 contig 数量
# 输出文件：
# ${OUTDIR}/final.contigs.fa → 这就是组装出的 contigs，后面所有步骤都基于它
多样本联合组装 vs 单样本组装？
单样本组装（推荐初学者）：
优点：简单、快、结果容易理解
缺点：低丰度菌可能组装不出来
多样本联合组装（co-assembly）：
优点：数据量大，低丰度菌也能拼出来
缺点：不同样本的同源菌可能被拼成嵌合体（chimeric contigs）
你的 T2D 项目建议：
先单样本组装，熟悉流程后再尝试联合组装
4.2 比对回帖 —— 计算覆盖度¶
# ============================================================
# 第二步：把 clean reads 比对回 contigs，计算每条 contig 的覆盖度
# 目的：覆盖度是 binning 最重要的特征之一
# 白话：看每条 contig 被多少 reads 覆盖 → 推断这条 contig 来自哪个菌
# ============================================================
CONTIGS="assembly/${SAMPLE}/final.contigs.fa" # 组装结果
BAM_DIR="mapping/${SAMPLE}" # 比对结果目录
mkdir -p ${BAM_DIR}
# 第 2.1 步：建索引
# 白话：给 contigs 建目录（就像给一本书编页码，方便快速查找）
bwa-mem2 index ${CONTIGS} # 给 contigs 文件建索引
# 第 2.2 步：比对
# 白话：把每条 read 放回 contigs 上，看它原来属于哪条 contig
bwa-mem2 mem \
-t 16 \ # 16 线程并行
${CONTIGS} \ # 参考序列（contigs）
${CLEAN_R1} ${CLEAN_R2} \ # 要比对的 reads
2> ${BAM_DIR}/bwa.log \ # 日志输出
| samtools sort \ # 管道直接排序（按基因组位置排序）
-@ 8 \ # 排序用 8 线程
-o ${BAM_DIR}/${SAMPLE}.sorted.bam # 输出排序后的 BAM 文件
# 第 2.3 步：建 BAM 索引
# 白话：给排序后的比对结果也建一个目录，方便后续工具快速读取
samtools index ${BAM_DIR}/${SAMPLE}.sorted.bam
# 第 2.4 步：生成覆盖度文件（MetaBAT2 专用格式）
# 白话：统计每条 contig 平均被多少 reads 覆盖（覆盖深度）
jgi_summarize_bam_contig_depths \
--outputDepth ${BAM_DIR}/${SAMPLE}.depth.txt \ # 输出覆盖度文件
${BAM_DIR}/${SAMPLE}.sorted.bam # 输入 BAM 文件
# 查看覆盖度文件前几行
head -5 ${BAM_DIR}/${SAMPLE}.depth.txt
# 输出示例：
# contigName contigLen totalAvgDepth T2D_sample01.sorted.bam T2D_sample01.sorted.bam-var
# k141_1 3567 45.2 45.2 120.3
# k141_2 12890 12.8 12.8 35.6
多样本覆盖度（提升 binning 效果的关键）：
# 如果你有多个样本（比如 T2D 组 + 对照组），每个样本都比对回同一套 contigs
# 白话：同一个菌在不同样本里的丰度变化模式是一致的
# 多样本覆盖度让 binning 工具更准确地判断"谁和谁是一家的"
for SAMPLE in T2D_01 T2D_02 T2D_03 HC_01 HC_02 HC_03; do
bwa-mem2 mem -t 16 ${CONTIGS} \
clean_data/${SAMPLE}_R1.fq.gz clean_data/${SAMPLE}_R2.fq.gz \
| samtools sort -@ 8 -o mapping/${SAMPLE}.sorted.bam
samtools index mapping/${SAMPLE}.sorted.bam
done
# 一次性计算所有样本的覆盖度
jgi_summarize_bam_contig_depths \
--outputDepth all_samples.depth.txt \
mapping/*.sorted.bam # 所有 BAM 文件一起算
4.3 Binning（核心步骤）¶
# ============================================================
# 第三步：Binning —— 把 contigs 分到不同的"箱子"里
# 白话：根据覆盖度 + 序列特征，把可能来自同一个菌的 contig 分到一组
# 策略：用三个工具分别 binning，后面用 DAS Tool 取精华
# ============================================================
BIN_DIR="binning/${SAMPLE}"
mkdir -p ${BIN_DIR}/{metabat2,maxbin2,concoct}
# ==================
# 工具1：MetaBAT2（最常用，速度快）
# ==================
metabat2 \
-i ${CONTIGS} \ # 输入 contigs
-a ${BAM_DIR}/${SAMPLE}.depth.txt \ # 覆盖度文件（上一步生成的）
-o ${BIN_DIR}/metabat2/bin \ # 输出 bin 的前缀
-m 1500 \ # contig 最小长度阈值（1500bp）
-t 16 \ # 线程数
--minClsSize 200000 # bin 最小总大小（200kb，太小的 bin 没意义）
# 产出：metabat2/bin.1.fa, bin.2.fa, bin.3.fa, ...
# 每个 .fa 文件就是一个 bin，里面包含被归到同一组的 contigs
ls ${BIN_DIR}/metabat2/*.fa | wc -l # 看产出了多少个 bin
# ==================
# 工具2：MaxBin2（用 marker gene 辅助判断）
# ==================
# 先生成 MaxBin2 需要的覆盖度格式（它不认 MetaBAT2 的格式）
# 白话：MaxBin2 挑食，要给它专门格式的"菜单"
cut -f1,4 ${BAM_DIR}/${SAMPLE}.depth.txt | tail -n +2 \
> ${BIN_DIR}/maxbin2_abundance.txt # 只取 contig 名 + 覆盖度
run_MaxBin.pl \
-contig ${CONTIGS} \ # 输入 contigs
-abund ${BIN_DIR}/maxbin2_abundance.txt \ # 覆盖度文件
-out ${BIN_DIR}/maxbin2/bin \ # 输出前缀
-thread 16 # 线程数
# 产出：maxbin2/bin.001.fasta, bin.002.fasta, ...
# ==================
# 工具3：CONCOCT（对低丰度菌更敏感）
# ==================
# CONCOCT 需要先把 contigs 切成等长片段（10kb），再分别计算覆盖度
# 第 3.1 步：切割 contigs 为 10kb 片段
cut_up_fasta.py ${CONTIGS} \
-c 10000 \ # 切成 10kb 的片段
-o 0 \ # 无重叠
--merge_last \ # 最后一段不足 10kb 的合并到前一段
-b ${BIN_DIR}/concoct/contigs_10K.bed \ # 输出切割坐标文件
> ${BIN_DIR}/concoct/contigs_10K.fa # 输出切割后的 fasta
# 第 3.2 步：计算切割片段的覆盖度
concoct_coverage_table.py \
${BIN_DIR}/concoct/contigs_10K.bed \ # 切割坐标
${BAM_DIR}/${SAMPLE}.sorted.bam \ # BAM 文件
> ${BIN_DIR}/concoct/coverage_table.tsv # 覆盖度表
# 第 3.3 步：运行 CONCOCT
concoct \
--composition_file ${BIN_DIR}/concoct/contigs_10K.fa \ # 切割后的 contigs
--coverage_file ${BIN_DIR}/concoct/coverage_table.tsv \ # 覆盖度表
-b ${BIN_DIR}/concoct/ \ # 输出目录
-t 16 # 线程数
# 第 3.4 步：把切割片段的 bin 归属映射回原始 contigs
merge_cutup_clustering.py \
${BIN_DIR}/concoct/clustering_gt1000.csv \ # CONCOCT 的聚类结果
> ${BIN_DIR}/concoct/clustering_merged.csv # 合并后的结果
# 第 3.5 步：提取每个 bin 的 fasta
extract_fasta_bins.py \
${CONTIGS} \ # 原始 contigs
${BIN_DIR}/concoct/clustering_merged.csv \ # 聚类结果
--output_path ${BIN_DIR}/concoct/ # 输出目录
# 产出：concoct/1.fa, 2.fa, 3.fa, ...
4.4 Bin 精炼（DAS Tool）¶
# ============================================================
# 第四步：用 DAS Tool 整合三个 binner 的结果，取其精华
# 白话：三个老师各自改了一份试卷，现在找一个总监来取每个老师改得最好的部分
# 原理：DAS Tool 用单拷贝基因评分，从三套结果中挑出最优的 bin 组合
# ============================================================
REFINE_DIR="binning_refined/${SAMPLE}"
mkdir -p ${REFINE_DIR}
# 第 4.1 步：生成 DAS Tool 需要的输入格式
# 格式：contig名\tbin名（制表符分隔，每行告诉 DAS Tool 这条 contig 属于哪个 bin）
# MetaBAT2 的结果 → contigs2bin 格式
Fasta_to_Contig2Bin.sh \
-i ${BIN_DIR}/metabat2/ \ # MetaBAT2 输出目录
-e fa \ # bin 文件的后缀名
> ${REFINE_DIR}/metabat2_contigs2bin.tsv # 输出
# MaxBin2 的结果
Fasta_to_Contig2Bin.sh \
-i ${BIN_DIR}/maxbin2/ \
-e fasta \
> ${REFINE_DIR}/maxbin2_contigs2bin.tsv
# CONCOCT 的结果
Fasta_to_Contig2Bin.sh \
-i ${BIN_DIR}/concoct/ \
-e fa \
> ${REFINE_DIR}/concoct_contigs2bin.tsv
# 第 4.2 步：运行 DAS Tool
DAS_Tool \
-i ${REFINE_DIR}/metabat2_contigs2bin.tsv,\
${REFINE_DIR}/maxbin2_contigs2bin.tsv,\
${REFINE_DIR}/concoct_contigs2bin.tsv \ # 三个 binner 的结果（逗号分隔）
-l metabat2,maxbin2,concoct \ # 对应的标签名
-c ${CONTIGS} \ # 原始 contigs 文件
-o ${REFINE_DIR}/das_tool \ # 输出前缀
--write_bins \ # 输出精炼后的 bin fasta 文件
--write_bin_evals \ # 输出每个 bin 的评估分数
-t 16 \ # 线程数
--score_threshold 0.5 # 质量分阈值（低于 0.5 的 bin 丢弃）
# 产出：
# das_tool_DASTool_bins/ → 精炼后的 bin（这些才是你的"候选 MAGs"）
# das_tool_DASTool_summary.tsv → 每个 bin 来自哪个 binner、评估分数
ls ${REFINE_DIR}/das_tool_DASTool_bins/*.fa | wc -l # 看精炼后剩多少个 bin
# 查看总结报告
column -t ${REFINE_DIR}/das_tool_DASTool_summary.tsv | head -20
4.5 质量评估（CheckM2）¶
# ============================================================
# 第五步：用 CheckM2 评估每个 bin 的质量
# 目的：判断完整性（拼了多完整）和污染度（混了多少其他菌）
# 白话：给每个"拼图"打分——拼了百分之几？混进了几块别人的碎片？
# ============================================================
CHECKM_DIR="checkm2/${SAMPLE}"
mkdir -p ${CHECKM_DIR}
# 运行 CheckM2（比 CheckM1 快 10 倍以上，推荐）
checkm2 predict \
--input ${REFINE_DIR}/das_tool_DASTool_bins/ \ # DAS Tool 精炼后的 bin 目录
--output-directory ${CHECKM_DIR} \ # 输出目录
-x fa \ # bin 文件后缀
--threads 16 # 线程数
# 查看质量报告
column -t ${CHECKM_DIR}/quality_report.tsv | head -20
# 输出示例：
# Name Completeness Contamination Completeness_Model_Used
# bin.1 96.5 1.2 Neural Network
# bin.2 78.3 3.5 Neural Network
# bin.3 45.1 12.8 Neural Network ← 这个质量太差，不算 MAG
# ============================================================
# 筛选高质量和中质量 MAG（按 MIMAG 标准）
# 高质量：完整性 ≥ 90%，污染度 ≤ 5%
# 中质量：完整性 ≥ 50%，污染度 ≤ 10%
# ============================================================
# 用 awk 筛选（第 2 列 = 完整性，第 3 列 = 污染度）
awk -F'\t' 'NR==1 || ($2>=50 && $3<=10)' ${CHECKM_DIR}/quality_report.tsv \
> ${CHECKM_DIR}/medium_quality_MAGs.tsv # 中质量以上的 MAG 列表
awk -F'\t' 'NR==1 || ($2>=90 && $3<=5)' ${CHECKM_DIR}/quality_report.tsv \
> ${CHECKM_DIR}/high_quality_MAGs.tsv # 高质量 MAG 列表
echo "高质量 MAG 数量："
tail -n +2 ${CHECKM_DIR}/high_quality_MAGs.tsv | wc -l
echo "中质量 MAG 数量："
tail -n +2 ${CHECKM_DIR}/medium_quality_MAGs.tsv | wc -l
4.6 物种注释（GTDB-Tk）¶
# ============================================================
# 第六步：用 GTDB-Tk 对 MAG 进行物种分类
# 目的：搞清楚每个 MAG 到底是什么菌
# 白话：拿着拼好的"拼图"去菌种图鉴里比对，找到最像的已知菌
# 注意：GTDB-Tk 内存需求大（~210 GB），需要在服务器上跑
# ============================================================
GTDB_DIR="gtdbtk/${SAMPLE}"
mkdir -p ${GTDB_DIR}
# 只对中质量以上的 MAG 做注释（低质量的不浪费时间）
# 先把合格的 MAG 复制到一个目录
QUALIFIED_DIR="qualified_MAGs/${SAMPLE}"
mkdir -p ${QUALIFIED_DIR}
# 根据 CheckM2 结果筛选合格 MAG 并复制
while read -r name comp cont rest; do
[ "$name" == "Name" ] && continue # 跳过表头
cp ${REFINE_DIR}/das_tool_DASTool_bins/${name}.fa ${QUALIFIED_DIR}/
done < ${CHECKM_DIR}/medium_quality_MAGs.tsv
# 运行 GTDB-Tk 分类
gtdbtk classify_wf \
--genome_dir ${QUALIFIED_DIR} \ # MAG 文件目录
--out_dir ${GTDB_DIR} \ # 输出目录
-x fa \ # 文件后缀
--cpus 16 \ # CPU 数
--pplacer_cpus 4 # pplacer 用的 CPU（吃内存大户，控制数量）
# 查看分类结果
column -t ${GTDB_DIR}/classify/gtdbtk.bac120.summary.tsv | head -10
# 输出示例：
# user_genome classification fastani_reference
# bin.1 d__Bacteria;p__Firmicutes;c__Clostridia;...;s__Faecalibacterium prausnitzii
# bin.2 d__Bacteria;p__Bacteroidota;c__Bacteroidia;...;s__Bacteroides vulgatus
4.7 MAG 功能注释（Prokka / DRAM）¶
# ============================================================
# 第七步：功能注释 —— 搞清楚每个 MAG 里有什么基因、能干什么
# 白话：知道了"这是什么菌"之后，还要知道"这个菌有什么本事"
# ============================================================
ANNOT_DIR="annotation/${SAMPLE}"
mkdir -p ${ANNOT_DIR}
# ==================
# 方法1：Prokka（快速基因预测 + 基础注释）
# 适合：快速出结果、看基因数量和基本功能
# ==================
for MAG_FILE in ${QUALIFIED_DIR}/*.fa; do
MAG_NAME=$(basename ${MAG_FILE} .fa) # 提取文件名（不含后缀）
prokka \
${MAG_FILE} \ # 输入 MAG 文件
--outdir ${ANNOT_DIR}/prokka/${MAG_NAME} \ # 输出目录
--prefix ${MAG_NAME} \ # 输出文件前缀
--metagenome \ # 宏基因组模式（关闭 tRNA/rRNA 检测的严格限制）
--cpus 8 \ # CPU 数
--kingdom Bacteria \ # 物种域（细菌）
--force # 覆盖已有结果
done
# Prokka 输出的关键文件：
# *.gff → 基因坐标和注释（GFF3格式）
# *.faa → 预测的蛋白质序列
# *.ffn → 预测的核酸基因序列
# *.tsv → 注释统计表
# *.gbk → GenBank 格式文件（可导入 Artemis 查看）
# ==================
# 方法2：DRAM（深度代谢注释，更全面）
# 适合：分析代谢通路、碳源利用、抗性基因等
# ==================
# DRAM 需要先设置数据库（首次使用需要 prepare）
# DRAM-setup.py prepare_databases --output_dir /path/to/dram_db
# 注释
DRAM.py annotate \
-i "${QUALIFIED_DIR}/*.fa" \ # 输入 MAG 文件（支持通配符）
-o ${ANNOT_DIR}/dram/ \ # 输出目录
--threads 16 # 线程数
# 生成代谢汇总报告（DRAM 的杀手锏功能）
DRAM.py distill \
-i ${ANNOT_DIR}/dram/annotations.tsv \ # 上一步的注释结果
-o ${ANNOT_DIR}/dram/distill/ \ # 输出目录
--trna_path ${ANNOT_DIR}/dram/trnas.tsv \ # tRNA 结果
--rrna_path ${ANNOT_DIR}/dram/rrnas.tsv # rRNA 结果
# DRAM distill 输出：
# metabolism_summary.xlsx → 代谢通路热图（可直接用于论文图）
# genome_stats.tsv → 每个 MAG 的基因组统计
# product.html → 可视化交互报告
5. 工具对比：MetaBAT2 vs MaxBin2 vs CONCOCT vs SemiBin¶
| 维度 | MetaBAT2 | MaxBin2 | CONCOCT | SemiBin |
|---|---|---|---|---|
| 算法 | 密度聚类 | EM + marker gene | 高斯混合模型 | 深度学习（自监督对比学习） |
| 速度 | 最快 | 中等 | 最慢 | 中等（GPU 加速快） |
| 内存需求 | 低 | 中 | 高 | 中 |
| 对低丰度菌 | 一般 | 较好 | 较好 | 最好 |
| 多样本支持 | 支持（推荐） | 不直接支持 | 支持 | 支持 |
| 是否需要标记基因 | 否 | 是 | 否 | 否 |
| 使用难度 | 简单 | 简单 | 复杂（步骤多） | 简单 |
| 发表年份 | 2019 | 2016 | 2014 | 2022 |
| 推荐场景 | 日常首选 | 辅助验证 | 辅助验证 | 追求最优效果 |
建议策略：
初学者：MetaBAT2 一个工具跑完就行
进阶用法：MetaBAT2 + MaxBin2 + CONCOCT 三个都跑 → DAS Tool 精炼（本教程的方法）
最新方案：SemiBin2（2023年更新），单工具效果已接近多工具整合
SemiBin2 快速示例：
# SemiBin2 安装
mamba install -c bioconda semibin -y
# 一步完成 binning（内置了覆盖度计算）
SemiBin2 single_easy_bin \
-i ${CONTIGS} \ # contigs
--input-bam ${BAM_DIR}/${SAMPLE}.sorted.bam \ # BAM 文件
-o ${BIN_DIR}/semibin/ \ # 输出目录
--environment human_gut \ # 预训练模型（人肠道，适合你的T2D项目！）
-t 16 # 线程数
# --environment 可选值：human_gut, dog_gut, ocean, soil, cat_gut,
# human_oral, mouse_gut, pig_gut, built_environment, wastewater, global
# SemiBin2 的优势：有针对人肠道的预训练模型，不需要多个工具再整合
6. 质量标准：MIMAG（Minimum Information about MAGs）¶
MIMAG 是 2017 年 Nature Biotechnology 发表的 MAG 质量分级标准，现在是领域通用标准。
| 质量等级 | 完整性 | 污染度 | 额外要求 | 你能用它做什么 |
|---|---|---|---|---|
| 高质量 | ≥ 90% | ≤ 5% | 23S/16S/5S rRNA + ≥18 个 tRNA | 发论文、做精细功能分析 |
| 中质量 | ≥ 50% | ≤ 10% | 无 | 物种分类、初步功能注释 |
| 低质量 | < 50% | > 10% | 无 | 基本没用，丢弃 |
常见结果预期：
一般宏基因组测序（10-20 Gb/样本）：
- 总 bin 数：20-100 个
- 中质量以上 MAG：10-50 个
- 高质量 MAG：5-20 个
你的 T2D 项目如果有 6 个样本，多样本联合 binning 后：
- 预期能获得 30-80 个中质量以上 MAG
- 其中可能 10-30 个是高质量 MAG
7. 和你 T2D 项目的关联¶
你现有的数据和分析¶
已有数据：
├── T2D 组 + 健康对照组的宏基因组测序数据
├── 已完成：质控（fastp）→ 去宿主（BWA 去人基因组）→ Kraken2 物种分类
└── 已知结果：T2D 组和对照组的菌群组成差异
接下来可以做 Binning：
├── 直接用质控+去宿主后的 clean reads 作为输入
├── 从 MEGAHIT 组装开始，走完本教程的全流程
└── 最终得到每个样本的 MAGs + 功能注释
具体操作路径¶
# 你的 T2D 项目目录结构（参考）
cd /home/pweaz/t2d_ai_project/
# 假设 clean reads 在这个目录：
ls clean_data/
# T2D_01_R1.fq.gz T2D_01_R2.fq.gz
# T2D_02_R1.fq.gz T2D_02_R2.fq.gz
# HC_01_R1.fq.gz HC_01_R2.fq.gz
# ...
# 方案1：单样本独立组装 + binning（简单，推荐先试）
for SAMPLE in T2D_01 T2D_02 HC_01 HC_02; do
# 1. 组装
megahit -1 clean_data/${SAMPLE}_R1.fq.gz \
-2 clean_data/${SAMPLE}_R2.fq.gz \
-o assembly/${SAMPLE} --min-contig-len 1000 -t 16
# 2-7. 后续步骤按本教程执行
done
# 方案2：多样本联合组装（效果更好，但更复杂）
# 把所有样本的 reads 合并后一起组装
cat clean_data/*_R1.fq.gz > all_R1.fq.gz
cat clean_data/*_R2.fq.gz > all_R2.fq.gz
megahit -1 all_R1.fq.gz -2 all_R2.fq.gz -o assembly/coassembly \
--min-contig-len 1000 -t 32
# 然后每个样本分别比对回 coassembly 的 contigs，获得多样本覆盖度
MAGs 能给你的论文/面试加什么分¶
当前分析：
"T2D 患者肠道中 Faecalibacterium 丰度显著降低"
加上 MAGs 后：
"我们提取了 F. prausnitzii 的高质量 MAG（完整性 95.3%，污染度 1.1%），
功能注释发现其携带完整的丁酸合成通路（but-bcd-etfAB），
而 T2D 组该通路关键基因 bcd 的覆盖度显著低于健康组（p=0.003），
提示 T2D 患者肠道中的 F. prausnitzii 丁酸产生能力可能减弱。"
面试加分点：
- 说明你不只停留在"物种注释"层面，还能做深入的功能分析
- 展示你掌握了宏基因组下游分析的核心技能
- 证明你理解从数据到生物学意义的完整闭环
8. 常见报错与解决¶
报错1：MEGAHIT 内存不足被 kill¶
错误信息：
Signal 9 (SIGKILL) - 进程被系统强制终止（通常是 OOM Killer）
原因：宏基因组组装非常吃内存，默认用全部内存可能被系统 kill
解决方案：
# 方法1：限制内存使用
megahit ... --memory 0.5 # 只用 50% 内存
# 方法2：用更大的服务器或申请更多资源
salloc -n 1 --mem=200G # SLURM 集群申请 200GB 内存
# 方法3：利用 MEGAHIT 的断点续跑
megahit ... --continue # 从上次中断的地方继续
报错2：MetaBAT2 输出 0 个 bin¶
错误信息：
没有报错，但 bin 目录是空的
常见原因：
1. contig 太少或太短 → 降低 --min-contig-len（但不建议低于 1000）
2. 覆盖度文件格式不对 → 检查 depth.txt 是否有内容
3. 样本测序深度不够 → 数据量太少，组装出的 contig 本身就少
排查步骤：
# 检查 contigs 数量和长度分布
grep -c ">" final.contigs.fa
# 如果只有几百条 contig，binning 效果本来就差
# 检查覆盖度文件
head -5 depth.txt
# 确认第 4 列（覆盖度）不是全 0
# 降低阈值再试
metabat2 -i contigs.fa -a depth.txt -o bin -m 1000 --minClsSize 50000
报错3：CheckM2 数据库找不到¶
错误信息：
FileNotFoundError: CheckM2 database not found
解决方案：
# 下载 CheckM2 数据库
checkm2 database --download --path /path/to/checkm2_db
# 设置环境变量
export CHECKM2DB=/path/to/checkm2_db/uniref100.KO.1.dmnd
# 或在命令中指定
checkm2 predict ... --database_path /path/to/checkm2_db
报错4：GTDB-Tk pplacer 内存溢出¶
错误信息：
pplacer: Fatal error: out of memory
原因：pplacer 是 GTDB-Tk 里最吃内存的步骤（需要 ~210 GB）
解决方案：
# 方法1：减少 pplacer 线程（每个线程都独立占内存！）
gtdbtk classify_wf ... --pplacer_cpus 1 # 只用 1 个线程
# 方法2：用 GTDB-Tk 的 ani_screen 跳过 pplacer（v2.3+）
gtdbtk classify_wf ... --skip_ani_screen False
# 方法3：改用 GTDB-Tk v2 + --mash_db 预筛选，减少 pplacer 负担
报错5：DAS Tool "No bins passed the quality threshold"¶
错误信息：
WARNING: No bins passed the quality threshold
原因：所有 bin 的质量分都低于默认阈值（0.5）
解决方案：
# 方法1：降低阈值（但质量也会下降）
DAS_Tool ... --score_threshold 0.3
# 方法2：检查各 binner 的结果是否正常
ls binning/metabat2/*.fa | wc -l # MetaBAT2 产出了多少 bin
ls binning/maxbin2/*.fasta | wc -l # MaxBin2 产出了多少 bin
# 方法3：数据本身测序深度不够 → 考虑多样本联合 binning 提升效果
报错6：CONCOCT 报错 "ValueError: could not convert string to float"¶
错误信息：
ValueError: could not convert string to float: 'contigName'
原因：覆盖度文件包含了表头行，CONCOCT 不认
解决方案：
# 检查覆盖度文件的第一行
head -1 coverage_table.tsv
# 如果有文字表头，用 tail 跳过
tail -n +2 coverage_table.tsv > coverage_table_noheader.tsv
# 或者确认用的是 concoct_coverage_table.py 生成的文件（不会有此问题）
# 手动从 jgi_summarize_bam_contig_depths 转格式时容易出这个错
报错7：Prokka 在宏基因组模式下漏掉基因¶
错误信息：
不是报错，但注释的基因数明显少于预期
原因：Prokka 默认会跳过"太短的 contig"上的部分基因预测
解决方案：
# 加 --metagenome 标记（放宽预测限制）
prokka ... --metagenome --compliant
# 或者改用 Prodigal 做基因预测后，再用 eggNOG-mapper 注释
prodigal -i mag.fa -a mag.faa -p meta # -p meta 是宏基因组模式
9. 速查表¶
MAG 提取流程速查¶
| 步骤 | 工具 | 核心命令 | 输入 | 输出 |
|---|---|---|---|---|
| 组装 | MEGAHIT | megahit -1 R1 -2 R2 -o out |
clean reads | contigs.fa |
| 建索引 | bwa-mem2 | bwa-mem2 index contigs.fa |
contigs.fa | 索引文件 |
| 比对 | bwa-mem2 + samtools | bwa-mem2 mem ref R1 R2 \| samtools sort -o out.bam |
reads + contigs | sorted.bam |
| 覆盖度 | jgi_summarize | jgi_summarize_bam_contig_depths --outputDepth depth.txt *.bam |
BAM 文件 | depth.txt |
| Binning | MetaBAT2 | metabat2 -i contigs.fa -a depth.txt -o bin |
contigs + depth | bin.*.fa |
| 精炼 | DAS Tool | DAS_Tool -i c2b1,c2b2 -c contigs.fa -o out |
多个 binner 结果 | 精选 bin |
| 质控 | CheckM2 | checkm2 predict --input bins/ -o out |
bin 目录 | quality_report.tsv |
| 物种注释 | GTDB-Tk | gtdbtk classify_wf --genome_dir bins/ -o out |
MAG 目录 | 分类结果 |
| 功能注释 | Prokka | prokka mag.fa --outdir out --metagenome |
MAG 文件 | GFF/FAA/GBK |
| 代谢注释 | DRAM | DRAM.py annotate -i 'mags/*.fa' -o out |
MAG 文件 | 代谢通路报告 |
MIMAG 质量标准速查¶
| 等级 | 完整性 | 污染度 | rRNA | tRNA |
|---|---|---|---|---|
| 高质量 | ≥ 90% | ≤ 5% | 23S+16S+5S | ≥ 18 |
| 中质量 | ≥ 50% | ≤ 10% | - | - |
| 低质量 | < 50% | > 10% | - | - |
关键参数速查¶
| 参数 | 推荐值 | 为什么 |
|---|---|---|
| contig 最小长度 | 1000-1500 bp | 太短的 contig 特征信息不够，分不准 |
| MEGAHIT --memory | 0.5 | 防止内存溢出被 kill |
| MetaBAT2 -m | 1500 | 和 contig 最小长度一致 |
| DAS Tool --score_threshold | 0.5 | 默认值，低质量 bin 会被过滤 |
| GTDB-Tk --pplacer_cpus | 1-4 | 每个线程独立占内存，开太多会炸 |
| CheckM2 完整性阈值 | ≥ 50% | MIMAG 中质量标准 |
| CheckM2 污染度阈值 | ≤ 10% | MIMAG 中质量标准 |
资源需求速查¶
| 步骤 | CPU | 内存 | 磁盘 | 时间（6样本） |
|---|---|---|---|---|
| MEGAHIT | 16-32 | 50-200 GB | 输入的 2-5 倍 | 2-8 小时/样本 |
| BWA 比对 | 16 | 10-20 GB | BAM 约 = 输入大小 | 1-2 小时/样本 |
| MetaBAT2 | 16 | 5-10 GB | 很小 | 10-30 分钟 |
| DAS Tool | 16 | 10-20 GB | 很小 | 30-60 分钟 |
| CheckM2 | 16 | 10-30 GB | ~1 GB | 10-30 分钟 |
| GTDB-Tk | 16 | 210+ GB | ~100 GB（数据库） | 2-6 小时 |
| Prokka | 8 | 2-4 GB | 很小 | 5-15 分钟/MAG |
10. 延伸学习资源¶
必读论文¶
| 论文 | 要点 | 推荐度 |
|---|---|---|
| Bowers et al., 2017 (Nature Biotechnology) | MIMAG 标准的定义论文，MAG 质量分级的出处 | 必读 |
| Kang et al., 2019 (PeerJ) | MetaBAT2 原始论文 | 推荐 |
| Sieber et al., 2018 (Nature Microbiology) | DAS Tool 原始论文 | 推荐 |
| Parks et al., 2015 (Genome Research) | CheckM 原始论文 | 推荐 |
| Chklovski et al., 2023 (Nature Methods) | CheckM2 论文（机器学习改进版） | 推荐 |
| Pan et al., 2023 (Nature Communications) | SemiBin2 论文（深度学习 binning） | 了解前沿 |
实操教程¶
- Happy Belly Bioinformatics (astrobiomike.github.io)：最友好的宏基因组入门教程，有完整的 binning 实操
- Metagenomics Wiki (metagenomics-wiki.readthedocs.io)：涵盖从测序到 MAG 分析的完整知识体系
- ATLAS pipeline (metagenome-atlas.github.io)：一键式宏基因组分析流程，内置 binning + 注释
进阶方向¶
| 方向 | 工具/方法 | 适用场景 |
|---|---|---|
| 长读长 binning | SemiBin2 + HiFi reads | PacBio/Nanopore 数据，更完整的 MAG |
| 菌株级分析 | inStrain / StrainPhlAn | 从 MAG 中分辨菌株差异 |
| 代谢网络重建 | CarveMe + SMETANA | 从 MAG 构建代谢模型，预测菌群互作 |
| 病毒 MAG | vRhyme / VIBRANT | 从宏基因组中提取病毒基因组 |
| 泛基因组分析 | Roary / PPanGGOLiN | 比较同物种不同 MAG 的基因组内容 |
| 宏基因组流程管理 | nf-core/mag, ATLAS | 用 Nextflow/Snakemake 自动化整个流程 |
与你已有知识的衔接¶
你已经会的： 接下来可以学的：
Kraken2 物种分类 → MAGs 功能注释（本教程）
fastp 质控 → 已经是 binning 的前置步骤
BWA 去宿主 → BWA 比对回帖（同一个工具，不同用法）
随机森林分类 → 用 MAG 功能特征做机器学习（进阶）
Snakemake 流程管理 → 把整个 binning 流程写成 Snakefile（强烈推荐！）
最后更新：2026-05 | 审核状态：自审通过 | 预估字数：~4500 字