宏基因组分析全流程

从原始测序数据到生物学发现 — 每一步的输入、输出和结果

一图看懂全流程

1. 质控 QC
→
2. 去宿主
→
3. 物种注释
→
4. 功能注释
→
5. 组装分箱
→
6. 统计分析

原始 FASTQ → 干净 reads → 物种丰度表 → 功能通路表 → MAGs → 差异分析/可视化

01 质量控制 (Quality Control) 必做 ▶

去掉低质量碱基、接头序列和过短的 reads，保证后续分析准确。

输入文件

sample1_R1.fastq.gz（正向reads，~5GB）
sample1_R2.fastq.gz（反向reads，~5GB）

测序公司给的原始数据，通常双端 150bp PE reads

→

输出文件

sample1_clean_R1.fastq.gz（质控后正向）
sample1_clean_R2.fastq.gz（质控后反向）
fastp.html（质控报告）
fastp.json（质控统计JSON）
fastp（推荐，快且全） Trimmomatic FastQC（只看不切）
# fastp 一行搞定质控（去接头+去低质量+去短reads）
fastp \
  -i sample1_R1.fastq.gz \
  -I sample1_R2.fastq.gz \
  -o sample1_clean_R1.fastq.gz \
  -O sample1_clean_R2.fastq.gz \
  -q 20 -l 50 \
  --detect_adapter_for_pe \
  -h fastp.html -j fastp.json
📖 参数逐行讲解
-i / -I 输入文件：小写 i = 正向 R1，大写 I = 反向 R2（双端测序必须成对）
-o / -O 输出文件：对应的质控后 R1/R2 文件
-q 20 质量阈值：碱基质量低于 Q20（错误率 1%）会被修剪掉
-l 50 最短长度：修剪后不足 50bp 的 reads 直接丢弃
--detect_adapter_for_pe 自动检测双端接头序列并去除（不需要手动指定接头文件）
-h / -j 输出质控报告：HTML 可视化报告 + JSON 机器可读报告
fastp JSON 输出（真实格式）：
"before_filtering": {
  "total_reads": 6510970, ← R1+R2 总 reads 数
  "total_bases": 976645500, ← 总碱基数 (bp)
  "q20_rate": 0.942172, ← Q20 占比 (错误率≤1%的碱基)
  "q30_rate": 0.873813, ← Q30 占比 (错误率≤0.1%的碱基)
  "gc_content": 0.440688 ← GC含量 (正常40-60%)
},
"filtering_result": {
  "passed_filter_reads": 6312238, ← 通过质控的 reads (97%)
  "low_quality_reads": 172442, ← 低质量被丢弃 (>40%碱基低于Q15)
  "too_many_N_reads": 2, ← N碱基太多被丢弃 (N>5个)
  "too_short_reads": 24358, ← 太短被丢弃 (默认<15bp)
  "low_complexity_reads": 1930 ← 低复杂度被丢弃 (重复序列)
},
"adapter_cutting": {
  "adapter_trimmed_reads": 124850, ← 切掉接头的 reads 数
  "adapter_trimmed_bases": 1562340 ← 切掉的接头碱基数
},
"duplication": {
  "rate": 0.12453 ← PCR重复率 ~12.5%
}
质控前后 reads 数量对比
原始 reads
48,523,678
质控后 reads
45,891,234 (94.6%)
去掉的接头
1.2%
去掉的低质量
3.1%
去掉的太短
1.1%
碱基质量分布（Q 值标准）
指标质控前质控后合格标准判定
Q20 (%)95.1%97.3%> 95%✅ 合格
Q30 (%)88.7%93.1%> 85%✅ 合格
GC 含量48.5%48.5%40-60%✅ 正常
reads 通过率-94.6%> 80%✅ 良好
  • GC 含量异常：偏离 40-60% 说明可能有污染（接头残留、其他物种混入）
  • 通过率过低 (<70%)：测序质量差，需联系测序公司重测或加大测序量
  • 不要过度质控：Q 值阈值不要设太高（>25），否则会丢太多数据
  • 双端 reads 必须配对：质控后 R1 和 R2 数量必须一致，fastp 自动处理
  • 多样本批量处理：用 for 循环或 Snakemake，不要一个个手动跑
Q20 = 错误率 1%，Q30 = 错误率 0.1%。宏基因组数据 Q30 > 85% 就算合格。fastp 比 Trimmomatic 快 2-5 倍且自动检测接头。
02 去宿主污染 (Host Removal) 人/动物样本必做 ▶

把人的 DNA 序列去掉，只留微生物的序列。人肠道样本宿主 DNA 可达 90%+！

输入文件

sample1_clean_R1.fastq.gz
sample1_clean_R2.fastq.gz
hg38 人类参考基因组索引
→

输出文件

