4. 数据质量评估(FastQC)

测序数据下来第一件事是什么?不是急着跑比对,也不是算变异。

先看看质量。我见过不少新手,直接拿原始FASTQ怼进流程,跑了一半报错,回头才发现是接头污染严重。白白浪费半天时间。

所以质量评估这步,我建议每个人都养成习惯:拿到FASTQ,先跑一遍FastQC。你想想看,就几分钟的事,能帮你躲开后面一堆坑。

4.1 为什么非做不可?

FastQC 干的事很简单——给你一份“测序报告”。

它会告诉你:

  • 每条read的质量分布
  • 碱基组成是否均匀
  • 有没有接头残留
  • GC含量是否异常
  • 有没有过量的重复序列
  • N的比例高不高

说白了,它就像体检单上的血常规。单项异常可能没事,但如果好几个指标同时飘红,就得警惕了。

核心原则:先评估,后处理。不要跳过这一步直接去做质控(trimming)。你没有基准线,怎么知道切没切干净?

4.2 快速上手:跑一次FastQC

安装很简单。如果用的 conda,一条命令搞定:

conda install -c bioconda fastqc

运行更简单。一个 FASTQ 文件就这样:

fastqc sample_R1.fastq.gz

一次跑多个文件:

fastqc sample_R1.fastq.gz sample_R2.fastq.gz

或者用通配符处理整个目录:

fastqc *.fastq.gz -o ./fastqc_results/

输出是两个文件:一个 .html 报告和一个 .zip 压缩包。我一般只看 .html,双击就能在浏览器打开。

提速技巧:如果你的数据量很大,比如几十个样本,可以用多线程:

fastqc -t 8 *.fastq.gz

-t 参数指定线程数。不过别开太多,IO 反而会打架。我之前在 32 核服务器上试过 -t 16,结果还不如 -t 8 快。

4.3 报告逐项解读:别被“红叉”吓住

FastQC 报告分十来个模块,每个用绿勾、黄叹号、红叉标识“通过”、“警告”、“失败”。

我第一次看到满屏红叉时,心里咯噔一下。后来发现,有些红叉其实是“正常现象”。比如 RNA-seq 数据在“Per base sequence content”上几乎必挂——因为随机引物带来的偏差。嗯,这里要注意区分“真问题”和“假警报”。

4.3.1 最重要的几个模块

模块名称 看什么 典型问题
Per base sequence quality 每个碱基位置的Q值箱线图 尾部Q值快速下降 → 需要切除低质量碱基
Per sequence quality scores 所有read的平均Q值分布 出现双峰 → 可能混入低质量batch
Per base sequence content 每个位置的A/T/G/C比例 前几个碱基偏移 → 可能是引物偏倚或接头残留
GC content 整个文库GC含量的分布 偏离理论曲线 → 污染或文库扩增偏好
Adapter Content 已知接头序列的累积比例 >5% 就亮黄灯 → 接头污染,必须切除

我曾经在一个WGS项目里,看到Adapter Content从第20bp开始飙升到80%。原始数据没跑FastQC直接去比对,结果比对率只有60%。后来花了一整天重新处理。别学我,先看报告。

4.3.2 那些“挂掉也没事”的指标

  • Per base N content:如果N比例整体很低(<5%),即使某段位置跳红,也往往是随机噪音。不用怕。
  • Sequence Duplication Levels:RNA-seq 中高表达基因会产生大量重复reads,这是生物学信号,不是技术问题。但WGS中出现>50%的重复率,就要怀疑PCR duplication。
  • Overrepresented sequences:偶尔几个高丰度序列可能是rRNA残留或常见接头。如果全是未知序列且比例很高,再查。

怎么判断是真的有问题?我的经验是:看趋势,不看单个点。比如质量箱线图,如果尾部掉到Q20以下,而且整条reads都偏低,那就得切。如果只是最后的几个碱基拉低,有时可以容忍——取决于后续分析对末端质量的敏感度。

