数据预处理：别让脏数据混进你的分析流程

大家好。上一章我们聊了批次效应的全局概念，今天直接上手干活。数据预处理这块，说白了就是在正式分析前，先把原始数据的“底细”摸清楚。我见过太多人跳过这一步，结果跑完流程才发现问题出在源头——那真是白忙活好几天。

预处理主要看三件事：原始数据格式、质量分数编码、接头污染。这三样搞不定，后续所有分析都是空中楼阁。咱们一个一个来拆。

数据预处理核心流程 原始数据（FASTQ/FASTA） 格式检查 质量分数编码确认 接头污染识别 必检项 三项通过才能往下走 常见故障 格式错乱 / 编码混淆 / 接头残留

原始数据格式检查：FASTQ 不是随便写的

你拿到的原始数据，大概率是 FASTQ 格式。这是一种四行一组的文本格式。我当初刚接触测序数据时，第一次看到 FASTQ 文件真有点懵——乱七八糟的字母和符号。

标准 FASTQ 记录长这样：

@EAS139:136:FC706VJ:2:2104:15343:197393
GTGCTTAGCCTCGAGTAGAGAGGTGATAATCATGACCTCATC
+
EEEEEEEEEEEEEEEEEEEEEEEEEEDDDDDDDDDDDDDDDD

四行不能乱：

  • 第1行：以 @ 开头，是序列标识符（包含仪器、流道、坐标等信息）。
  • 第2行：原始碱基序列（A/T/C/G/N）。
  • 第3行：以 + 开头，通常可空着或重复标识符。
  • 第4行：质量分数字符串，长度必须与第2行碱基序列长度一致。

格式问题常见于：

  • 行尾换行符混用（Windows \r\n vs Linux \n）。
  • 质量行与碱基行长度不一致——这会导致解析器崩溃。
  • 标识符行用 > 开头（那是 FASTA 格式，不是 FASTQ）。

快速验证方法很简单，用 Linux 命令行看一眼前几条记录：

zcat sample_R1.fastq.gz | head -n 20

你看，头四条记录是不是规整的四行结构？我习惯再用 awk 检查一下所有序列长度是否一致：

zcat sample_R1.fastq.gz | awk 'NR%4==2 {len=length($0)} NR%4==0 {if(length($0)!=len) print "长度异常，行号：" NR-1}'

嗯，这个命令能快速揪出行数不对齐的坏记录。有次我在项目中遇到的批次，200 万条序列里有 3 条质量行少了一个碱基——就是这个命令抓出来的。

质量分数编码：Phred+33 还是 Phred+64？

质量分数是测序仪给每个碱基打的分，告诉你这个碱基叫对的概率有多大。但坑就坑在编码方式上。

目前主流是 Phred+33 编码，质量值 Q 从 0 到 41 或更高，映射到 ASCII 码 33-126。老式 Solexa 或早期 Illumina 机器用过 Phred+64，ASCII 从 64 开始。

编码方式 ASCII 范围 质量值范围 常见平台
Phred+33 33-126 0-93 Illumina 1.8+ 几乎所有现代数据
Phred+64 64-126 0-62 Illumina 1.3-1.7, Solexa

判断方法很简单。我个人习惯用 FastQC 的 Basic Statistics 模块，它会在页面顶部直接告诉你编码类型。如果没条件跑 FastQC，也可以用命令行手动判断。

举个例子，你已经知道 Phred+33 的质量字符最小是 !（ASCII 33）。那可以看看文件里出现的最小 ASCII 值是什么：

zcat sample.fastq.gz | awk 'NR%4==0 {for(i=1;i<=length($0);i++) {ascii=ord(substr($0,i,1)); if(ascii<min || min=="") min=ascii; if(ascii>max || max=="") max=ascii}} END {print "最小ASCII:", min, "最大ASCII:", max}'

不过这个命令要有个 ord() 函数，需要提前写个小脚本。更稳妥的办法是看着质量字符的分布：如果大部分是 E、F、G（ASCII 69-71），那很可能是 Phred+64。

我个人的经验： 现在 99% 的公共数据库数据都是 Phred+33。只有老项目或某些特定的 Illumina 旧版本才会碰到 +64。但我仍然坚持每个批次都确认一遍——因为一旦编码判断错，整个质量过滤阶段就全错了。

接头污染识别：别让接头序列浑水摸鱼

接头污染，就是测序时没切干净的 adapter 序列。它们会像影子一样粘在 reads 的末端，严重影响比对和组装质量。

常见接头类型：

接头名称 序列（5'→3'） 常用于
Illumina TruSeq Universal Adapter AGATCGGAAGAGC 单端或双端通用
Illumina Small RNA 3' Adapter TGGAATTCTCGG 小 RNA 测序
Nextera Transposase Sequence CTGTCTCTTATACACATCT Nextera 文库构建

我用 FastQC 的 Overrepresented sequences 模块来快速筛查。如果某个序列出现频率异常高，而且能匹配到公共接头库，那就是污染没跑了。

手动验证也可以用 grep 搜索：

zcat sample_R1.fastq.gz | grep -c "AGATCGGAAGAGC"

但要注意，接头可能被部分截断，所以我更喜欢用 cutadapt 先做一轮检测，它会输出接头匹配的详细统计。

cutadapt --detect-adapter-for-pe -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o /dev/null -p /dev/null sample_R1.fastq.gz sample_R2.fastq.gz 2>&1 | grep "Total reads"

我曾经踩过一个坑： 有批数据看起来接头污染很严重，但仔细一看，其实是文库构建时引入了一段基因特异性序列，被 FastQC 误判为接头。类似的假阳性情况不少，尤其在高重复区域。所以自动检测的结果一定要手动核对几条序列，别直接全删了。

实战操作：三步完成预处理核查

下面是我自己项目中用的标准流程，你可以直接拿来用。

  1. 格式检查：先用 zcat | head 看结构，再用 awk 判长度一致性。
  2. 质量编码确认：运行 FastQC，看 Basic Statistics 页面的 Encoding 字段。如果是 Phred+64，需要转换。
  3. 接头识别：跑 FastQC 的 Overrepresented sequences，配合 cutadapt --detect-adapter-for-pe 做验证。

遇到异常怎么办？别慌，先把问题归类：

  • 格式问题：用 seqtk 或自定义脚本修复行结构。
  • 编码错误：用 seqtk seq -Q64 -V sample.fastq.gz > fixed.fastq 强制转换为 Phred+33。
  • 接头残留：用 cutadapt 或 Trimmomatic 切除。

行了，预处理这关就这么过。记住一个原则：宁可多花半小时检查原始数据，也别在分析跑完后才发现源头有问题。这是我在一次次翻车中学到的教训。

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