sample1_nohost_R1.fastq.gz（微生物reads）
sample1_nohost_R2.fastq.gz（微生物reads）
host_mapping.log（比对统计）
Bowtie2（最常用） BWA KneadData（整合质控+去宿主）
# Step1: Bowtie2 比对到人类基因组
bowtie2 -x /db/hg38/hg38 \
  -1 sample1_clean_R1.fastq.gz \
  -2 sample1_clean_R2.fastq.gz \
  --very-sensitive -p 8 \
  -S sample1_mapped.sam

# Step2: 提取没比对上的（=微生物的）reads
samtools view -bS -f 12 sample1_mapped.sam | \
samtools sort -n | \
samtools fastq -1 sample1_nohost_R1.fastq.gz \
  -2 sample1_nohost_R2.fastq.gz -
📖 参数逐行讲解
-x /db/hg38/hg38 宿主参考基因组索引路径（hg38 = 人类基因组第38版）
-1 / -2 双端输入文件：R1 和 R2
--very-sensitive 最高灵敏度模式：尽量把所有宿主 reads 都找出来（速度慢但去得干净）
-p 8 使用 8 个线程并行
-S 输出 SAM 格式比对文件
samtools view -bS -f 12 -bS SAM转BAM；-f 12 只保留双端都没比对上的 reads（flag 12 = 4+8，即 R1 和 R2 都 unmapped = 微生物 reads）
samtools sort -n 按 read name 排序（配对提取需要同名排列）
samtools fastq -1 ... -2 ... 从 BAM 提取回 FASTQ 格式（分别输出 R1/R2）
Bowtie2 比对输出（真实 stderr 格式）：
45891234 reads; of these: ← 输入的 reads 总数（R1+R2 每条各算一条）
  45891234 (100.00%) were paired; of these: ← 双端配对的 reads
    7633343 (16.63%) aligned concordantly 0 times ← 双端都没比上 = 微生物 reads ✅
    35128456 (76.55%) aligned concordantly exactly 1 time ← 双端唯一比对 = 宿主 DNA
    3129435 (6.82%) aligned concordantly >1 times ← 多处比对 = 宿主重复区域
  ----
  7633343 pairs aligned concordantly 0 times; of these:
    234567 (3.07%) aligned discordantly 1 time ← 不协调比对（方向/距离异常）
  ----
  7398776 pairs aligned 0 times concordantly or discordantly; of these:
    14797552 mates make up the pairs; of these:
      13245678 (89.51%) aligned 0 times ← 完全没比上的单端 reads
      1234567 (8.34%) aligned exactly 1 time
      317307 (2.14%) aligned >1 times
83.37% overall alignment rate ← 总比对率 = 宿主占比（越高宿主越多）
─────────────────────────────
关键：用 samtools -f 12 提取"concordantly 0 times"的 reads = 微生物
-f 12 = FLAG 位第3位(mate未比对)+第4位(自身未比对) = 双端都没比上
不同样本类型宿主 DNA 占比
人肠道粪便
~83% 宿主
微生物 ~17%
唾液
~92% 宿主
微生物 ~8%
皮肤
~95% 宿主
微生物 ~5%
土壤
~0%
微生物 ~100%
海水
~0%
微生物 ~100%
  • 环境样本可跳过此步：土壤、水体等无宿主样本直接进入物种注释
  • 参考基因组要匹配：人用 hg38，小鼠用 mm10，猪用 susScr11，用错了白跑
  • 微生物 reads 太少要警惕：去宿主后 <100 万条 reads，后续分析可能不可靠
  • SAM 文件巨大：中间产生的 .sam 文件可达几十 GB，分析完记得删除
  • -f 12 参数含义：双端都没比对上 = 非宿主的微生物序列（面试常问！）
-f 12 表示"双端都没比对上"，即非宿主的微生物序列。环境样本（土壤、水）不需要去宿主。临床样本宿主 DNA 占比高，这一步直接决定后续数据量。
03 物种注释 (Taxonomic Profiling) 核心步骤 ▶

回答"样本里有哪些微生物？各占多少？" — 生成物种丰度表，是宏基因组最核心的结果之一。

输入文件

sample1_nohost_R1.fastq.gz
sample1_nohost_R2.fastq.gz
Kraken2/MetaPhlAn 数据库
→

输出文件（Kraken2）

sample1.kraken（每条read的分类）
sample1_report.txt（物种丰度报告）

输出文件（MetaPhlAn）

sample1_profile.txt（物种丰度表）
sample1.bowtie2.bz2（比对中间文件）
Kraken2（k-mer，快） MetaPhlAn4（marker gene，准） Bracken（Kraken2 后续丰度校正）
# 方法A: Kraken2（速度快，适合大数据）
kraken2 --db /db/kraken2_std \
  --paired sample1_nohost_R1.fastq.gz sample1_nohost_R2.fastq.gz \
  --output sample1.kraken \
  --report sample1_report.txt \
  --threads 8

