📚 课程目录第3章 BLAST数据库

BLAST 数据库：格式说明、公共资源与自定义构建

做 BLAST 分析，说白了，你查的不是原始序列文件，而是经过特殊处理的数据库。我刚开始干这行的时候也犯过迷糊——拿着 FASTA 文件就直接往上怼，结果程序报错，半天找不到原因。后来才明白，BLAST 需要的是“索引化”的数据，就跟查字典前得先有拼音索引一样。

这一章，咱们就把 BLAST 数据库的里里外外说清楚。从格式原理，到怎么下公共数据库，再到自己动手建库，一条龙讲明白。

3.1 数据库格式与索引原理

BLAST 用的数据库，不是简单的文本文件，而是经过 makeblastdb 工具处理后的二进制索引文件。为什么要这么折腾？为了快。你想想看，一次比对可能要搜索几亿条序列，挨个读原始文本得等到猴年马月。

构建完成后，你会看到一堆后缀不同的文件。下面这张表，把常见后缀和它们的作用说清楚了：

后缀	文件类型	说明
.phr	头部索引	存储序列头信息的索引
.pin	信息索引	存储数据库元数据，如序列数量、长度
.psq	序列数据	压缩后的序列碱基/氨基酸数据
.pal / .nal	累积索引	用于多卷数据库的累计索引
.pos / .nos	偏移索引	记录每条序列在 .psq 中的位置
.pal / .nal	对数索引	加速高斯比对查找

嗯，这么多文件看着挺吓人，但你只需要记住一点：它们是一整个整体，缺一不可。移动或复制数据库时，必须打包带走所有同前缀的文件。

小贴士：我个人习惯在项目目录下建一个 db/ 文件夹，专门放数据库文件。所有索引文件保持原始命名，不要手动改名。否则 BLAST 会找不到对应的索引。

3.2 常用公共数据库下载

做生物信息学，免不了跟几个巨无霸数据库打交道。NCBI 的 BLAST 数据库家族是最常用的。我来给你盘一盘最核心的几类：

nr（non-redundant）——最通用的蛋白数据库，去冗余后的综合数据库。几乎所有物种的比对都能用它。
nt（nucleotide）——核酸版 nr，覆盖范围同样极广，但数据量比蛋白版大得多。
refseq_protein / refseq_rna——NCBI 人工审核的参考序列数据库，质量高、注释全。我建议做精准分析时优先选这个。
swissprot——UniProt 提供的经过手工注释的蛋白数据库，准确率极高，但覆盖度不如 nr。
env_nr——专门用于环境样本（宏基因组）的蛋白数据库，包含了大量未培养微生物的序列。

公共数据库的下载，推荐用 update_blastdb.pl 这个 Perl 脚本，或者直接用 curl/wget 从 NCBI 的 FTP 拉。我自己更常用脚本，因为它自动处理压缩包和解压，省事儿很多。

# 使用 update_blastdb.pl 下载 nr 数据库
# 先下载脚本
wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/update_blastdb.pl
chmod +x update_blastdb.pl

# 下载 nr 数据库（蛋白）
./update_blastdb.pl --decompress nr

# 下载 refseq_protein 数据库
./update_blastdb.pl --decompress refseq_protein

# 下载 16S_ribosomal_RNA 数据库（核酸）
./update_blastdb.pl --decompress 16S_ribosomal_RNA

注意：nr 数据库压缩包解压后有 50GB 以上，nt 更是超过 100GB。下载之前一定确认磁盘空间足够。我曾经在服务器上下了一半才发现 /home 分区炸了，白白浪费半天流量。

如果你只是想快速验证流程，不需要全量数据，可以试试小一些的专用库，比如 pdbaa（PDB 蛋白）或 landmark（模式物种）。

# 下载小库用于测试
./update_blastdb.pl --decompress pdbaa
./update_blastdb.pl --decompress landmark

3.3 构建自定义数据库

公共数据库虽好，但很多时候不够灵活。你的分析可能只关注某一个物种、某一条信号通路，或者自己测序得到了新的序列。这时候，就需要动手建自己的库。

建库的核心命令就是 makeblastdb。下面我给出两个最常用的场景：

场景一：蛋白序列建库

假设你有一个 my_proteins.fasta，包含了你关注的基因家族序列。

makeblastdb -in my_proteins.fasta \
            -dbtype prot \
            -out my_proteins_db \
            -title "My custom protein database" \
            -parse_seqids

场景二：核酸序列建库

对于基因组或转录组数据，你只需把 -dbtype 改成 nucl。

makeblastdb -in my_genome.fasta \
            -dbtype nucl \
            -out my_genome_db \
            -title "My custom genome database" \
            -parse_seqids

参数解释一下：

-in：输入的 FASTA 文件
-dbtype：指定为 prot（蛋白）或 nucl（核酸）
-out：输出数据库文件的前缀名，后面所有索引文件都以此开头
-title：为数据库设置一个可读的标题，会在 BLAST 报告里显示
-parse_seqids：解析序列 ID，方便后续用 blastdbcmd 按 ID 检索序列

  核心要点：自定义数据库是保证可重复计算的关键一环。把数据库的构建命令写入 Snakefile 或 shell 脚本，配合 CI/CD，就能做到“一键重建、结果可溯”。我自己团队的所有分析流程，都要求必须从原始 FASTA 开始重建数据库，绝不用别人建好的二进制库。

还有一件事要提醒你。FASTA 文件的序列 ID 必须唯一且不包含特殊字符。否则 -parse_seqids 会报错，或者检索时找不到对应序列。避坑指南我写在这里：

避坑指南：我曾经有一个项目，序列 ID 里包含了竖线（|）和空格，结果 makeblastdb 执行成功了，但后续 blastdbcmd 死活提取不出来。检查了好久才发现是 ID 格式的问题。建议序列头格式统一为：>序列标识符|描述信息，并用下划线代替空格。

3.4 数据库验证与信息查看

数据库建好之后，怎么确认它没问题？用 blastdbcmd 这个工具。

# 查看数据库基本信息
blastdbcmd -db my_proteins_db -info

# 查看所有序列 ID 列表
blastdbcmd -db my_proteins_db -entry all -outfmt %i

# 提取某一条序列
blastdbcmd -db my_proteins_db -entry "sequence_ID_123"

-info 的输出会告诉你序列数量、长度分布、字母表类型等。我每次建完库都会跑一遍这个命令，作为日志记录下来。这样后面如果发现比对结果异常，可以回溯到数据库本身排查故障。

说到这，我把本章的内容框架用一张图总结一下：

这张图把数据库的整个生命周期串起来了：从原始 FASTA 出发，经过 makeblastdb 加工成索引文件，然后可以来自公共资源或自己构建，最后供 BLAST 程序使用。理解了这个流程，你就能在可重复计算的框架下，把数据库构建这个环节牢牢控制住。