4. 读取FASTQ文件：用Python把数据“抓”出来

拿到一个FASTQ文件，第一件事是什么？当然是把它读进去。我在刚入行时，以为读个文本文件很简单，等真上了几百万条 reads 的数据，才发现 效率 和 正确性 缺一不可。这节我就带你用 Biopython 和 pysam 这两个库，把 FASTQ 里的序列、质量和 ID 整整齐齐地取出来，再做点基本统计。

4.1 为什么不用 Python 自带 open？

你想想看，FASTQ 每 4 行为一条记录，如果用 open() 按行读，还得自己处理换行符、分隔符、质量字符串里的“@”和“+”混淆问题。我很久以前写过一次“纯手工”解析，结果一不小心把序列行里的“@”当成了标题行…… 从那以后我就老老实实用现成的库了。

我曾经踩过的坑：千万别自己用正则去解析 FASTQ 的 4 行结构。FASTQ 的质量行可以出现任何 ASCII 字符，包括 '@'。Biopython 和 pysam 已经帮你处理了这些边界情况，直接用就好。

4.2 用 Biopython 读取 FASTQ

Biopython 的 SeqIO.parse 是解析 FASTQ 最常用的入口。它返回一个迭代器，每条记录是一个 SeqRecord 对象，包含 id、seq 和 letter_annotations['phred_quality']。

from Bio import SeqIO
import os

fastq_path = "sample.fastq"

records = []
for record in SeqIO.parse(fastq_path, "fastq"):
    records.append(record)
    # 只取前 1000 条测试
    if len(records) >= 1000:
        break

print(f"成功读取 {len(records)} 条记录")
print(f"第一条序列长度: {len(records[0].seq)}")
print(f"第一条质量分数: {records[0].letter_annotations['phred_quality'][:10]}")

我个人习惯在遍历时加上 count 计数器，防止文件太大导致内存溢出。SeqIO.parse 返回的是延迟加载的迭代器，不转成列表就不会把整个文件塞进内存。

提示：如果文件是 gzip 压缩的（.fastq.gz），Biopython 会自动识别，无需额外解压。传入 SeqIO.parse(gzip.open("sample.fastq.gz", "rt"), "fastq") 就可以。

4.3 用 pysam 读取 FASTQ（更贴近生产环境）

pysam 本来是处理 BAM 文件的，但它也提供了 FastxFile 类来读取 FASTQ、FASTA。如果你的下游流程需要和 BAM 文件交互，pysam 的接口和它自带的 C 扩展能带来稍快一点的速度。

import pysam

with pysam.FastxFile("sample.fastq") as fh:
    for i, entry in enumerate(fh):
        if i >= 1000:
            break
        # entry 对象包含 name, sequence, quality
        print(entry.name, len(entry.sequence), entry.quality[:10])
        # 注意：pysam 的 quality 是字符串，不是数值列表

对比一下：Biopython 返回的是 SeqRecord，质量是数值列表；pysam 返回的是简易命名元组，质量是原始字符串。我通常在分析流程的中间环节用 pysam，因为它不需要额外转换就能直接对接 samtools 命令。

特性 Biopython SeqIO pysam FastxFile
质量格式 list of int (Phred 分数) str (原始 ASCII)
内存占用 迭代器式，低（小心转列表） 迭代器式，低
压缩支持 需手动用 gzip.open 自动识别 .gz 文件
并行解析 可选 multiprocessing 有内置的线程池支持
额外功能 丰富的序列操作 与 BAM/CRAM 索引交互方便

4.4 基本统计：给你的数据“摸个底”

读完数据，下一步就是算算序列长度分布、碱基质量分布。我常用的是在读取过程中一边读一边统计，不额外多扫一遍文件。

from Bio import SeqIO
import numpy as np

lens = []
qual_sums = []
qual_counts = []

for record in SeqIO.parse("input.fastq", "fastq"):
    lens.append(len(record.seq))
    quals = record.letter_annotations['phred_quality']
    qual_sums.append(sum(quals))
    qual_counts.append(len(quals))

lens = np.array(lens)
qual_avg = np.array(qual_sums) / np.array(qual_counts)

print(f"序列长度: min={lens.min()}, max={lens.max()}, mean={lens.mean():.1f}")
print(f"平均质量: min={qual_avg.min():.1f}, max={qual_avg.max():.1f}, mean={qual_avg.mean():.1f}")
我建议：在统计前先检查文件是否为空，或者有没有超过合理长度（比如 1e6 bp）的异常条目。我在实际项目中遇到过一条 200 kbp 的“污染”序列，如果不做过滤，后面的比对效率会非常低。

4.5 归纳总结：解析 FASTQ 的核心流程

不管用哪个库，套路都一样：打开文件 → 迭代记录 → 提取 ID、序列、质量 → 聚合统计或逐条处理。下面这张图帮你把关系理清楚。

FASTQ 解析与基本统计流程 FASTQ 文件 ( .fastq / .fastq.gz ) Biopython SeqIO pysam FastxFile 迭代记录：name, sequence, quality 统计序列长度分布 统计平均/各位置质量 输出：min / max / mean / 直方图等
一句话归纳：Biopython 适合做序列分析的前期探索（质量转数字、序列操作丰富），pysam 适合嵌入生产流程（速度稍快、与 BAM 交互友好）。选哪个？看你的下游需要什么。如果只是读 FASTQ 做统计，两个随便挑，我偏爱 Biopython 因为调试时能看到更多细节。

4.6 避坑指南

  • 文件名大小写：FASTQ 可能叫 .fq 或 .fastq，Biopython 的 format 参数必须填 "fastq"，不能写 "fq"。
  • 质量编码偏移：绝大多数现代 FASTQ 使用 Sanger 偏移（Phred+33）。如果遇到老数据（如 Solexa），Biopython 会自动处理，pysam 需要手动设置 quality_scale="solexa"。
  • 迭代不重置：迭代器是一次性的，如果先统计了长度，需要重新读取文件才能再算质量。我建议将统计逻辑合并到同一个循环里，或者把重要字段存为结构化数据（如 pandas DataFrame）。
小技巧：当你手上有几十个 FASTQ 文件时，用 glob.glob("*.fastq") 配合循环，再把统计结果汇总到一张 CSV 表里。我在一个 200 个样本的项目里就是这样干的，后续质控报告直接由 CSV 生成。

好了，读取并统计 FASTQ 的基本功就这些。下节课我们看看怎么把这些原始数据跟元数据（样本名、条件、批次）挂上钩——那才是真正让数据“活起来”的一步。

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