3. 人源基因组参考序列与常用数据库
做BLAST分析，第一步就是搞清楚你要比对的“尺子”是什么。
人源基因组参考序列，就是那把尺子。没有它，你测出来的序列就是一堆乱码。我早些年刚接触生信时，觉得随便找个基因组版本用就行。结果项目做完，文章被拒，理由是对照的参考基因组版本过旧。嗯，这个教训挺深刻的。
3.1 人源基因组参考序列的版本演变
人源基因组参考序列，说白了就是科学家们拼接出来的一条“标准答案”。但这条答案一直在更新。目前主流版本有两个：
- GRCh38/hg38：2013年发布，目前最常用。
- GRCh37/hg19：2009年发布，仍在大量旧文献和工具中使用。
- T2T-CHM13：2022年发布的完整版，首次补齐了所有着丝粒和端粒区域。
我个人习惯，新项目一律用GRCh38。只有涉及老数据复现或特定工具兼容性时，才会用hg19。
3.2 常用数据库分类与特点
你想想看，做BLAST找注释，总得知道去哪儿查吧。以下是我常用且推荐的几个核心数据库：
| 数据库 | 数据内容 | 典型用途 | 我推荐的访问方式 |
|---|---|---|---|
| NCBI RefSeq | 人工审核的非冗余参考序列（mRNA、蛋白） | 基因功能注释、保守性分析 | NCBI FTP / BLAST+ 本地库 |
| Ensembl | 自动化注释+部分审校，包含变异信息和基因树 | 多物种比较基因组学、变异影响预测 | Ensembl REST API / 基因组浏览器 |
| UCSC Genome Browser | 整合多种公共数据的可视化平台 | 查看基因组区间、保守性轨道、调节元件 | UCSC Table Browser / 命令行工具 |
| GENCODE | 高精度人源和鼠源基因注释（ENCODE项目官方） | RNA-seq定量、可变剪接分析 | GENCODE FTP |
| dbSNP / ClinVar | 单核苷酸多态性和临床相关变异 | 致病性评估、群体频率检查 | NCBI dbSNP / ClinVar FTP |
我在做群体遗传分析时，最常用的是RefSeq和GENCODE。前者序列准，后者注释全。
3.3 如何获取与构建本地BLAST数据库
线上BLAST虽然方便，但一旦数据量大或者需要对私密数据做分析，还是得本地跑。下面是我最常用的下载和建库流程：
以GRCh38参考基因组和RefSeq mRNA序列为例：
# 1. 下载参考基因组（FASTA格式）
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.26_GRCh38/GCF_000001405.26_GRCh38_genomic.fna.gz
gunzip GCF_000001405.26_GRCh38_genomic.fna.gz
# 2. 使用makeblastdb建索引
makeblastdb -in GCF_000001405.26_GRCh38_genomic.fna \
-dbtype nucl \
-parse_seqids \
-title "GRCh38_genomic" \
-out GRCh38_genomic
# 3. 下载RefSeq mRNA序列（FASTA）
wget https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/mRNA_Prot/human.1.rna.fna.gz
gunzip human.1.rna.fna.gz
# 4. 同样建库
makeblastdb -in human.1.rna.fna \
-dbtype nucl \
-parse_seqids \
-title "human_refseq_mRNA" \
-out human_refseq_mRNA
-parse_seqids参数。这样后续用blastdbcmd按Accession号提取单条序列时会快很多。我刚开始没加这个参数，结果每次提取序列都得全库扫描，一个操作等半小时。
验库是否成功，用这个命令：
blastdbcmd -db GRCh38_genomic -info
输出应该类似：
Database: GRCh38_genomic
Number of sequences: 2398
Number of letters: 3099756071
Number of sequences with undefined ids: 0
...
3.4 不同用途下的数据库选择策略
做BLAST之前，先问自己一句：我到底想找什么？
- 找新基因或同源基因：用NCBI nr/nt库或RefSeq。
- 做RNA-seq定量：用GENCODE或Ensembl的转录本注释。
- 分析变异位点：搭配dbSNP和ClinVar。
- 查调控元件：用UCSC的ENCODE轨道。
- 做比较基因组学：Ensembl的Compara数据库最好用。
我曾在项目中遇过一个情况：用RefSeq做BLAST找出来的基因，在RNA-seq定量时却对不上。后来排查发现，RefSeq的某些转录本版本比GENCODE旧，剪接模式不一致。所以从那以后，做转录组相关分析我就只用GENCODE了。
3.5 本章知识体系结构
下面的SVG总结了人源基因组参考序列与常用数据库间的逻辑关系，以及BLAST应用时如何选择数据源：
3.6 避坑指南与个人经验总结
grch38.fa，避免路径中有空格或非ASCII字符。
还有一点要说：别只看Accession号。有一次我发现BLAST结果里明明有某个基因，但用基因名却搜不到。后来才搞清楚，dbSNP数据库里同一个位点可能对应多个rs号，不同版本之间还有合并或拆分的情况。所以做变异注释时，最好再拿ClinVar做一次交叉验证。
最后，给个实际操作建议：
# 用blastdbcmd检查数据库是否包含期望的序列
blastdbcmd -db human_refseq_mRNA -entry NM_001101.5 -range 1-100
# 输出应该能看到ACTB基因的前100bp
如果输出为空，那就要检查数据库版本是否包含这个Accession号了。