第三章：输入数据格式要求——FASTQ、FASTA、BAM的适用场景与预处理注意事项

说到MetaPhlAn的输入格式，我见过不少人在第一步就栽跟头。其实选对格式，你的分析就已经成功了一半。我个人习惯在项目启动前，先理清楚手里的数据到底是什么来路。

三种格式的核心区别

咱们先看个直观的对比表，心里有个底：

格式 包含信息 典型来源 MetaPhlAn是否原生支持
FASTQ 序列 + 碱基质量分数 Illumina / MGI 测序仪下机数据 是，最推荐
FASTA 纯序列，无质量信息 组装结果、公共数据库下载 是，但丢失质控维度
BAM 比对结果（序列 + 比对位置 + 质量） 先比对到参考基因组后的输出 是，需要排序并建索引
我的一点心得： 如果你是直接从NCBI的SRA下载数据，fastq-dump出来的通常是FASTQ。但如果你拿到的数据是别人预处理好的，尤其是一些“整洁”的FASTA文件，那就要多留个心眼——这些数据可能已经丢掉了质量信息，后续无法进行有效的低质量过滤。

FASTQ 格式：最推荐的起点

FASTQ是MetaPhlAn处理的首选。为什么？因为它包含了碱基质量分数。这玩意儿在后续的质控和比对中非常关键。你想想看，如果一个碱基的测序质量是10分（错误率10%），那它和另一个质量30分（错误率0.1%）的碱基，在比对时的可信度能一样吗？

适用场景：

  • 原始或轻度预处理（去接头、去低质量）的鸟枪法宏基因组测序数据。
  • 需要做完整的质控流程（FastQC + Trimmomatic/Trimmomatic）的数据。
  • 你希望保留最大的分析灵活性。

预处理注意事项：

嗯，这里要注意，不是拿到FASTQ就直接丢给MetaPhlAn。我建议先做一轮轻量级质控：

# 1. 检查原始质量
fastqc sample_R1.fastq.gz sample_R2.fastq.gz

# 2. 去除低质量碱基和接头（以Trimmomatic为例）
java -jar trimmomatic.jar PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_R1_paired.fastq.gz sample_R1_unpaired.fastq.gz \
  sample_R2_paired.fastq.gz sample_R2_unpaired.fastq.gz \
  ILLUMINACLIP:adapters.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

我曾经遇到一个项目，用户直接拿原始FASTQ跑MetaPhlAn，结果菌种鉴定率极低。排查后发现测序数据里有大量低质量read，全被Bowtie2比对阶段过滤掉了。所以，质控这一步不能省。

特别提醒： MetaPhlAn 4.0.6及以上版本要求FASTQ文件必须为未压缩或gzip压缩的格式。如果你遇到`.bz2`或`.xz`后缀，记得先转成`.gz`。

FASTA 格式：当原始数据不可得时

FASTA文件是纯序列文件，没有每个碱基的质量分数。它的适用场景比较窄：

  • 你从公共数据库获取的是已组装的contigs（比如MEGAHIT、metaSPAdes的输出）。
  • 你手里只有他人分析完后分享的FASTA文件，原始FASTQ已经丢失。
  • 你只想快速验证某条序列是否存在，不做定量分析。

预处理注意事项：

用FASTA直接跑MetaPhlAn，说白了就是放弃了质控的可能性。我个人习惯在做之前加一步：

# 建议：过滤掉长度小于 500bp 的序列
awk '/^>/ {if (seq_len >= 500) print header; printf "%s", seq; print ""} 
     {if (/^[A-Za-z]/) seq_len = length($0); else seq_len = 0; header = $0}' 
     input.fasta > filtered.fasta

为什么过滤短序列？因为短序列（尤其是 <200bp）在MetaPhlAn的比对过程中，很容易出现多匹配或匹配到非特异的区域，导致假阳性。我在做小鼠肠道宏基因组时，就因为没过滤短contigs，结果检出了几十个低丰度的环境菌，实际上那些都是比对错误。