# 方法B: MetaPhlAn4（精度高，更常发文章）
metaphlan sample1_nohost_R1.fastq.gz,sample1_nohost_R2.fastq.gz \
  --input_type fastq \
  --bowtie2db /db/metaphlan \
  -o sample1_profile.txt \
  --nproc 8
📖 参数逐行讲解
▸ Kraken2（k-mer 精确匹配，速度快）
--db 数据库路径：标准库含细菌/古菌/病毒/人类（~70GB）
--paired 声明双端输入（两个文件用空格隔开）
--output 每条 read 的分类结果（逐行）
--report 汇总报告：每个物种的 reads 数和占比（分析主要看这个）
▸ MetaPhlAn4（标记基因比对，精度高）
--input_type fastq 输入格式声明（也支持 bowtie2out、sam）
--bowtie2db MetaPhlAn 专用标记基因数据库路径（~15GB，含 clade-specific markers）
-o 输出物种丰度表（相对丰度 %，各分类层级）
--nproc 线程数
Kraken2 report（真实 6 列格式，来自官方手册）：
%reads clade_reads taxon_reads rank taxID name
↑ ↑ ↑ ↑ ↑ ↑
占比% 该分支下 仅本节点 分类等级 NCBI 物种名
所有reads 直接匹配 D/P/C/ 税号
的reads O/F/G/S
─────────────────────────────────────────────
25.31 1932045 0 P 1239 Firmicutes ← 门级：0表示没有reads直接匹配到门，都在下级
8.23 628451 312056 G 816 Bacteroides ← 属级：312056条直接匹配到属，其余在种级
5.67 432891 198234 S 821 B. vulgatus ← 种级：198234条精确匹配到这个种
4.12 314567 156789 S 33039 F. prausnitzii ← 产丁酸关键菌！
3.89 296987 142345 S 1680 B. longum ← 益生菌代表
─────────────────────────────────────────────
⚠ rank 列：D=域 P=门 C=纲 O=目 F=科 G=属 S=种 S1=株
⚠ clade_reads ≥ taxon_reads（因为包含所有子节点）
⚠ 必须跑 Bracken 才能得到准确的种级丰度估计！
MetaPhlAn4 profile（真实格式，来自官方 wiki）：
#clade_name NCBI_tax_id relative_abundance additional_species
↑ 管道符分隔的分类层级 ↑ ↑
k=界 p=门 c=纲 o=目 f=科 g=属 s=种 NCBI税号 相对丰度(%)
─────────────────────────────────────────────
k__Bacteria 2 100.0 ← 界级加起来=100%
k__Bacteria|p__Firmicutes 1239 42.35 ← 门级也各自加起来=100%
k__Bacteria|p__Bacteroidetes 976 31.28
k__Bacteria|p__Proteobacteria 1224 12.67
k__Bacteria|...|g__Bacteroides|s__B_vulgatus 821 8.92 ← 种级丰度，所有种也加起来=100%
k__Bacteria|...|s__F_prausnitzii 33039 6.45
k__Bacteria|...|t__F_prausnitzii_SGB15044 33039 6.45 ← t__前缀=strain级(MetaPhlAn4新增)
─────────────────────────────────────────────
⚠ 每个分类层级独立加到100%（门级总和=100%，种级总和=100%）
⚠ MetaPhlAn4 新增 t__ 前缀表示 strain 级别的亚种分辨率
⚠ 只能检测数据库中有 marker gene 的物种，未知物种会被漏掉
门水平物种丰度（典型肠道样本）
Firmicutes
42.3%
Bacteroidetes
31.3%
Proteobacteria
12.7%
Actinobacteria
7.8%
Verrucomicrobia
3.5%
Others
2.4%
种水平 TOP10 物种丰度
排名物种丰度功能特征
1Bacteroides vulgatus8.9%多糖降解
2Faecalibacterium prausnitzii6.5%产丁酸，抗炎
3Prevotella copri4.2%纤维发酵
4Bacteroides uniformis3.8%多糖降解
5Roseburia intestinalis3.5%产丁酸
6Bifidobacterium longum3.1%益生菌
7Eubacterium rectale2.8%产丁酸
8Akkermansia muciniphila2.5%黏蛋白降解，改善代谢
9Escherichia coli1.9%条件致病菌
10Alistipes putredinis1.7%蛋白质发酵
  • 数据库版本要一致：同一项目所有样本必须用同版本数据库，否则结果不可比
  • Kraken2 标准库需 64GB 内存：内存不够用 minikraken2 (8GB)，但精度会下降
  • MetaPhlAn 只能检测数据库中有的物种：未知物种会被漏掉
  • Bracken 校正很重要：Kraken2 原始结果的丰度不准，必须跑 Bracken 再校正
  • unclassified 比例：Kraken2 通常 20-40% reads 无法分类，这是正常的
  • 物种名要规范：结果里用拉丁名（斜体），面试时能说出常见菌的中文名加分
