第4节 宏基因组分析基础：鸟枪法测序数据质控、去宿主污染、组装、基因预测与丰度定量

宏基因组分析的第一步，说白了就是把测序仪吐出来的原始数据，变成我们能用的基因列表和丰度表。这个过程绕不开几个关键环节：质控、去宿主、组装、基因预测和丰度定量。今天我就把这几个步骤拆开揉碎了讲，里面有不少我在项目中踩过的坑，你注意别重蹈覆辙。

宏基因组分析核心流程 原始数据 质控 去宿主 组装 基因预测 丰度定量 产出：基因丰度矩阵

4.1 鸟枪法测序数据质控

原始测序数据直接拿来做分析？千万别急。我见过太多新手上来就跑组装，结果一半是低质量序列和接头污染，组装出来的contig惨不忍睹。质控这步省不了，也急不得。

常用的工具有 FastQC 和 fastp。FastQC 用来做质量评估，fastp 能一站式完成过滤、裁剪和报告生成。我个人习惯先用 FastQC 扫一眼全局质量，再用 fastp 做精细化处理。

实用技巧：fastp 默认参数对大多数宏基因组数据已经够用。但如果你的数据来自极端环境（比如高GC菌群），建议调整 --cut_front 和 --cut_tail 的阈值。

下面是我常用的 fastp 命令，你直接拿去用就行：

fastp -i R1.fastq.gz -I R2.fastq.gz \
      -o R1.clean.fastq.gz -O R2.clean.fastq.gz \
      --detect_adapter_for_pe \
      --cut_front --cut_tail \
      --cut_front_mean_quality 20 \
      --cut_tail_mean_quality 20 \
      --qualified_quality_phred 20 \
      --length_required 50 \
      --thread 8 \
      --json fastp.json --html fastp.html

跑完之后，重点看 Q20/Q30 比例 和 clean reads 占比。我一般要求 Q30 不低于 85%，clean reads 占比低于 70% 就要回头检查是不是建库出了问题。

4.2 去宿主污染

宏基因组样本里经常混有宿主DNA，尤其是临床样本或肠道内容物。去宿主不干净，后续组装出来的contig里会掺大量宿主序列，直接拉低基因预测的准确率。

标准的做法是用 Bowtie2 把 reads 比对到宿主基因组，然后保留未比对的 reads。我经历过一次惨痛教训——用了一个不完整的宿主参考基因组，结果去宿主不彻底，下游差异分析全乱了。

避坑指南：宿主参考基因组一定要选最新版本，并且建议用 combined 模式（比如人类基因组用 T2T-CHM13 + GRCh38 的补丁）。我用的是 Ensembl 或 NCBI 的官方 release。

构建索引和比对命令示例：

# 构建宿主索引
bowtie2-build host_genome.fasta host_index

# 去宿主比对
bowtie2 -x host_index \
        -1 R1.clean.fastq.gz \
        -2 R2.clean.fastq.gz \
        --very-sensitive \
        -p 16 \
        --un-conc-gz host_removed \
        -S aligned.sam

比对结束后，host_removed.1.fastq.gz 和 host_removed.2.fastq.gz 就是去宿主后的 clean 数据。你可以在日志里看一眼比对率，如果宿主污染超过 40%，说明采样环节可能需要优化。

4.3 宏基因组组装

组装的目的是把短 reads 拼成更长的 contig，方便后续基因预测。宏基因组组装比单菌基因组组装难得多，因为群体里物种多样性高、丰度不均一。目前社区里主流工具就两个：MEGAHIT 和 metaSPAdes。

我个人习惯这么选：数据量大、计算资源有限时用 MEGAHIT；追求组装质量、不差机器时用 metaSPAdes。这不是绝对的，但我在 50+ 个项目里试下来，这个策略最稳。

工具 特点 适用场景 内存需求
MEGAHIT 极快，内存优化好，适合复杂群落 大数据量、普通服务器 中等（~64 GB）
metaSPAdes 组装质量高，scaffold 更长，但慢 小数据量、高精度需求 高（~256 GB）

MEGAHIT 使用示例：

megahit -1 host_removed.1.fastq.gz \
        -2 host_removed.2.fastq.gz \
        -o megahit_output \
        --num-cpu-threads 16 \
        --memory 0.8 \
        --k-min 27 --k-max 127 --k-step 10

metaSPAdes 使用示例：

metaspades.py -1 host_removed.1.fastq.gz \
              -2 host_removed.2.fastq.gz \
              -o metaspades_output \
              -t 16 \
              -m 200

关键指标：组装完成后评估 N50 和 总组装长度。宏基因组组装 N50 通常 1-10 Kb 都算正常，低于 1 Kb 说明数据质量或复杂度有问题。我见过最差的 N50 只有 300 bp，最后排查发现是去宿主时参数设太严把菌群序列也滤掉了。

4.4 基因预测与丰度定量

拿到 contig 之后，下一步是找出里面的编码基因。这一步我用 Prodigal，没有之一。它专门为细菌和古菌基因组设计，支持宏基因组的 metagenomic 模式，能自动处理不完整的基因片段。

基因预测命令很简单：

prodigal -i final.contigs.fa \
         -o genes.gff \
         -a proteins.faa \
         -d genes.fna \
         -p meta \
         -q

输出文件里，genes.fna 是核苷酸序列，proteins.faa 是氨基酸序列。我个人建议把 genes.gff 也保留下来，后面做功能注释时能用上坐标信息。

基因预测完成后，就到了定量环节。这里我用 Salmon，它基于转录本定量思路，速度快、内存低，而且不依赖比对，直接通过 k-mer 计数估算丰度。

经验之谈：Salmon 的索引构建非常快，但要注意 —k 参数默认 31，如果你的 reads 长度只有 100 bp 或更短，建议降到 23-27，否则敏感度会下降。我吃过这个亏，换了参数后基因检出率提升了 12%。

先建索引，再做定量：

# 构建索引
salmon index -t genes.fna \
             -i salmon_index \
             -k 31

# 定量
salmon quant -i salmon_index \
             -l A \
             -1 host_removed.1.fastq.gz \
             -2 host_removed.2.fastq.gz \
             -o salmon_output \
             -p 16 \
             --validateMappings

定量结果在 salmon_output/quant.sf 里。核心列是 TPM 和 NumReads。TPM 用于样本间比较，NumReads 用于绝对丰度计算。

注意：Salmon 的 -l A 参数会自动判断文库类型，但如果你对文库类型有把握，手动指定 -l ISF 或 -l IU 能提高准确性。我遇到过自动判断出错的情况，尤其当数据存在 3' 偏倚时。

到此为止，你已经从原始测序数据走到了基因丰度矩阵。这一整套流程，我在不同项目里跑过上百次。每次出问题，80% 都出在质控和去宿主这两个前期环节上——前期的数据质量不过关，后面的组装、预测、定量再努力也是白费。嗯，这个道理你记住了，后面的分析就不会太折腾。

← 上一章 📖 返回目录 下一章 →