第4章：参考基因组索引问题——缺失.fai或.dict

运行GATK时，你是不是也碰到过这样的错误：

## ERROR: Could not find index file for reference genome
## ERROR: The genome file must have a .fai index and a .dict dictionary.

别急，这不是你代码写错了。说白了，就是参考基因组缺了两个“身份证”文件。我个人习惯管它们叫索引和字典，少了哪一个，GATK都不干活。

4.1 这两文件到底干啥的？

GATK处理序列比对、变异检测时，得快速定位参考基因组上的位置。.fai是FASTA索引，让程序能随机跳到染色体的某个位置；.dict是序列字典，记录染色体的名字、长度、顺序。你可以想象没有目录的百科——翻死你。

文件	全称	作用	创建工具
.fai	FASTA index	支持随机读取，记录每条序列的偏移量和碱基行宽	samtools faidx
.dict	Sequence dictionary	记录染色体名、长度、MD5校验值，GATK用它做一致性检查	Picard CreateSequenceDictionary 或 gatk CreateSequenceDictionary

我在项目里见过有人只建了.fai就开跑，结果GATK报错说找不到.dict。嗯，这两个都得有。

4.2 快速生成两个索引文件

生成方法其实就两行命令。先给FASTA文件建.fai：

samtools faidx reference.fasta

运行完，同级目录下多出一个reference.fasta.fai。接着建.dict：

gatk CreateSequenceDictionary -R reference.fasta -O reference.dict

或者如果你用的Picard（老版本）：

java -jar picard.jar CreateSequenceDictionary R=reference.fasta O=reference.dict

记得把reference.dict和reference.fasta放在同一个文件夹，名字要匹配。GATK默认会找.fasta.dict或者.fa.dict，所以reference.dict就够了。

4.3 避坑：索引文件必须与FASTA一致

注意！索引文件不是通用的。如果你换了参考基因组版本（比如从hg19换成hg38），或者修改过FASTA文件（比如重命名了染色体），索引必须重新生成。

我曾经遇到一个case：同事从公共数据库下载了一个已经带.fai和.dict的参考基因组，但他用sed改了染色体的>行（去掉chr前缀），结果所有下游分析都报错。排查了两天才发现是索引里的染色体名和实际FASTA不匹配。这坑我替你们踩过了。

4.4 检查索引是否存在、是否过期

你可以在运行GATK之前写个简单的shell检查：

if [ ! -f "${REF}.fai" ]; then
    echo "ERROR: .fai not found for $REF"
    exit 1
fi
if [ ! -f "${REF%.*}.dict" ] && [ ! -f "${REF}.dict" ]; then
    echo "ERROR: .dict not found for $REF"
    exit 1
fi

更严谨一点，检查FASTA文件和.fai的修改时间是否一致。如果FASTA比.fai新，说明源文件变过，得重建。我习惯用find的-newer比较，但这属于进阶操作。初期你只要保证每次换参考基因组后，手动跑一遍索引命令就ok。

4.5 常见错误日志对照

我把常见的几个报错和解决方案列一下，方便你快速对号：

报错片段	原因	解法
Could not find index file for reference	缺少.fai	运行samtools faidx重建
No such file: reference.dict	缺少.dict	运行CreateSequenceDictionary
Index not found for file	.fai文件名不对	检查.fai是否与FASTA同名且同一目录
Sequence dictionary does not match FASTA	.dict内容与FASTA不一致	删除旧的.dict，重新生成

4.6 索引文件关系图

下面这张图帮你理解FASTA、.fai、.dict三者的关系：

从图中可以看到，FASTA是源文件，.fai和.dict是两个衍生文件，缺一不可。GATK启动时会把FASTA载入索引，然后检查.dict里的序列列表，如果不匹配就会报错退出。

4.7 批量处理多个参考基因组

如果你手头有多个参考基因组（比如不同物种或者补丁版本），可以写个循环一次性搞定：

for file in *.fasta; do
    base="${file%.*}"
    echo "Processing $file..."
    samtools faidx "$file"
    gatk CreateSequenceDictionary -R "$file" -O "${base}.dict"
done

注意：如果某个.fasta文件很大（比如人类基因组3G），samtools faidx很快，但CreateSequenceDictionary可能慢一些。我一般在服务器后台跑，用nohup挂起来，不耽误干活。

4.8 总结一下（就两句）

搞丢了.fai或.dict，别慌，两条命令搞定。重点是养成好习惯：每次下载或修改参考基因组，第一时间补全这两个文件。你想想看，要是跑了一天的变异检测，最后发现是索引问题，多亏。