Kraken2 用 k-mer 精确匹配（快但假阳性多），MetaPhlAn 用 marker gene 比对（慢但准）。发文章多用 MetaPhlAn，大规模筛选用 Kraken2。两者可互补验证。Bracken 可对 Kraken2 结果做丰度再校正。
04 功能注释 (Functional Annotation) 核心步骤 ▶

回答"微生物能做什么？" — 把基因/reads 映射到功能数据库（KEGG、GO、COG），看代谢通路有没有差异。

输入文件

sample1_nohost_R1.fastq.gz（基于reads）
或 genes.faa（基于组装后的基因）
MetaPhlAn 物种丰度（HUMAnN需要）
→

输出文件（HUMAnN）

sample1_genefamilies.tsv（基因家族丰度）
sample1_pathabundance.tsv（通路丰度）
sample1_pathcoverage.tsv（通路覆盖度）

输出文件（eggNOG-mapper）

sample1.emapper.annotations（COG/GO/KEGG）
HUMAnN3（通路，基于reads） eggNOG-mapper（COG/GO/KEGG，基于基因） DIAMOND（快速蛋白比对） 专用库: CAZy, VFDB, CARD
# HUMAnN3: reads → 通路丰度（自动调用MetaPhlAn）
humann --input sample1_nohost_cat.fastq.gz \
  --output humann_out/ \
  --threads 8

# eggNOG-mapper: 预测基因 → COG/GO/KEGG 注释
emapper.py -i genes.faa \
  --output sample1 \
  -m diamond \
  --cpu 8
📖 参数逐行讲解
▸ HUMAnN3（reads → 代谢通路丰度）
--input 输入：合并的 R1+R2 FASTQ（需要先 cat 合并双端）
--output 输出目录：包含 genefamilies.tsv、pathabundance.tsv、pathcoverage.tsv 三个核心文件
--threads 内部自动调用 MetaPhlAn → Diamond → MinPath 三步流程
▸ eggNOG-mapper（基因 → 功能注释）
-i genes.faa 输入：预测的蛋白序列文件（来自 Prodigal/MetaGeneMark）
-m diamond 比对方法：Diamond 蛋白比对（比 BLAST 快 500x）
--output 输出前缀：生成 .annotations 文件含 COG/GO/KEGG 注释
--cpu 线程数
HUMAnN3 通路丰度表（真实 pathabundance.tsv 格式）：
# Pathway RPKs
↑ MetaCyc 通路ID: 通路名称 ↑ Reads Per Kilobase（标准化丰度）
─────────────────────────────────────────────
UNMAPPED 245891.2 ← 没比对到任何基因的reads（正常30-50%）
UNINTEGRATED 89234.5 ← 比对到基因但无法归入任何通路
PWY-5100: pyruvate fermentation to acetate 3456.7 ← 通路总丰度（所有菌的贡献之和）
PWY-5100|g__Bacteroides.s__B_vulgatus 2134.5 ← 分层！这个菌贡献了 62% 的该通路
PWY-5100|g__Faecalibacterium.s__F_prausnitzii 891.2 ← 另一个菌贡献了 26%
PWY-5100|unclassified 431.0 ← 无法确定是哪个菌贡献的部分
GLYCOLYSIS: glycolysis I (from glucose 6-P) 2987.3
PWY-6609: adenine/adenosine nucleotide salvage 1234.5
─────────────────────────────────────────────
⚠ 分层格式 = 通路|g__属名.s__种名，HUMAnN 独有特色！
⚠ RPKs 已经按基因长度标准化，可跨样本比较
⚠ 跨样本比较前需要再做 CPM 标准化（humann_renorm_table）
eggNOG-mapper v2 注释结果（真实 .annotations 格式）：
#query seed_ortholog evalue score ... COG_cat Description GOs KEGG_ko KEGG_Pathway
↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑
基因ID 最佳同源 E-value 比对 COG 功能描述 GO术语 KEGG编号 KEGG通路
(越小越好) 得分 分类
─────────────────────────────────────────────
gene_001 1038710.MER 1.2e-45 234 .. G Hexokinase GO:0005975,... K00844 ko00010 ← G=碳水代谢 K00844=己糖激酶
gene_002 562.b1234 3.4e-78 456 .. E Glutamine synthetase GO:0006520,... K01915 ko00250 ← E=氨基酸代谢
gene_003 820.c5678 7.8e-23 123 .. C Succinate dehydrogenase GO:0055114,... K00239 ko00020 ← C=能量代谢 (TCA循环)
gene_004 1234.d910 2.1e-56 345 .. K Transcription factor GO:0006351,... K02529 ko02020 ← K=转录调控
gene_005 5678.e111 9.9e-12 089 .. S Function unknown - - - ← S=功能未知（约20%基因）
─────────────────────────────────────────────
⚠ COG单字母分类：G碳水 E氨基酸 C能量 K转录 L复制修复 S未知
⚠ eggNOG-mapper 输入必须是蛋白序列(.faa)，需先 Prodigal 预测基因
⚠ 约 20-30% 基因注释为"功能未知"，这是微生物组的常态
KEGG 通路丰度 TOP8（HUMAnN3 输出示例）
糖酵解
2987 RPKs
丙酮酸发酵
3457 RPKs
支链氨基酸合成
988 RPKs
腺嘌呤补救通路
1235 RPKs
TCA 循环
1923 RPKs
丁酸合成
1487 RPKs
叶酸合成
765 RPKs
脂多糖合成
634 RPKs
COG 功能分类统计（eggNOG-mapper 输出示例）
COG 类别含义基因数占比
G碳水化合物代谢3,24515.2%
E氨基酸代谢2,89113.5%
C能量代谢2,15610.1%
K转录调控1,9879.3%
LDNA 复制与修复1,6547.7%
S功能未知4,12319.3%
  • HUMAnN3 需先跑 MetaPhlAn：它依赖物种丰度来选择参考库，不能跳过
  • UNMAPPED 和 UNINTEGRATED 很正常：通常 30-60% 的 reads 无法映射到已知通路
  • 双端 reads 要合并：HUMAnN 输入需要 cat R1.fq.gz R2.fq.gz > combined.fq.gz
  • eggNOG-mapper 需要预测基因：必须先组装 + Prodigal 预测基因，再输入 .faa 文件
  • RPKs 单位：Reads Per Kilobase，已标准化的丰度，可跨样本比较
  • 专用数据库按需使用：CAZy（碳水酶）、CARD（抗性基因）、VFDB（毒力因子）不是每个项目都要跑
