2. 基因组数据获取：从NCBI Genome抓取细菌序列

好，上一章我们理顺了BLAST注释的整体逻辑。现在要动手了——第一步，先把细菌的基因组序列搞到手。

我个人习惯用NCBI Genome数据库，这玩意儿是生物信息领域的“粮仓”。大肠杆菌、沙门氏菌、结核分枝杆菌……不管你要注释哪个菌，序列都得到这里拿。说白了，数据源不对，后面分析全是白搭。

2.1 怎么找到目标基因组

打开NCBI Genome页面，直接搜菌名就行。但有个坑——同名菌株可能有一堆。比如大肠杆菌K-12，光MG1655和W3110就是两个不同的版本。我刚开始做项目时，就错过一次。

具体搜索时，我一般用这个组合查询：

("Escherichia coli"[Organism] OR "E. coli"[All Fields])
AND ("complete genome"[Assembly Level])
AND ("latest"[Version])

为什么要加“complete genome”？因为NCBI里还有scaffold、contig级别的组装，那些碎片化的序列，做BLAST注释时你会哭的。

搜索条件	作用
[Organism]	限定物种，避免名字混淆
complete genome	只拿完整基因组，不要碎片
latest	确保拿到最新版组装

2.2 下载FASTA格式的序列文件

找到目标基因组后，NCBI显示的是一个Assembly页面。你要找的按钮是“Download Assemblies”。点开，选FASTA格式。

注意了：NCBI提供两种FASTA——genomic.fna（整个基因组）和cds_from_genomic.fna（编码序列）。我们做BLAST注释需要前者，因为要扫全基因组水平的信息。

命令行用户可以用curl直接拉，省去手动点击的麻烦：

curl -O https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz

这段命令看着长，其实结构很清晰：

• GCF_000005845.2——NCBI的基因组登录号，每个菌株唯一
• ASM584v2——组装版本名
• _genomic.fna.gz——文件后缀名，FASTA格式加gzip压缩

2.3 检查文件完整性

下载完了别急着解压。先检查MD5校验和——NCBI在每个基因组目录下都放了md5checksums.txt。我遇到过一个坑，直接下载后以为是完整的，结果BLAST跑到一半报错，原来是文件头部被截断了。

# 比对MD5
md5sum GCF_000005845.2_ASM584v2_genomic.fna.gz

# 查看NCBI的MD5列表
curl -s https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/md5checksums.txt \
  | grep genomic.fna

嗯，这里要注意：NCBI的MD5是对原始文件算的，不是对压缩包。所以你得先解压再算，或者用unzip -p管道传入。

2.4 解压与基本查看

确认文件没事后，解压：

gunzip GCF_000005845.2_ASM584v2_genomic.fna.gz

然后用head看一眼文件内容，确保它是标准的FASTA格式：

head -n 5 GCF_000005845.2_ASM584v2_genomic.fna

正常看到的应该是：

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTG
GTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGAC
AGATAAAAATTACAGAGTACACAACATCCATGAAACGCATTAGCACCACCATTACCACCACCATCACCATTACCACAGGT
AACGGTGCGGGCTGACGCGTACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAG

第一行是>开头的描述行，后面是序列。没有乱七八糟的空格或特殊字符，这就对了。

2.5 快速统计基因组基本参数

我个人做完这些步骤后，会顺手统计几个关键指标：

# 统计序列总长、GC含量等
seqkit stats GCF_000005845.2_ASM584v2_genomic.fna

输出类似：

file                           format  type   num_seqs    sum_len  min_len  avg_len  max_len
GCF_000005845.2_ASM584v2_genomic.fna  FASTA   DNA          1  4,641,652  4,641,652  4,641,652  4,641,652

为什么做这步？你想想看，如果下载的序列长度和文献报道的基因组大小差太多，八成是找错了菌株或文件不完整。这时候回头排查，比费半天劲跑BLAST要省事儿得多。

2.6 整理一下全流程

说了这么多，其实核心流程就这几步：

1. 定位——NCBI Genome页面搜菌名，加好查询限制条件
2. 确认——检查Assembly页面，看看是完整基因组还是片段，下载genomic.fna.gz
3. 校验——比对MD5，防止传输损坏
4. 解压——gunzip解压，检查FASTA格式是否正确
5. 验证——用seqkit统计基本参数，和数据表核对

下面这张图，我把整个数据获取的逻辑画在了一起。你看一眼，心里就有谱了。

2.7 拿到数据后，先确认这几点

文件放桌面、解压完、seqkit跑完——你是不是觉得可以关网页了？别急，再问自己三个问题：

• 登录号准确吗？ ——我见过有人把RefSeq的GCF和GenBank的GCA搞混，前者是参考基因组，后者是原始提交，有时候序列不完全一样。
• 线环都标注清楚了吗？ ——细菌基因组常有质粒，如果你只拿了主染色体，质粒上的抗性基因就漏了。
• 文件解压后是纯文本吗？ ——有些软件对换行符敏感，Windows的CRLF可能带来麻烦。我在Linux下统一用dos2unix洗一遍，求个安心。