数据获取与准备：下载公开 RNA-seq 数据与理解质量编码
做 RNA-seq 分析，第一步就是拿到原始测序数据。市面上主流的数据都存放在公共数据库里，最常用的是 NCBI 的 Sequence Read Archive (SRA)，还有 ENA、DDBJ 这些镜像站。我个人习惯优先用 SRA，因为它覆盖面最广，工具链也最成熟。
但你千万别一上来就去下载几十个样本的全量数据。我在早期项目里干过这事——下载了 30 个 GB 的 FASTQ，结果发现编码格式不对，重新下载又花了一整天。后来我学乖了，先拿一个小数据跑通流程，再批量操作。
从 SRA 下载 FASTQ
SRA 里的数据以 .sra 格式存储，我们需要用 NCBI 提供的 sratoolkit 工具转成 FASTQ。最核心的两个命令是 prefetch 和 fastq-dump。我的习惯是先用 prefetch 把 .sra 文件拉到本地，再用 fastq-dump 拆成 FASTQ。
安装 sratoolkit 很简单，直接从 NCBI 官网下载预编译包解压就行。我一般放在 /opt/bioinfo/sratoolkit 目录下，然后加到 PATH 里。
vdb-config --interactive 配置一下下载路径和代理，尤其是你在国内，需要设置合适的镜像，不然速度会慢到怀疑人生。
假设我们要下载一个公共数据集 SRR390728（来自人类 RNA-seq 测试数据），命令如下：
# 先下载 .sra 文件
prefetch SRR390728
# 再转成 FASTQ（默认单端）
fastq-dump SRR390728
# 如果是双端测序，加上 --split-files
fastq-dump --split-files SRR390728
# 或者直接一步到位，用 fasterq-dump（速度更快）
fasterq-dump SRR390728fasterq-dump 是 sratoolkit 里的新工具，多线程处理，速度能快 5-10 倍。我强烈推荐用它。但要注意，fasterq-dump 默认不压缩输出，你最好自己加上 --outfile 控制输出名，或者用管道配合 gzip。
# 输出压缩的 FASTQ
fasterq-dump --outdir ./fastq --temp ./tmp SRR390728
gzip ./fastq/SRR390728.fastqfastq-dump 不加 --split-files 只会输出所有 reads 混在一起，后续比对全部乱掉。所以，双端数据一定要显式指定 --split-files，养成肌肉记忆。
下载完 FASTQ 后，先看一眼文件大小和行数，确认没有损坏。我常用 wc -l 检查行数，FASTQ 每条记录 4 行，所以行数应是 4 的倍数。
wc -l SRR390728.fastq理解测序质量编码体系
拿到 FASTQ 文件之后，很多人会直接扔进比对流程。但我觉得先花 5 分钟看一眼质量编码体系，能帮你避免后面一大堆麻烦。
FASTQ 中每条记录的最后一行是质量字符串，每个字符对应一个 Phred 质量分数。Phred 分数 Q 与碱基错误率 P 的关系是：Q = -10 log10(P)。比如 Q30 表示错误概率 0.001，准确率 99.9%。
但问题在于，不同时期、不同测序平台使用的 ASCII 编码偏移量不一样。早期 Solexa/Illumina 1.0 用 Phred+64，后来 Illumina 1.8+ 统一成 Phred+33。如果你用 Phred+33 的解析器去读 Phred+64 的文件，质量值会全部偏高，造成假阳性。
| 编码体系 | ASCII 偏移量 | 质量范围 | 常见平台 |
|---|---|---|---|
| Phred+33 | 33 | 0–41 (Sanger) | Illumina 1.8+, Sanger, PacBio |
| Phred+64 | 64 | 0–40 (Illumina 1.0–1.7) | 早期 Illumina, Solexa |
怎么快速判断？看 FASTQ 质量字符串中出现的 ASCII 字符范围。Phred+33 的最小字符是 ! (ASCII 33)，最大常见是 J (ASCII 74)；Phred+64 的最小字符是 @ (ASCII 64)，常见是 h 或以上。我写了个小脚本，从文件前 1000 行里提取质量字符串，统计字符的 ASCII 编码范围，一眼就能看出编码体系。
import sys
from collections import Counter
def detect_phred_offset(fastq_file, n=1000):
"""从 FASTQ 前 n 条记录中检测 Phred 偏移量"""
with open(fastq_file, 'r') as f:
lines = []
for i, line in enumerate(f):
if (i % 4) == 3: # 质量行
lines.append(line.strip())
if len(lines) >= n:
break
all_chars = ''.join(lines)
if not all_chars:
raise ValueError("文件为空或记录数不足")
min_ord = min(ord(c) for c in all_chars)
max_ord = max(ord(c) for c in all_chars)
if min_ord >= 64:
return 64
elif min_ord >= 33:
return 33
else:
raise ValueError(f"无法识别编码，min_ord={min_ord}")
print(detect_phred_offset("SRR390728.fastq"))数据准备流程一览
到这里，我们的数据获取与准备工作可以总结成下面这张图。我把关键步骤和每个节点的注意事项都标出来了，方便你快速对照。
你看，图上标注了两个并行检查：质量编码检测和完整性检查。这两个步骤我建议在下载完 FASTQ 后立即执行，哪怕只是简单跑个脚本。记住，原始数据一错，后面所有分析都是垃圾进垃圾出。
一条龙示例：从 SRA 到可用的 FASTQ
最后我把整个流程串成一个 bash 脚本，你拿过去稍微改改 SRR 号就能直接用。这也是我项目里的标准 pipeline 第一步。
#!/bin/bash
# 批量下载并检查 FASTQ
SRA_LIST="SRR390728 SRR390729 SRR390730"
OUTDIR="./fastq"
TMPDIR="./tmp"
TOOLKIT_PATH="/opt/bioinfo/sratoolkit"
export PATH="$TOOLKIT_PATH/bin:$PATH"
mkdir -p "$OUTDIR" "$TMPDIR"
for acc in $SRA_LIST; do
echo "===== Processing $acc ====="
# 1. 下载
prefetch "$acc" --output-directory "$TMPDIR"
# 2. 转为 FASTQ（双端）
fasterq-dump --outdir "$OUTDIR" --temp "$TMPDIR" --split-files "$acc"
# 3. 压缩
gzip -f "$OUTDIR/${acc}_1.fastq"
gzip -f "$OUTDIR/${acc}_2.fastq"
# 4. 质量编码检测（Python 脚本）
python3 detect_phred.py "$OUTDIR/${acc}_1.fastq.gz"
# 5. 行数检查
zcat "$OUTDIR/${acc}_1.fastq.gz" | wc -l | awk '$1%4!=0{print "ERROR: not multiple of 4 in",acc}'
done
echo "All done."这个脚本里用到了前面写的 detect_phred.py。你只要把 SRA_LIST 换成自己的样本列表，就能自动完成从下载到准备的过程。
嗯，数据获取与准备就讲到这里。你拿到手的 FASTQ 文件，现在应该是干净的、知道编码体系的。下一件事，就是对这些原始数据进行质量评估——我会用 FastQC 和 MultiQC 来生成质控图表，到第四章我们再详细展开。