HUMAnN 特色是"物种贡献分层"——能看到每条通路由哪些菌贡献，适合研究菌群功能。eggNOG-mapper 适合组装后的基因注释。做碳水化合物酶用 CAZy，抗性基因用 CARD/ResFinder，毒力因子用 VFDB。
05 组装与分箱 (Assembly & Binning) 进阶分析 ▶

把短 reads 拼成长 contigs，再按"哪些 contig 来自同一个菌"分组成 MAG（Metagenome-Assembled Genome）。相当于从碎纸条还原出一本书。

输入文件

sample1_nohost_R1.fastq.gz
sample1_nohost_R2.fastq.gz
→

组装输出

contigs.fa（拼接好的连续序列）
assembly_stats.txt（组装统计）

分箱输出

bin.1.fa, bin.2.fa ...（每个MAG一个文件）

质控输出

checkm_results.tsv（MAG质量评估）
MEGAHIT（组装，省内存） metaSPAdes（组装，更准） MetaBAT2（分箱） MaxBin2（分箱） CONCOCT（分箱） DAS Tool（合并多个分箱结果） CheckM（MAG质量评估）
# Step1: MEGAHIT 组装
megahit -1 sample1_nohost_R1.fastq.gz \
  -2 sample1_nohost_R2.fastq.gz \
  -o megahit_out -t 8

# Step2: Bowtie2 把 reads 比对回 contigs（算覆盖度）
bowtie2-build megahit_out/final.contigs.fa contigs_index
bowtie2 -x contigs_index -1 *_R1.fq.gz -2 *_R2.fq.gz -S mapped.sam -p 8
samtools sort mapped.sam -o mapped.sorted.bam

# Step3: MetaBAT2 分箱
jgi_summarize_bam_contig_depths --outputDepth depth.txt mapped.sorted.bam
metabat2 -i megahit_out/final.contigs.fa -a depth.txt -o bins/bin

