第二讲：宏基因组测序数据处理基础——FASTQ格式、质量控制、接头去除

咱们上一讲已经搭好了Docker环境，现在得聊聊原始数据长什么样。说白了，你从测序仪拿到的就是一堆FASTQ文件，里面藏着微生物群落的“密码”。但直接分析这些原始数据？我劝你别这么干——坑太多了，我当年第一次跑Kraken2时就被垃圾数据坑过。

2.1 FASTQ格式：读懂测序仪的“四行诗”

每个测序read在FASTQ文件里占四行。结构其实很简单：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
GTGGCTGCGGTGGCTCGTGCTTGGTGATCTGACCTGCGTGCCGTCGTCGCCGCCGAGCCG
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

• 第1行：以@开头，是read的标识符（仪器ID、流道坐标、barcode等）。
• 第2行：碱基序列，A/T/C/G/N。
• 第3行：以+开头，通常只写个+，老版本会重复第1行。
• 第4行：质量评分，每个字符对应一个Phred质量分数。

那这个质量评分咋看的？我用个表给你说清楚：

2.2 质量控制：别让垃圾数据毁了你的分类结果

宏基因组数据里常见的“脏东西”有三类：

• 低质量碱基：测序仪读不准的区域，Q值低于20（错误率≥1%）的碱基。
• 重复序列：PCR扩增产生的duplicate reads，会扭曲丰度估计。
• 接头序列：建库时加的adaptor，没切除干净就混在read里。

我常用FastQC做初步检查，再用Trimmomatic或cutadapt做精细处理。运行一下你就知道数据有多糙：

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o fastqc_raw/

打开HTML报告，主要看这几个指标：

• Per base sequence quality——箱线图在Q30以下超过30%碱基的区域要小心。
• Adapter Content——如果出现接头富集，必须切除。
• Overrepresented sequences——占比超0.1%的未知序列可能是污染。

2.3 接头去除：精准切除“人工尾巴”

接头序列通常是已知的，比如Illumina的通用接头：

Universal Adapter (P5-P7): AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC
Index Adapter: 依index不同而异，常见GATCGGAAGAGCACACGTCTGAACTCCAGTCAC

用cutadapt切除，我一般这样写：

cutadapt \
  -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
  -A AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
  -o trimmed_R1.fastq.gz \
  -p trimmed_R2.fastq.gz \
  -m 50 \
  sample_R1.fastq.gz sample_R2.fastq.gz

参数解释：

跑完切接头的流程，你还可以再做一次质量控制检查，看看接头残留率是否降到1%以下。嗯，这个步骤我一般在容器里用一条pipeline搞定：

# Docker容器内一键质控+切接头
docker run --rm -v /data:/data biocontainers/trimmomatic:0.39 \
  trimmomatic PE \
  /data/sample_R1.fastq.gz /data/sample_R2.fastq.gz \
  /data/trimmed_R1.fastq.gz /data/unpaired_R1.fastq.gz \
  /data/trimmed_R2.fastq.gz /data/unpaired_R2.fastq.gz \
  ILLUMINACLIP:/opt/conda/share/trimmomatic-0.39-1/adapters/TruSeq3-PE.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50

2.4 本章核心流程概览

下面我把整个原始数据处理流程画成一张图，你看一遍就能记住：

2.5 常见问题与避坑

• Q：双端测序的R1和R2质量不一致怎么办？
  A：我曾经遇到过R1质量很好但R2一片糊的情况。我当时检查了文库制备记录，发现是R2测序通道堵塞。我的处理方式是：单独对R2设置更严格的滑动窗口参数（比如SLIDINGWINDOW:4:20），甚至丢弃R2中大量低质量片段。
• Q：切完接头后read变短了，还能用于Kraken2吗？
  A：可以，但要注意Kraken2的最小k-mer长度是默认31bp，如果你的read长度小于31bp，k-mer就完全匹配不上。所以我建议设置MINLEN:50或更高（如75）。
• Q：质控报告里过表达序列占比0.5%算异常吗？
  A：对宏基因组来说，0.5%不算高，但如果该序列比对到人类基因组或phiX对照，就需要排除。你可以用--contaminants文件把污染序列过滤掉。