3. 数据准备

群体联合基因分型，数据质量是命根子。我做这个领域这些年，踩过的坑十有八九都跟数据质量有关。所以本章咱们老老实实把原始数据格式搞清楚，再把质量控制这关做扎实。

3.1 FASTQ：测序仪的「原生语言」

你拿到手里的下机数据，绝大多数是 FASTQ 格式。它长得像这样：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

一个记录就四行：

  • 第1行：以 @ 开头，是序列标识符，里面藏了仪器、流道、坐标、Index 等信息。我个人习惯用 zcat 先看一眼标识符规律，能快速判断数据来源。
  • 第2行：碱基序列，全部大写，A/C/G/T/N。注意 N 代表未识别碱基，比例高了说明测序质量差。
  • 第3行：以 + 开头，通常跟第1行相同标识，但也可以为空。嗯，这里不重要。
  • 第4行：质量分数，每个字符对应一个碱基的 Q 值。字符转 Q 值用 Phred 质量公式：Q = -10 * log10(P)，其中 P 是错误概率。
小贴士：质量字符的编码有历史遗留问题。Sanger 标准（Phred+33）是目前绝大多数平台使用的，字符范围是 !（Q=0）到 ~（Q=93）。你常看到的 Illumina 1.8+ 就是这种。我曾经碰见过老数据用 Solexa+64 编码，FastQC 会自动检测，但跑 Trimmomatic 时要小心，参数要写对。

3.2 FASTQ 的质量分数，你得会看

Q30 这个概念你肯定听过。Q30 意味着错误率千分之一。我一般拿 Q20 当作质检的「及格线」，Q30 是「良好线」。下表是常见 Q 值与错误率的对应：

Q值 错误概率 碱基正确率 我通常的解读
Q10 10% 90% 极差，基本不能用
Q20 1% 99% 及格，看项目要求
Q30 0.1% 99.9% 良好，群体分析推荐水平
Q40 0.01% 99.99% 优秀，但我很少见整条 read 都是 Q40

3.3 FastQC：质量报告的「体检单」

拿到 FASTQ 的第一步，跑 FastQC。说白了它就是给每条 FASTQ 文件做一次全面体检，输出一个 HTML 报告。安装很简单：

# 直接用 conda 装，跨平台无烦恼
conda install -c bioconda fastqc

# 或者下载二进制，解压即用
wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.12.1.zip
unzip fastqc_v0.12.1.zip
chmod +x FastQC/fastqc

运行命令：

# 对单端数据
fastqc sample1_R1.fastq.gz -o ./fastqc_reports/

# 双端数据一起跑
fastqc sample1_R1.fastq.gz sample1_R2.fastq.gz -o ./fastqc_reports/

# 加 -t 参数指定线程数
fastqc -t 4 *.fastq.gz -o ./fastqc_reports/

报告里我最关注的几个模块：

  • Per base sequence quality：看每个碱基位置的 Q 值箱线图。我一般要求整个 read 前 90% 的区域中位数不低于 Q30，否则直接报警。
  • Per sequence quality scores：观察平均 Q 值的分布。峰值左移说明整体质量差。
  • Adapter Content：接头污染情况。测序读长越长，越容易跑到接头区域。我在一个 RNA-seq 项目里就见过 30% 的 read 尾部被接头污染，Trimmomatic 直接切掉了 20bp。
  • Overrepresented sequences：看看有没有异常的富集序列。有时候是引物二聚体或 rRNA 残留。
注意： FastQC 的「PASS」标记只是参考阈值，并非绝对合格。比如 GC 含量这条，人类基因组是 41% 左右，细菌可能 50%+，项目不同容忍度不同。不要太依赖红绿灯，自己读图更重要。

3.4 Trimmomatic：数据清洗的「瑞士军刀」

Trimmomatic 是我最喜欢的清洗工具。它功能全、速度快、参数灵活。安装同样推荐 conda：

conda install -c bioconda trimmomatic

基础命令太长了，咱们拆开看参数。双端清洗的例子：