# Step4: CheckM 质量评估
checkm lineage_wf bins/ checkm_out -t 8 -x fa
📖 参数逐行讲解
▸ Step1: MEGAHIT 组装
-1 / -2 双端输入（去宿主后的 clean reads）
-o 输出目录：核心文件 final.contigs.fa（组装好的连续序列）
-t 8 线程数（MEGAHIT 内存友好，适合服务器资源有限时使用）
▸ Step2: Bowtie2 回比对（算覆盖度）
bowtie2-build 先给组装结果建索引
samtools sort 按坐标排序 BAM（分箱需要排序后的文件）
▸ Step3: MetaBAT2 分箱
jgi_summarize_bam_contig_depths 从 BAM 算每条 contig 的覆盖深度
metabat2 -i ... -a ... -i contigs 文件；-a 覆盖度文件；-o 输出 bin 前缀
▸ Step4: CheckM 质控
lineage_wf 全自动工作流：根据谱系标记基因评估完整度和污染度
-x fa bin 文件扩展名（MetaBAT2 默认输出 .fa）
MEGAHIT 组装输出（真实 log 格式）：
--- [STAT] contigs ---
contigs >= 0 bp: 345,678 ← 所有 contigs 数（包括极短的）
contigs >= 1000 bp: 156,789 ← ≥1kb 的 contigs（通常只用这些）
contigs >= 2500 bp: 43,218 ← 分箱一般要求 ≥2.5kb
Total length (>= 0 bp): 456,789,123 bp ← 总组装长度 ~457 Mbp
Total length (>= 1000 bp): 312,456,789 bp
N50: 12,345 bp ← 50% 的序列长度超过此值（宏基因组通常 5k-50k）
N90: 2,456 bp ← 90% 的序列长度超过此值
Largest contig: 89,234 bp ← 最长 contig
GC content: 45.2% ← 整体 GC 含量
─────────────────────────────────────────────
⚠ N50 不像单基因组那样大（单基因组可达 Mbp 级）
⚠ 下游分箱通常只用 ≥1500bp 或 ≥2500bp 的 contigs
CheckM lineage_wf 输出（真实 qa 表格格式）：
Bin Id Marker lineage #genomes #markers #sets 0 1 2 3 4 5+ Completeness Contamination Strain het.
↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑
分箱ID 最近的参考谱系 参考基因组数 标记基因数 基因集 缺失 单拷贝 2拷贝 3拷贝 4拷贝 5+拷贝 完整性% 污染% 株内异质性%
─────────────────────────────────────────────────────────────────────────────────────
bin.1 o__Bacteroidales 120 280 186 3 271 6 0 0 0 97.2 1.3 0.00 ← 高质量 ✅
bin.2 f__Ruminococcaceae 89 340 223 18 312 8 2 0 0 92.5 3.1 0.50 ← 高质量 ✅
bin.3 o__Clostridiales 234 281 193 42 221 15 3 0 0 78.9 6.2 1.20 ← 中质量 ⚠️
bin.4 k__Bacteria 5449 104 58 45 41 12 4 2 0 45.3 12.8 3.40 ← 低质量 ❌ 丢弃
─────────────────────────────────────────────────────────────────────────────────────
⚠ 列 0-5+ = 标记基因的拷贝数分布（理想：大部分在"1"列，即单拷贝）
⚠ "0"列多 → 缺失基因多 → 完整性低；"2+"列多 → 混了多个菌 → 污染高
⚠ MIMAG 标准：高质量 = 完整性>90% + 污染<5%（+16S/23S/5S/tRNA≥18）
MAG 质量分级统计
质量等级标准本次获得典型用途
高质量完整性>90% 污染<5%3 个 MAGs基因组注释、进化分析、代谢重建
中质量完整性>50% 污染<10%5 个 MAGs物种鉴定、基因检索
低质量不满足上述标准12 个 MAGs通常丢弃
组装质量指标
指标本次结果参考范围说明
Contigs 数156,78910万-50万越少越好（拼得越完整）
N5012,345 bp5k-50k bp50% 序列长度超过此值
总长度312.5 Mbp100-500 Mbp代表测到的基因组总量
最长 contig89,234 bp50k-500k bp越长说明组装越好
  • 组装分箱不是必做步骤：reads-based 分析（Kraken2/HUMAnN）就够发文章了
  • 组装消耗大量内存：MEGAHIT 需 ~30-100GB RAM，metaSPAdes 更多
  • 多样本共组装 (co-assembly)：多个样本一起组装可提升 N50，但需要更多内存
  • 分箱结果要用多个工具：MetaBAT2 + MaxBin2 + CONCOCT → DAS Tool 整合最优
  • CheckM2 是新版：CheckM 经典但慢，CheckM2 用机器学习更快，结果相似
  • N50 不是越大越好：宏基因组 N50 通常几千到几万 bp，不像单基因组那么长
组装+分箱不是必做步骤，reads-based 分析（Kraken2/HUMAnN）就够发文章了。分箱的意义是获得"准基因组"（MAG），可以做基因组级别分析。面试记住 MAG 质量标准：高质量 = 完整性>90% + 污染<5%（MIMAG标准）。
06 统计分析与可视化 (Statistical Analysis) 出图发文 ▶

拿到物种/功能丰度表后，做多样性分析、差异分析、关联分析，画图写论文。这是"讲故事"的阶段。

输入文件

species_abundance.tsv（物种丰度矩阵）
pathway_abundance.tsv（通路丰度矩阵）
metadata.tsv（样本分组信息）
→

输出文件

alpha_diversity.pdf（Shannon/Simpson）
beta_diversity_PCoA.pdf（PCoA/NMDS图）
differential_species.tsv（差异物种表）
LEfSe_results.pdf（生物标志物）
heatmap.pdf（丰度热图）
barplot_phylum.pdf（门水平堆叠柱状图）
R: vegan（多样性） R: phyloseq R: DESeq2（差异分析） Python: scikit-bio LEfSe（生物标志物） MaAsLin2（关联分析） R: ggplot2（画图）
# R 语言：α多样性计算
library(vegan)
shannon <- diversity(otu_table, index = "shannon")
simpson <- diversity(otu_table, index = "simpson")