关键点： 如果你必须用FASTA，请确保序列描述行（以`>`开头）是干净的——不要包含空格或特殊字符，否则MetaPhlAn在解析时可能会报错。最好只保留序列ID，比如：
>scaffold_1 ✅
>scaffold_1 length=1234 coverage=5.6x ❌

BAM 格式：进阶用户的提速利器

BAM是已经比对到参考基因组或标记基因集的结果文件。用BAM作为输入，MetaPhlAn可以跳过最耗时的比对步骤，直接从profile丰度估计开始。适合跑大规模重复分析（比如重跑不同版本的数据库）。

适用场景：

  • 你之前已经用相同的数据库跑过一次比对，保存了BAM文件。
  • 你想快速测试不同的参数（比如--add_viruses、--unknown_estimation）。
  • 你需要做降采样分析（只用部分数据快速验证）。

预处理注意事项：

BAM文件不是拿来就能用的。MetaPhlAn要求：

# 1. 排序（按染色体位置）
samtools sort input.bam -o sorted.bam

# 2. 建索引
samtools index sorted.bam

# 3. （可选）过滤掉未比对上或低质量的read
samtools view -b -F 260 -q 30 sorted.bam > filtered.bam

# 4. 重新索引
samtools index filtered.bam

这里有个容易踩的坑：BAM文件必须是成对end的，并且read1和read2要成对出现。如果BAM里混入了未成对的singleton，MetaPhlAn会报错。我之前处理一个来自HMP项目的公开数据时，就因为BAM里混了少量singleton，折腾了大半天才排查出来。

避坑指南： 千万别直接用samtools view -bS把SAM转成BAM就丢给MetaPhlAn。我曾经这么干过，结果跑了20分钟报错，一看日志是“BAM is not sorted”。排序+建索引是硬性要求，缺一不可。

格式选择决策指南

说了这么多，到底该选哪个？我总结了一张决策图，你看一眼就明白了。

MetaPhlAn 输入格式选择决策流程 你手里的原始数据 是否同时有 R1 和 R2 FASTQ？ 推荐：FASTQ 质控 > MetaPhlAn 是否有现成的 BAM？ 推荐：BAM 排序+索引 > MetaPhlAn 只有 FASTA 组装结果？ 可用：FASTA 过滤短序列 > MetaPhlAn 常见误区 ❌ BAM 不做排序索引直接用 ❌ FASTA 不做长度过滤 ❌ 以为所有格式跑出来结果一样

最严苛的预处理检查清单

不管你选哪种格式，这几点最好落地：

  1. 文件名统一：不要用中文、空格、特殊符号。我习惯用项目名_样本编号_R1.fastq.gz这种模式。
  2. 文件完整性校验：下载或传输后，务必做md5sum检查。我在HPC集群上遇到过传输中断导致文件损坏的案例，跑了两天报错，重来。
  3. 环境变量问题：如果你在集群上使用module load加载samtools/bowtie2，记得确认版本一致。不同版本的索引格式可能不兼容。
  4. 数据存储路径：MetaPhlAn运行过程中会生成临时文件（默认在/tmp）。如果/tmp空间不足（比如只有10GB），会直接导致OOM。我建议在命令里指定：
# 指定临时目录到数据盘
export TMPDIR=/scratch/$USER/tmp
mkdir -p $TMPDIR
metaphlan sample.fastq.gz --input_type fastq -o profile.txt --tmp_dir $TMPDIR
一句话总结： 优先用FASTQ，因为它保留了最完整的原始信息。BAM适合快速复跑，但排序和索引必须做。FASTA是备选方案，使用时需额外留意序列质量。不管什么格式，预处理阶段多花10分钟，能省下后续几小时甚至几天的排错时间。
← 上一章 📖 返回目录 下一章 →