第四讲:BLAST 数据库介绍与构建(NCBI nt/nr 等)

为什么要搞明白数据库?

做 BLAST 本质上就一句话:拿你的序列去跟一个已知的数据库比对。说白了,数据库就是你的「参考字典」。字典好不好,直接决定你比对出来的结果有没有价值。我刚入行那会儿,觉得 BLAST 数据库就是个索引文件,随便下个 nt 就能用。结果项目跑了两天,发现数据库版本错了,整批分析结果全废……嗯,从那以后,我每次构建数据库前都会先做三件事:确认版本、校验完整性、检查索引状态。

NCBI 提供两类核心数据库:nt 和 nr。你想想看,nt 是所有核酸序列的集合,nr 是非冗余蛋白库。它们的目标完全不同,使用场景也各有讲究。下面我把常见的几类数据库列出来,方便你做选择。

数据库名称全称数据内容常用场景
ntNucleotide collectionGenBank/EMBL/DDBJ 所有核酸序列(含冗余)核酸比对、物种鉴定、引物设计验证
nrNon-redundant protein非冗余蛋白序列(去冗余后)蛋白功能注释、同源搜索、结构域识别
refseq_rnaRefSeq RNARefSeq 精选 RNA 序列(高置信度)转录组注释、基因结构验证
refseq_proteinRefSeq ProteinRefSeq 精选蛋白序列高质量蛋白功能预测
swissprotUniProt Swiss-Prot人工校验的蛋白序列高置信度功能注释
我个人习惯:做大规模基因组注释时,先拿 nr 做一轮初筛,再用 Swiss-Prot 做精细验证。这样既快又准。

数据库构建:从下载到索引

你从 NCBI 下载的是 FASTA 文件,BLAST 不能直接用它搜索。你需要先通过 `makeblastdb` 构建成二进制索引。为什么?我举个例子:文本文件里找一条序列要逐行扫描,索引文件则像字典的拼音目录,直接定位到目标位置。

下面是下载和构建 nr 数据库的完整流程。注意:nr 今天已经超过 100 GB,网速慢的话别全量下载。我建议先用 `update_blastdb.pl` 脚本增量更新,或者只下你需要的子集。

1. 使用 update_blastdb.pl 自动下载

# 下载 nt 数据库(核酸)
perl update_blastdb.pl --decompress nt

# 下载 nr 数据库(蛋白)
perl update_blastdb.pl --decompress nr

# 只下载特定部分(比如只下载环境样本库)
perl update_blastdb.pl --decompress env_nt

这个脚本会自动解压 `.tar.gz` 文件,最终得到 `nr.00.tar.gz`、`nr.01.tar.gz` …… 以及里面的 FASTA 文件和索引文件。它还会自动跳过已下载且完整的文件,你只管跑就行。

我曾经遇到过:直接在集群上用 `wget` 下载全部 tar 包,结果网络中断导致解压失败。后来我改用脚本加 `--decompress` 参数,它会校验每个包的 MD5,省心多了。

2. 手动构建索引

如果因为网络限制,你只能拿到 FASTA 文件(比如从第三方镜本站下载),那就需要手动构建。

# 核酸库索引
makeblastdb -in nt.fa -dbtype nucl -out nt_db

# 蛋白库索引
makeblastdb -in nr.fa -dbtype prot -out nr_db

# 给索引添加标题(方便后续识别)
makeblastdb -in my_genome.fa -dbtype nucl -title "MyGenome v1.0" -out my_genome

参数说明:

  • -dbtype:必须跟原始序列类型一致,核酸选 nucl,蛋白选 prot。搞反了会报错 "BLAST engine error"。
  • -out:输出前缀,后续 BLAST 命令的 -db 参数就填这个前缀。
  • -title:非必须,但建议加上,尤其在管理多个版本时。

3. 索引文件结构

构建成功后你会看到一堆同名前缀的文件,比如 nr_db.phr、nr_db.pin、nr_db.psq…… 它们各自存了不同的索引信息:

后缀内容
.pbd / .nbdBLAST 数据库头文件(版本、统计信息)
.phr / .nhr序列头部信息(ID、描述行)
.pin / .nin序列编号到位置映射
.psq / .nsq序列数据本身(二进制压缩)