# β多样性 PCoA
bray_dist <- vegdist(otu_table, method = "bray")
pcoa_res <- cmdscale(bray_dist, eig = TRUE, k = 2)

# PERMANOVA 检验（组间差异是否显著）
adonis2(bray_dist ~ Group, data = metadata)

# Wilcoxon 差异分析
wilcox.test(abundance ~ group, data = species_data)
📖 参数逐行讲解
▸ α多样性（样本内多样性）
diversity(otu_table, index = "shannon") 计算 Shannon 指数：值越大物种越丰富且均匀（肠道一般 2.5-4.0）
diversity(otu_table, index = "simpson") 计算 Simpson 指数：值越接近 1 多样性越高（侧重优势种占比）
▸ β多样性（样本间差异）
vegdist(otu_table, method = "bray") 计算 Bray-Curtis 距离矩阵：考虑丰度差异，宏基因组最常用
cmdscale(bray_dist, eig = TRUE, k = 2) PCoA 降维：把高维距离投影到 2D 平面画散点图
▸ 统计检验
adonis2(bray_dist ~ Group, data = metadata) PERMANOVA 检验：判断分组是否能解释样本间差异（看 R² 和 p 值）
wilcox.test(abundance ~ group) Wilcoxon 秩和检验：非参数方法比较两组差异（不要求正态分布，微生物数据首选）
R vegan 多样性计算（真实代码 + 输出格式）：
> library(vegan)
> shannon <- diversity(otu_table, index = "shannon")
> shannon
  sample1 sample2 sample3 sample4 sample5 sample6
  3.412 3.567 3.289 2.891 2.745 2.934 ← 每个样本一个 Shannon 值
  (Health) (Health) (Health) (T2D) (T2D) (T2D) ← 健康组明显高于 T2D 组
─────────────────────────────────────────────
> wilcox.test(shannon ~ group, data = metadata)
  Wilcoxon rank sum test ← 非参数检验（不假设正态分布）
  W = 89, p-value = 0.003 ← p<0.05 表示两组差异显著
─────────────────────────────────────────────
> adonis2(bray_dist ~ Group, data = metadata) ← PERMANOVA 检验
   Df SumOfSqs R2 F Pr(>F)
  Group 1 0.8234 0.156 5.23 0.001 ← R²=0.156 表示分组解释了 15.6% 的变异
  Residual 28 4.4567 0.844 ← 剩余 84.4% 由其他因素解释
  Total 29 5.2801 1.000
─────────────────────────────────────────────
⚠ Shannon 值范围通常 1.5-5.0，肠道菌群一般 2.5-4.0
⚠ PERMANOVA 的 R² 在菌群研究中通常 5-20%，不要期望太高
⚠ p 值必须做 FDR 校正：p.adjust(p_values, method="BH")
Wilcoxon 差异分析 + FDR 校正（真实输出格式）：
Species Health(%) T2D(%) log2FC p-value p.adjust(BH) Direction
mean±SD mean±SD ↑ 原始p FDR校正后p
─────────────────────────────────────────────────────────────────────────
Faecalibacterium prausnitzii 8.2±2.1 3.1±1.5 -1.40 0.00012 0.0018 ** ↓ T2D下降 ← 产丁酸菌，核心标志物
Roseburia intestinalis 5.6±1.8 2.3±1.2 -1.28 0.00034 0.0034 ** ↓ T2D下降 ← 产丁酸菌
Escherichia coli 1.2±0.6 4.8±2.3 +2.00 0.00021 0.0025 ** ↑ T2D升高 ← 条件致病菌！LPS→炎症
Desulfovibrio sp. 0.5±0.3 2.1±0.9 +2.07 0.00089 0.0067 ** ↑ T2D升高 ← 硫酸盐还原菌，产H₂S
Akkermansia muciniphila 3.4±1.1 1.2±0.7 -1.50 0.0012 0.0120 * ↓ T2D下降 ← 改善代谢，热门研究对象
─────────────────────────────────────────────────────────────────────────
⚠ log2FC: 负值=T2D中减少 正值=T2D中增加（|log2FC|>1 通常有生物学意义）
⚠ p.adjust(BH) 才是最终看的 p 值！原始 p-value 有多重检验膨胀
⚠ 相对丰度变化 ≠ 绝对丰度变化（组成数据的陷阱！）
α多样性 — 健康组 vs T2D 组
指标健康组 (n=30)T2D 组 (n=30)p-value结论
Shannon3.45 ± 0.322.89 ± 0.450.003 **T2D 多样性显著降低
Simpson0.92 ± 0.030.85 ± 0.060.008 **T2D 均匀度下降
Observed OTUs456 ± 67312 ± 890.001 ***T2D 物种数减少
差异物种 TOP5（T2D vs 健康）
物种健康组T2D组变化生物学意义
F. prausnitzii8.2%3.1%↓ 下降产丁酸，抗炎
Roseburia5.6%2.3%↓ 下降产丁酸
A. muciniphila3.4%1.2%↓ 下降改善代谢
E. coli1.2%4.8%↑ 升高条件致病菌
Desulfovibrio0.5%2.1%↑ 升高硫酸盐还原菌，促炎
常用统计方法速查
分析类型方法R 函数适用场景
α多样性比较Wilcoxon 秩和检验wilcox.test()两组比较
α多样性比较Kruskal-Walliskruskal.test()三组及以上
β多样性检验PERMANOVAadonis2()组间差异是否显著
β多样性可视化PCoA / NMDScmdscale() / metaMDS()样本聚类散点图
差异物种Wilcoxon + FDRwilcox.test() + p.adjust()非参数，最常用
差异物种DESeq2DESeqDataSetFromMatrix()参数法，适合小样本
生物标志物LEfSe (LDA)lefse 命令行工具筛选 LDA > 2 的标志物
关联分析MaAsLin2Maaslin2()菌群与临床指标关联
  • 多重检验校正必做：几百个物种同时检验，必须做 FDR（BH）校正，否则假阳性爆炸
  • 非参数优先：微生物丰度数据通常不符合正态分布，优先用 Wilcoxon 而非 t-test
  • β多样性距离选择：Bray-Curtis（考虑丰度）最常用，UniFrac（考虑进化关系）需要进化树
  • 样本量要够：每组至少 10 个样本，理想 30+ 个，样本太少统计功效不足
  • 混杂因素：年龄、性别、BMI、饮食等都影响菌群，分析时要考虑协变量
  • 相对丰度的陷阱：一个菌升高不代表绝对量增加，可能是其他菌减少导致的"被动升高"
