BLAST数据库构建：从NCBI下载到自定义索引

构建BLAST数据库，说白了就是给比对工具准备好一套“参考词典”。你搜索的序列就是词，数据库就是词典。没有好词典，再厉害的搜索引擎也白搭。

我记得刚做基因组注释那会儿，直接拿NCBI的nr库跑，机器跑了三天没停。后来才发现，根本用不着那么大的库。选对数据库，效率能差一个数量级。

NCBI官方数据库的下载策略

NCBI提供了多层级的数据源。选择哪个，取决于你比对的目的。我个人习惯按这个优先级来挑：

下载NCBI库：update_blastdb与wget

NCBI官方推荐用 update_blastdb 脚本。它会自动解压、校验文件完整性。省心。当然也可以直连FTP，不过我建议你优先用脚本。

方式一：update_blastdb（推荐）

# 下载nt库（核苷酸）
perl update_blastdb.pl --decompress nt

# 下载refseq蛋白库
perl update_blastdb.pl --decompress refseq_protein

# 只下载指定物种（比如人类）
perl update_blastdb.pl --decompress human_genome

方式二：直接wget（适合批量或脚本集成）

# NCBI FTP基础路径
BASE="https://ftp.ncbi.nlm.nih.gov/blast/db"

# 下载nr库的分卷文件
for i in `seq 1 60`; do
    wget ${BASE}/nr.${i}.tar.gz
    tar -xzf nr.${i}.tar.gz
    rm nr.${i}.tar.gz
done

构建自定义索引：makeblastdb实战

很多时候，你不需要全库。比如我只关心某个科的植物基因。那就自己建一个。流程三步走：预处理FASTA → 格式化 → 校验。

第一步：FASTA预处理

# 去除短序列（<50bp），避免索引污染
awk '/^>/ {if(seq) print seq; print; seq=""} !/^>/ {seq=seq $0} END {if(seq) print seq}' raw.fa \
  | seqkit seq -m 50 -o clean.fa

# 统一标识符格式（NCBI要求不含空格和特殊符号）
sed -E 's/[^A-Za-z0-9_.|>]/_/g' clean.fa > target.fa

第二步：运行makeblastdb

# 构建核苷酸数据库
makeblastdb -in target.fa \
            -dbtype nucl \
            -title "MyPlant_CDS_v1.0" \
            -out myplant_db \
            -hash_index

# 构建蛋白数据库
makeblastdb -in proteins.fa \
            -dbtype prot \
            -title "MyPlant_Proteome_v1.0" \
            -out myplant_prot \
            -hash_index

第三步：验证数据库完整性

# 方法1：blastdbcmd查看数据库信息
blastdbcmd -db myplant_db -info

# 方法2：统计序列数量和长度分布
blastdbcmd -db myplant_db -entry all -outfmt "%l" | sort | uniq -c | head

索引参数调优：什么时候该调整？

默认参数对大部分场景够用。但如果你要频繁做短序列比对（比如100bp以下），或者数据库非常大（超过10G），调参能明显提升性能。

版本管理与一致性检查

基因组注释是长周期项目。不同版本之间可能差别很大。我建议你养成几个习惯：

• 建库时在数据库描述里加构建日期和参数版本
• 用 md5sum 记录原始FASTA的校验值
• 把构建脚本和日志放在和数据库相同的目录下

# 示例：建库脚本附带版本信息
#!/bin/bash
DB_DATE=$(date +%Y%m%d)
DB_VERSION="v1.0"
makeblastdb -in genome.fa \
            -dbtype nucl \
            -title "Rice_Genome_${DB_VERSION}_${DB_DATE}" \
            -out rice_${DB_VERSION}_${DB_DATE} \
            -logfile build_${DB_DATE}.log
md5sum genome.fa > genome.fa.md5

# 后续检查时
blastdbcmd -db rice_${DB_VERSION}_${DB_DATE} -info | grep "Title"