3、输入数据准备：FastQ文件的格式与质量控制
MetaPhlAn跑出来的结果靠不靠谱，很大程度取决于你喂进去的数据。我见过不少同学拿着原始测序数据直接开跑，结果物种丰度图里冒出来一堆奇奇怪怪的微生物……说白了，就是低质量数据和接头污染在捣乱。
这一章，我们把数据准备这件事聊透。从FastQ的格式细节，到质量控制的标准流程，再到我踩过的那些坑，一并说清楚。
3.1 FastQ格式：不止是四行文字
你可能觉得FastQ就是四行一组的数据，没什么好讲的。但我们在实际项目中，至少有30%的分析异常跟对格式的理解偏差有关。
一条标准的FastQ记录长这样：
@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
CTCAAGGTTGTTGCAAGACGGGAGGTAGGTGTCTGGGTAGGGTAGAGTT
+
JJJJJIIJJJJJIHHHHHFFFFFEEEEEEEEEDDDDDDDDDDDDDDDD四行分别对应什么？
| 行号 | 标识 | 内容说明 | 常见坑点 |
|---|---|---|---|
| 第1行 | @ | 序列标识符及描述信息 | 部分平台的描述字段包含空格，解析时容易出错 |
| 第2行 | - | 碱基序列（A/T/C/G/N） | 长读长数据可能出现换行截断 |
| 第3行 | + | 分隔符，有时重复描述信息 | 不能为空行，否则下游工具直接报错 |
| 第4行 | - | 质量评分字符串（ASCII编码） | 长度必须等于第2行序列长度 |
质量值是怎么编码的？ 每个字符对应一个Phred质量分数，公式是 Q = -10 * log10(P)，其中P代表该碱基测序错误的概率。字符与分数的映射关系取决于版本：
- Sanger / Illumina 1.8+：ASCII 33~126，对应Q值0~93
- Illumina 1.3~1.7：ASCII 64~126，对应Q值0~62
- Illumina < 1.3：ASCII 59~126，对应Q值0~62
我个人习惯在拿到新数据后，先用 head -n 40 看两眼质量值那行。如果看到一大堆 B、? 这类低质量字符占主导，就要当心了——要么是测序质量不好，要么是编码格式搞错了。
为什么强调这个？我曾在项目里拿到一批公共数据库的RNA-seq数据，用FastQC怎么检查都提示质量异常。后来发现是数据用了Illumina 1.3格式，但我们的质控工具默认用的是Sanger偏移量。嗯，这里要注意——编码格式不匹配会导致质量分数量级错误，过滤阈值完全失效。
3.2 质控流程：从原始数据到可靠输入
MetaPhlAn对输入数据的要求很简单——高质量、无接头、尽可能保留有效序列。但实际操作起来，每个环节都有门道。
完整的质控流程我一般拆成三步：
- 质量评估：先用FastQC摸清数据底细
- 过滤与修剪：切除低质量碱基和接头序列
- 数据核对：确认过滤后的数据是否满足下游要求
下面这张图是我经常用来给团队成员讲流程的，你一看就明白：
3.3 实战：用FastQC做质量评估
FastQC是我工具箱里的常客。运行起来很简单：
fastqc sample_R1.fastq sample_R2.fastq -o ./qc_output/ -t 4参数说明：
sample_R1.fastq和sample_R2.fastq：双端测序的两个文件-o ./qc_output/：指定输出目录-t 4：使用4个线程并行处理
运行结束后会生成一个HTML报告。我重点看这几个模块：
| 模块 | 关注指标 | 合格标准 |
|---|---|---|
| Per base sequence quality | 每条read上每个位置的Q值分布 | 整体Q>30，尾部Q>20 |
| Per sequence GC content | GC含量分布是否接近理论分布 | 无明显多峰或异常偏移 |
| Adapter Content | 接头序列在各位置的富集程度 | 任意位置<5% |
| Overrepresented sequences | 高度重复的异常序列比例 | 单个序列占比<0.1% |
曾经踩过的坑：有一次我看到Per base sequence quality的末尾一塌糊涂，原本以为是测序质量问题。后来发现是双端测序的R1和R2文件传反了——R2的质量在尾部天然比R1差。所以拿到数据先看一眼基本统计，别急着下结论。
3.4 过滤与修剪：Trimmomatic实战
我个人偏好用Trimmomatic做双端数据的过滤，因为它的ILLUMINACLIP功能可以同时处理多种接头。下面是一个常用的配置：
java -jar trimmomatic-0.39.jar PE \
sample_R1.fastq sample_R2.fastq \
sample_R1_trimmed.fastq sample_R1_unpaired.fastq \
sample_R2_trimmed.fastq sample_R2_unpaired.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
LEADING:3 \
TRAILING:3 \
SLIDINGWINDOW:4:20 \
MINLEN:36每个参数在做什么？
| 参数 | 含义 | 我建议的阈值 |
|---|---|---|
ILLUMINACLIP |
接头序列比对与切除 | 2:30:10（错配数:比对分数阈值:最小匹配长度） |
LEADING:3 |
切除5'端质量低于3的碱基 | 通常3~5，看数据质量调整 |
TRAILING:3 |
切除3'端质量低于3的碱基 | 与LEADING保持一致 |
SLIDINGWINDOW:4:20 |
4碱基滑窗，平均质量低于20则切除 | 对于Metagenomic数据，20~25比较合理 |
MINLEN:36 |
丢弃长度小于36的read | 太短会影响MetaPhlAn的k-mer比对 |
小技巧：不用每次都输这么长的命令。我习惯把参数写到一个配置文件里，或者用Makefile统一管理。你想一想，如果处理100个样本，每次手动敲命令，光检查参数就会疯掉。
3.5 过滤后的数据核对
跑完Trimmomatic别急着进MetaPhlAn，再做一轮Quick检验：
- 检查输出日志：Trimmomatic会打印过滤统计，看看丢弃比例是否合理（通常<20%算正常，超过30%要查数据来源）
- 对比FastQC前后报告：确认接头污染被清除，碱基质量分布改善
- 查看序列数量：确保双端文件行数一致（用
wc -l快速核对）
核对命令示例：
# 统计FastQ行数
echo "R1 lines:" && wc -l sample_R1_trimmed.fastq
echo "R2 lines:" && wc -l sample_R2_trimmed.fastq
# 查看前几条记录确认格式正确
head -n 8 sample_R1_trimmed.fastq曾经踩过的坑：有一回Trimmomatic跑完，我看两个trimmed文件的行数居然不匹配。排查了好久才发现是磁盘空间不够，导致输出文件写到一半中断了。从那以后，我给每批数据都加了行数校验，并且在脚本里嵌入磁盘空间检查——你想想看，跑了一整晚的作业因为磁盘满了而失败，多冤。
3.6 质控后的数据，到底什么样才算“好”？
MetaPhlAn对输入数据其实挺宽容的，但我对数据入手的门槛从不含糊。下面是我的验收清单：
- 每条read的读长：双端测序建议保留至少75bp以上；对于MetaPhlAn的k-mer匹配，太短的read容易产生假阳性
- 平均质量：Q≥30（也就是碱基错误率≤0.1%）
- 接头残留：FastQC的Adapter Content模块应显示为绿色（通过）
- 双端序列对数：经过过滤后，paired reads数应不低于原始数据的70%
如果你发现过滤比例太高，先别急着调低阈值。我在项目里遇到过样本质量差导致过滤比例高达40%的情况——强行保留低质量数据，只会让MetaPhlAn的谱系注释引入更多错误。