第三章 数据获取：公共数据库检索与FASTQ格式详解
公共数据库三巨头：GEO、SRA、ENA
做RNA-seq差异表达分析，第一步永远是找数据。公共数据库说白了就是全球实验室的"数据仓库"。我个人习惯把这三个库看作一套体系：GEO（Gene Expression Omnibus）是NCBI下的数据门户，SRA（Sequence Read Archive）是它的原始数据仓库，ENA（European Nucleotide Archive）是欧洲的兄弟版本。
| 数据库 | 全称 | 主要用途 | 我推荐的操作 |
|---|---|---|---|
| GEO | Gene Expression Omnibus | 查阅实验设计、样本信息、处理后的表达矩阵 | 先在这里搜实验，确认数据可用 |
| SRA | Sequence Read Archive | 存储原始FASTQ/FASTQ.gz，需要prefetch下载 | 从GEO的SRA链接跳转，拿SRR编号 |
| ENA | European Nucleotide Archive | 和SRA互为镜像，有时下载更快 | 碰上SRA太慢时，切到ENA试试 |
为什么会这样？三个库其实是同一个项目的不同入口。我记得刚入行时有个坑：在GEO上看到了漂亮的统计图，兴冲冲下载表达矩阵，结果发现原始测序数据不在里面。嗯，这里要注意——GEO只能告诉你"这个实验存在"和"最终结果长了什么样"，原始测序文件一定要去SRA或者ENA拿。
GEO检索实操
我在项目中遇到过最头疼的情况：搜索关键词太宽，溢出几千条数据；关键词太窄，又漏掉重要实验。这里分享一个我反复验证过的工作流。
- 第一步：明确生物问题和条件——比如"人乳腺癌细胞系MCF-7，DOX处理 vs 对照"，比搜"breast cancer"精准十倍
- 第二步：在NCBI GEO主页使用Advanced Search Builder，组合条件：
(breast cancer[MH]) AND (MCF-7[Text Word]) AND (Expression profiling by high throughput sequencing[DataSet Type]) - 第三步：筛选结果——看Series个数，一般5-50个是合理范围。点进去后重点看"Samples"表格，确认是否有对照和处理、生物学重复数目。
我个人习惯在GEO上先读那个"Summary"区块，如果实验设计写得清楚，我会继续下去；如果描述模棱两可，直接跳过，后面踩坑的概率极大。
SRA下载实战
拿到SRR编号后，下载工具最主流的是SRA Toolkit里的prefetch和fasterq-dump。我的标准流程如下。
# 下载SRA文件（.sra格式，压缩态）
prefetch SRR1234567
# 转换成FASTQ（默认生成两个文件：SRR1234567_1.fastq, SRR1234567_2.fastq）
fasterq-dump --split-files SRR1234567
# 想更精细控制：指定输出目录、线程数
fasterq-dump --split-files --outdir ./fastq/ --threads 4 SRR1234567你想想看，为什么我还要用prefetch再fasterq-dump？直接fastq-dump不就行了吗？其实prefetch下载的是SRA专有格式，可以断点续传，网络断开后不用从头开始。我经历过在服务器上下载20G数据到95%断网，没这个功能真的会崩溃。
prefetch了20个样本，结果因为写入磁盘太猛把inode用光了。现在我的习惯：每次同时下载不超过5个样本，并且用--max-size 50G限制单个文件大小。还有，千万别在/tmp里直接跑——空间不够的话，你会发现辛辛苦苦下到一半的文件就没了。
如果SRA官方的面镜像速度太慢，你可以切到ENA：ENA的下载地址格式通常是：ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR123/SRR1234567/SRR1234567_1.fastq.gz。我本人习惯优先从ENA下载，实测快不少，尤其是在欧洲的节点上。
- 关键参数：
--split-files：双端测序时把read1和read2分到两个文件，必加--threads：根据CPU核心数设置，我一般用4或8--outdir：指定输出目录，防止散落在当前目录
- 一次性下载多个样本：
# 准备SRR_Acc_List.txt，每行一个SRR号 cat SRR_Acc_List.txt | xargs -n 1 -P 3 prefetch
FASTQ格式：一字一句都不能错
FASTQ是RNA-seq的金标准格式。说白了，它就是给每个序列配上一个质量分数。一下是一个FASTQ记录的完整拆解。
@SRR1234567.1 HWI-ST1234:1:1101:1234:1234 length=100
ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII| 行号 | 标识符 | 内容 | 踩坑高频条目 |
|---|---|---|---|
| 第1行 | @ | 序列ID + 可选仪器/坐标信息 | ID中的空格可能导致下游工具解析错误，千万别自己添加空格 |
| 第2行 | 序列 | 碱基序列（ATCGN） | 如果N太多说明测序失败，这种read我会在质控后直接淘汰 |
| 第3行 | + | 通常只写一个"+"，有时会重复ID | 我遇到过某些工具要求第三行必须与第一行完全一致，否则报错，用seqtk可以标准化 |
| 第4行 | 质量 | 与第2行等长的ASCII字符，每个字符对应一个碱基的质量值 | 质量编码模式（Sanger/Illumina）不匹配是最常见的错误 |
为什么质量值用ASCII字符而不是数字？你想想看，如果直接用数字写，一行100个数值，那文件就膨了。ASCII字符相当于一个"查表式压缩"：每个字符的ASCII码值减去33（Sanger编码）就是质量分数。
# 用Python快速检查质量编码是不是Phred+33
from Bio import SeqIO
for record in SeqIO.parse("sample.fastq", "fastq"):
qual = record.letter_annotations["phred_quality"]
# 如果最小值小于0，那肯定是Phred+64（旧版本），需要转换
if min(qual) < 0:
print("检测到Phred+64编码，建议转换")
break
elif max(qual) > 40:
print("可能存在偏移")
break
else:
print("Sanger编码，一切正常")
break # 只用第一个read判断就够了实战流程：搜索→下载→格式确认
下面这张SVG图完整展示了从数据库检索到拿到合格FASTQ的路线。我一直觉得，看得见流程比闷头敲命令更重要。
看到这张图，你应该能理解：数据获取不是一条单线任务，而是一个需要反复验证的循环。我个人习惯在每一步都做一次"完整性确认"：GEO检索完确认实验设计，下载完确认文件不损坏，格式转换完确认质量编码正确。一次完整的RNA-seq分析，数据获取这部分如果做扎实了，后面80%的错误都可以避免。
SRR1234567_Control_Rep1，永远不会错。
格式细节：为什么只有4行却总出问题
FASTQ看起来简单，但我在教程里见过太多人踩坑了。
- 换行符问题：Windows的CRLF（\r\n）和Linux的LF（\n）混用会导致某些工具无法解析。我个人习惯用
dos2unix对所有FASTQ文件统一处理。 - 降碱基符号：除了A、T、C、G，还可能出现IUPAC模糊字符如R、Y、M、K。大部分比对工具都能处理，但如果你用自定义bin计数脚本，一定要考虑N和模糊字符。
- 质量值标度：Sanger（Phred+33）质量范围0-40；Illumina 1.8+（也是Phred+33，但上限可能到41或更）。我见过不少工具用
Phred+33当作固定假设，如果文件是Phred+64，所有质量值都会偏移33，导致错误低估。
# 用seqtk检查质量编码是否Sanger（Phred+33）
seqtk test sample.fastq | head -n 4这里面其实有一个隐藏的技巧：seqtk test命令会输出一个简化质量分布，你能看到第一个read的质量值ASCII范围。如果最大字符接近'I'（73 ASCII），可能是Phred+64；如果最大字符小于'K'，大概率是Phred+33。我在项目启动前10秒用这招，避免等跑完比对才发现编码不对。
fastqc输出的"Per base sequence quality"图。如果Quality scores over range显示的值完全不对，比如全部<10或>40，十有八九就是编码模式选错了。这种情况我会写一个简单脚本重新映射。
说到底，数据获取和格式理解是RNA-seq分析中最"苦"但又最关键的基线工作。很多人急于做差异表达，在这个阶段马虎了事，结果后面的分析全都建立在沙地上。我的建议：花40%时间在数据获取和质控上，这比花80%时间在调试下游分析省心多了。