数据质量控制：从原始数据到可信分析的第一步
数据质量控制，说白了就是给你的原始测序数据做一次全面体检。我早期在项目组里吃过亏——那时拿到一批RNA-seq数据，没仔细看QC就直接跑比对，结果发现样本间批次效应严重，回查才发现测序质量本来就参差不齐。
所以，我把QC视为整个验证方法论中最前置、最不可跳过的一环。你的pipeline跑得再快、再优雅，如果输入数据本身有问题，后续的输出全是错的。
FastQC：行业标准的质量显微镜
FastQC是目前使用最广泛的原始测序数据质量评估工具。它不依赖参考基因组，完全基于FASTQ文件自身的内容给出质量报告。
我建议每个生信工程师都把FastQC当作写pipeline的第一道检查站。什么时候该加？我个人习惯是——每一步数据变换（如trimming、deduplication）之前和之后都要跑一遍。
# 基础用法
fastqc sample_R1.fastq.gz sample_R2.fastq.gz --outdir ./qc_output你想想看，这个命令会生成什么？它会产出每个样本的HTML报告和ZIP压缩包。HTML可以直接用浏览器打开，ZIP包里是机器可读的原始数据，方便下游脚本做批量判断。
常用的参数我个人用得很稳定：
fastqc \
--outdir ./qc_reports \
--threads 8 \
--kmers 7 \
--contaminants contaminant_list.txt \
sample_*.fastq.gz- --threads：线程数，通常设为物理核心的50%
- --kmers：k-mer分析长度，默认7，如果需要更精细的overrepresented sequence，可以调大
- --contaminants：用户自定义的污染序列库，比如接头、rRNA序列
FastQC给出11项评估指标，但我认为最重要的几个需要我们持续关注：
| 模块名称 | 评估内容 | 核心判断依据 |
|---|---|---|
| Per base sequence quality | 每个碱基位置的质量得分分布 | Q30比例、中位数趋势线是否下沉 |
| Per sequence quality scores | 所有序列的平均质量分布 | 峰形是否左偏、尾部是否有低质量次峰 |
| Per base N content | 每个碱基位置N碱基的比例 | 任何位置超过5%即视为严重警告 |
| Sequence Length Distribution | 读长分布范围 | 单峰正常、多峰可能存在建库问题 |
| Sequence Duplication Levels | 序列重复率 | 富集文库在合理范围内取重复 |
| Overrepresented sequences | 超高丰度序列 | 接头污染、rRNA残留或扩增偏好 |
| Adapter Content | 接头残余量 | 3'端接头是否被有效去除 |
提示：FastQC的质量警告并不是绝对的「失败」。它在做极端保守的判断。比如WGBS的MspI酶切产物会出现特定碱基偏好，这是正常的——不要迷信红叉，要结合实验背景理解。
MultiQC：当样本量超过3个时的必需品
单个FastQC报告不难看。但当你手上有30个样本，每个样本两个端，一共60个HTML文件，你愿意一个个点开看吗？我不愿意。
MultiQC就是解决这个问题的。它能聚合多个FastQC报告，生成一份跨样本的质量总览。我早期做项目时，经常被几十个FastQC报告搞得头疼，自从用上MultiQC，这一步的效率提升了至少一个数量级。
# 安装方式
pip install multiqc# 多目录聚合
multiqc ./project_qc_reports/ \
--outdir ./multiqc_output \
--title "2025-03-21 项目QC总报告" \
--filename project_master_qc.html我个人习惯把 --force 加上，这样每次重新运行会自动覆盖旧报告，避免缓存干扰新版分析结果。
MultiQC输出的总报告非常直观：
- General Statistics：所有样本的关键指标并列，一眼看出离群样本
- Per Base Sequence Quality：所有样本的质量热图或叠加线图
- Sequence Duplication：样本间重复率的对比
- Adapter Content：接头污染在样本维度的分布
注意：MultiQC的聚合依赖于FastQC的原始数据（ZIP包里的fastqc_data.txt）。如果你只保留了HTML报告而删掉了ZIP，MultiQC读不到数据会报错。我建议在存储策略中设置规则——ZIP文件长期保留，HTML按项目清理。
QC指标的验证标准：何时通过、何时报警
现在坐下来，我们说说阈值的事。很多人争论Q30应该大于多少才算合格——我的观点是：没有绝对阈值，但有通用参考线。
| 质量指标 | 警戒线 | 严重问题 | 我的经验值 |
|---|---|---|---|
| Per base Q30比例 | < 80% | < 60% | HiSeq X: >90%; NovaSeq: >85% |
| GC含量偏差 | 偏离理论值 > 10% | 偏离 > 15% | WGS: 41% ± 5%; RNA-seq: 视物种 |
| 重复率 | 文库扩增 > 50% | > 70% | WGS通常 < 30%; RNA-seq可能 < 60% |
| 接头污染 | > 0.1% | > 1% | 大部分测序中心的目标是 < 0.01% |
| N比例 | > 5% 在任何碱基位置 | > 20% | 理想 < 1% |
这里我要提一个大家容易忽略的点：File completeness（文件完整性验证）。我曾经遇到过，一批数据所有FastQC指标都正常，但某个样本的FASTQ文件在传输过程中被截断——文件大小小了三分之一，却没有任何QC模块报警。为什么？因为FastQC只处理它读到的序列，截断后剩余的序列仍是合法格式。
所以，我建议在运行FastQC之前，先做两层基础检查：
# 1. 文件完整性校验
md5sum *.fastq.gz > md5list.txt
# 2. 行数验证（每个read占4行）
zcat sample.fastq.gz | wc -l验证方法论的核心：QC不是跑完FastQC就算完事。你需要将以下验证点写进pipeline的检查逻辑里：
- 格式验证：FASTQ四行结构完整、base质量值在有效范围内
- 统计验证：序列总数、总碱基数、GC含量在预期区间内
- 质量阈值验证：基于预定义的通过/警告/失败三级标准
- 批内一致性验证：同批样本的QC全貌不应出现严重离群
- 历史基线对比：与前几批次同类型数据的QC指标对比
我习惯把这些检查点做成一个独立的shell脚本或Python模块，放在pipeline的开头。如果某条检查没通过，pipeline直接挂掉，并发送日志给负责人。
# 一个简单的QC check函数，我个人用的版本
check_qc() {
local qc_report=$1
local q30_threshold=80
local q30_value=$(grep -A2 "Per base sequence quality" "$qc_report" \
| awk '{for(i=1;iNF;i++) if($i ~ /^[0-9]+/) print $i}' \
| head -1)
if (( $(echo "$q30_value < $q30_threshold" | bc -l) )); then
echo "FAIL: Q30 = $q30_value% < $q30_threshold%"
exit 1
fi
echo "PASS: Q30 = $q30_value%"
}避坑：我曾经犯过一个错误——在Q30通过后才去做接头去除。结果有一批数据的接头含量高达3%，但Q30居然没报警。为什么？因为接头区域的质量分本身就不低，FastQC的Q30统计不会因为序列内容异常而报警。所以，Q30和Adapter Content是两个独立的验证维度，不能互相替代。
好了，我们来看一个实际的质量控制流程概览。下面这张图展示了我日常项目中数据处理阶段与QC的关键节点：
这个流程看起来步骤不少，但实际操作起来很快。我通常把前两步（文件验证 + FastQC）放到集群的并行任务中，因为每个样本是独立计算的——互不依赖，理论上可以用N个线程同时处理N个样本。
测序质量控制在生信Pipeline的落地中就是个护栏——你不能指望次次拿到完美数据，但可以确保不合格数据不会悄无声息地流进下游分析。我记得在某个项目中，正是这个简单的FastQC + MultiQC组合，帮我揪出了一个经常被忽视的问题：样本标签的测序引物交叉污染，最终通过QC的overrepresented sequences模块准确锁定。
认真对待QC，它会保护你整个项目的数据质量。