你自己不需要手动改这些文件。它们全是二进制格式,BLAST 程序运行时自动加载。

核心原则:每次更新数据库后,必须重新构建索引,否则 BLAST 会报错 "Database seems corrupt" 或结果不正确。别问我怎么知道的……

SVG:数据库构建与比对流程

BLAST 数据库构建与使用流程 选择数据库 nt / nr / refseq / swissprot 下载 FASTA update_blastdb.pl / wget 构建索引 makeblastdb BLAST 比对 blastn / blastp / blastx … 结果解析 XML / 制表符 / 自定义格式 提示:每个数据库对应一套独立索引;更新任意数据库后必须重新构建索引

自定义数据库:构建自己的专属库

除了 NCBI 的公共库,你完全可以把自己的序列做成 BLAST 数据库。比如我手头有一个新测序的菌株,要跟已有的 20 个菌株比较,那就把它们的 CDS 或蛋白序列抽出来,做成一个私有库。

构建自定义库的步骤

  1. 将所有目标 FASTA 合并成一个文件(注意序列 ID 不能重复,否则会报 warning)。
  2. 运行 makeblastdb,指定 -dbtype 和 -out。
  3. 用 blastdbcmd 验证索引是否正常:blastdbcmd -db my_db -info
避坑提醒: 如果你的 FASTA 序列名包含空格或特殊字符,最好用下划线替换掉。否则 BLAST 输出结果里的 query 字段会断掉。
# 示例:构建一个自定义核酸库
cat genome1.fna genome2.fna > combined_cdna.fa

makeblastdb -in combined_cdna.fa -dbtype nucl -title "TwoGenomes_CDNA" -out two_genomes

# 验证
blastdbcmd -db two_genomes -info

输出类似这样:

Database: TwoGenomes_CDNA
    Number of sequences: 8572
    Number of letters: 12,345,678
    Number of letters in sequence data tree: 0
    Number of letters remaining after clustering: 0
    Total number of ambiguities residues: 0
    Number of isoelectric points: 0
    Number of molecular weights: 0
    Number of sequences with undefined pI/MW: 8572
    Date: 2025-03-15

看到 "Number of sequences: 8572" 就说明索引建好了。如果数字明显不对,可能是合并时漏了或者有空的序列头。

常用数据库维护命令

我整理了几条日常用得最多的命令,你们可以收藏起来:

操作命令
查看数据库统计信息blastdbcmd -db my_db -info
列出数据库中所有序列 IDblastdbcmd -db my_db -entry all -outfmt %i
提取指定序列(按 ID)blastdbcmd -db my_db -entry GI:12345 -out seq.fa
检查数据库是否完整blastdb_aliastool -db my_db -check
删除旧数据库直接删除所有同前缀文件即可:rm my_db.*
我常用 blastdbcmd -entry all -outfmt %i 快速提取所有序列 ID,再跟我的目标列表取交集,省掉自己写脚本解析 FASTA 头。

数据库选择策略总结

说到这,你可能会问:到底该用哪个库?我分享一个简单的决策树:

  • 快速筛选物种:用 nt,跑 blastn。因为核酸比对可以直接定种。
  • 功能注释:用 nr 跑 blastp(蛋白 vs 蛋白),同时打开 -evalue 1e-5,覆盖度要超过 70%。
  • 高置信度注释:用 Swiss-Prot 跑 blastp,结果几乎不用再过滤。
  • 跨物种比较:用 refseq 或 nr,别用 nt(不同物种的核酸差异太大)。
  • RNA 注释:优先用 refseq_rna,或者用 nt + 加 -task blastn-short(短序列模式)。
记住一件事:数据库磁盘占用通常比你想象的大 2~3 倍。索引文件加上解压后的 FASTA,几百 GB 稀松平常。我见过有人直接把 nt 下到 /tmp 下,跑了两天发现磁盘写满了……提前规划好存储路径和配额。

好了,数据库这部分算是讲透了。你要做的是动手跑一次:下载 nt 的一个子集(比如只下载 env_nt 或 refseq_rna),构建索引,然后用 blastdbcmd 看看里面有什么。亲自走一遍流程,比看十遍文档都管用。