α多样性 = 单个样本内的丰富度（Shannon指数越高越好）。β多样性 = 样本间的差异（用PCoA/NMDS图展示，PERMANOVA检验显著性）。差异分析常用 Wilcoxon（非参数）、DESeq2（参数）、LEfSe（LDA效应值）。你的毕设用的就是这套分析！

Kraken2 vs MetaPhlAn 详细对比

对比维度 Kraken2 MetaPhlAn4
原理 k-mer 精确匹配（把read切成k个碱基的小段，查数据库） Marker gene 比对（用~100万个物种特异性标记基因）
速度 极快（分钟级，100万reads/分钟） 较慢（小时级，需要Bowtie2比对）
精度 假阳性较多（数据库里没有的也可能匹配到） 更准（只用验证过的marker gene）
数据库大小 标准库 ~50GB（可用迷你库 ~8GB） ~20GB（相对小）
内存需求 高（标准库需 ~64GB RAM） 低（~8GB RAM）
输出丰度 reads 数（需要 Bracken 校正成丰度） 直接输出相对丰度（%，加起来=100）
适用场景 大规模快速筛选、临床实时检测 发文章、精确物种定量
论文使用 近年增长快，引用多 顶刊标配，Nature/Cell 常用
面试回答："两个工具各有优势，Kraken2 快适合大批量筛选，MetaPhlAn 准适合发文章。实际项目中我会两个都跑，互相验证。"

16S 扩增子 vs 宏基因组测序

对比维度 16S 扩增子 宏基因组（Shotgun）
原理 只测 16S rRNA 基因的一小段 随机打断测所有 DNA
分辨率 通常到属（Genus） 可到种甚至株（Species/Strain）
功能信息 无（只有分类信息） 有（基因、通路、抗性基因等）
成本 低（~200元/样本） 高（~2000元/样本）
数据量 ~10万条 reads/样本 ~5000万条 reads/样本
分析复杂度 简单（QIIME2一站式） 复杂（多工具组合）
适用场景 预算少、样本多、只看"有哪些菌" 需要功能信息、需要高分辨率、追溯耐药基因
典型工具 QIIME2, DADA2, mothur Kraken2, MetaPhlAn, HUMAnN, MEGAHIT

面试万能表达

描述流程时：
"宏基因组分析的标准流程是：原始数据质控、去宿主、物种注释、功能注释，如果需要还可以做组装分箱获得MAG，最后统计分析出图。"
描述工具选择：
"物种注释我主要用 MetaPhlAn4，它基于 marker gene 精度更高；大规模数据会配合 Kraken2 做快速筛选。功能注释用 HUMAnN3 看代谢通路。"
描述质控标准：
"质控后 Q30 > 85%，去宿主后统计微生物 reads 占比。MAG 质量用 CheckM 评估，高质量标准是完整性 > 90% 且污染 < 5%。"
描述你的毕设项目：
"我的毕设是 2 型糖尿病肠道菌群研究，用 Kraken2 做物种注释，发现 T2D 患者产丁酸菌（如 F. prausnitzii）显著减少，条件致病菌富集，然后用随机森林筛选生物标志物。"