4.4 应对策略:挂了的模块怎么办?

下面是我常用的“红框处理清单”:

  1. Per base quality 挂 → 用 Trimmomatic / fastp 做滑动窗口切除,通常设置质量阈值Q20,窗口大小4bp。
  2. Adapter Content 挂 → 先用 FastQC 识别接头类型(它会列出最可能的接头),再用 cutadapt 或 fastp 自动识别并切除。
  3. GC content 异常 → 先排除污染。如果GC曲线出现两个峰,很可能是混入了其他物种的DNA。跑一下Kraken2 或 BLAST。
  4. Sequence Duplication 过高 → WGS数据可以用Picard MarkDuplicates标记重复。RNA-seq 建议保留(除非做CNV分析)。

举个实际命令的例子。我用 fastp 做一步切除质量 + 去接头:

fastp -i raw_R1.fastq.gz -I raw_R2.fastq.gz \
      -o clean_R1.fastq.gz -O clean_R2.fastq.gz \
      --qualified_quality_phred 20 \
      --unqualified_percent_limit 40 \
      --cut_front --cut_front_window_size 1 --cut_front_mean_quality 20 \
      --detect_adapter_for_pe

我的习惯:处理完之后,再用 FastQC 跑一次 clean 数据,确认“红叉”变成了“绿勾”。这叫质控闭环,别省这一步。

4.5 一张图看懂FastQC核心流程

下面这张图概括了从拿到FASTQ到决定处理策略的完整逻辑。

FastQC 质量评估决策流程 FASTQ输入 运行FastQC 查看HTML报告 所有绿勾? 直接进入下游分析 针对红叉模块进行处理 再次FastQC验证 ⚠️ 注意: 并非所有“红叉”都意味着数据不可用。RNA-seq 的碱基含量偏倚、高深度测序的重复率等,需要结合实验设计判断。

4.6 批量处理与自动化

当样本量上来之后,一个个打开HTML报告看肯定不行。我通常会写一个简单的汇总脚本。

FastQC 的输出是 .zip 文件,里面包含 .txt 格式的统计数据。我们可以用 awk 或 Python 提取关键指标。

举个例子,提取所有样本的“基础质量”通过状态:

for zip in fastqc_results/*.zip; do
    unzip -q -o "$zip" -d /tmp/fastqc_tmp
    sample=$(basename "$zip" .zip)
    pass=$(grep "Per base sequence quality" /tmp/fastqc_tmp/*/fastqc_data.txt | awk '{print $NF}')
    echo "$sample $pass"
    rm -rf /tmp/fastqc_tmp/*
done

输出大概这样:

sample01 PASS
sample02 WARN
sample03 FAIL
sample04 PASS

看到 FAIL 的样本就优先处理。不用每个报告都看一遍。

进阶用法:我用 MultiQC 把所有 FastQC 报告合成一个网页,能一眼扫完所有样本的质量概览。后面章节会详细介绍 MultiQC,这里先记住它的名字就好。

4.7 避坑清单:我踩过的几个坑

  • 只看平均值不看分布:有的样本整体平均Q30看起来很漂亮,但50%的reads尾端质量掉到Q10。平均掩盖了问题。要盯着箱线图看。
  • 忽略读长信息:如果测序仪出了问题,可能出现某些read长度异常。FastQC的“Per base sequence quality”如果突然中断,可能是读长不一致。我曾经遇到过一个样本,fastp报告长度正常,但FastQC显示一半的reads只有50bp——原因是双端测序的第二端提前终止。还好FastQC及时发现了。
  • 接头切除后忘了重新评估:切完接头后的数据,再跑一次FastQC。有时候第一次接头识别不准,或者切过了头导致质量下降。我吃过一次亏:切接头时参数设得太激进,把碱基含量切偏了,后续突变检测出现大量假阳性。重跑一次FastQC才发现问题。

一句话总结:FastQC 是你对测序数据的第一道防线。花10分钟看报告,能省下后面10小时调试流程的时间。别跳过它。