java -jar trimmomatic-0.39.jar PE \
  sample1_R1.fastq.gz sample1_R2.fastq.gz \
  sample1_clean_R1.fastq.gz sample1_unpaired_R1.fastq.gz \
  sample1_clean_R2.fastq.gz sample1_unpaired_R2.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:15 \
  MINLEN:36

参数说明（我个人习惯的典型值）：

参数 含义 我推荐的值
ILLUMINACLIP 去除接头，后面分别是：接头文件、允许的最大错配数、比对阈值、最小匹配长度 2:30:10
LEADING 从 5' 端开始切除低于该 Q 值的碱基 3
TRAILING 从 3' 端开始切除低于该 Q 值的碱基 3
SLIDINGWINDOW 滑动窗口切除：窗口大小:平均质量阈值 4:15
MINLEN 切除后 read 最短长度，太短会被丢弃 36

你现在可能会想，为什么用 LEADING:3 而不是 LEADING:10？我解释一下：LEADING 是从头部开始切，每切一个碱基就重新检查下一个，直到遇到质量达标的。3 是 Q3，对应错误率 50%，看似很低，但真实场景中 5' 端质量通常很好，设太高反而切到正常碱基。TRAILING 同理。而 SLIDINGWINDOW:4:15 才是主体——窗口内4个碱基，平均 Q 值低于15就切掉窗口，然后往后滑动。这个组合是我在几十个项目中验证过的，对大多数 Illumina 数据都很稳健。

3.5 清洗后的复核：再跑一次 FastQC

千万别偷懒。清洗完一定要用 FastQC 再看一遍。我见过有人清洗完直接扔到 GATK 里，后来发现 adapter 根本没去掉，白白浪费半天跑流程。两次 FastQC 对比一下：

  • 接头污染曲线应该消失。
  • base quality 箱线图尾部明显上提。
  • reads 总数减少（因为 MINLEN 过滤了短片段），但整体质量提高了。

避坑指南： 我曾经在一个大型群体项目中，用 Trimmomatic 的默认参数去跑了100个样本，结果发现 ILLUMINACLIP 用的接头文件是单端版，双端数据的 adapter 根本没匹配上。后来重新跑了一遍，损失了3天。所以千万检查接头文件是否匹配你的文库构建试剂盒。TruSeq3-PE.fa、TruSeq2-PE.fa、NexteraPE-PE.fa 别搞混。

3.6 实战流程一览

我把整个过程画了个流程图，你跟着思路走就不会漏：

原始 FASTQ 数据质控与清洗流程 原始 FASTQ 文件 FastQC 质控报告 质量OK？ Trimmomatic清洗 再次 FastQC 清洗后FASTQ 可用的干净 FASTQ 清洗后若质量仍不佳，需调整Trimmomatic参数或检查原始数据来源

3.7 双端数据要特别注意：配对关系

群体项目基本都是双端测序。清洗后两个 FASTQ 文件的 read 数量必须一致，且顺序一一对应。Trimmomatic 会自动处理配对关系，它会输出 _forward_paired.fq 和 _reverse_paired.fq 以及 _unpaired_*.fq。但有时候脚本跑飞了，你检查一下行数：

# 统计压缩文件的行数（除以4得到read数）
zcat sample_clean_R1.fastq.gz | echo $((`wc -l` / 4))
zcat sample_clean_R2.fastq.gz | echo $((`wc -l` / 4))
# 如果两个数不等，后面 GATK 会直接报错

我给你的建议是：写一个简单的检查脚本，跑完 Trimmomatic 自动对比 read 数。别问我是怎么知道的——血泪史。

个人经验：当数据特别大（比如几百个样本），我会写个 Snakemake 或 Nextflow 流水线，把 FastQC→Trimmomatic→FastQC 串联起来，统一查看多报告。MultiQC 也是个好东西，能把所有 FastQC 报告汇总成一个页面，一眼就能看出哪些样本有问题。

3.8 本章落脚

数据准备就这些。FASTQ 格式不难理解，但很多人忽视质量分数和接头污染。FastQC 和 Trimmomatic 这对组合拳，打好了后面分析才能站得住。你把本章的指令跑一遍，记住两条：先看报告再动手，洗完后复检一遍。

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