📚 课程目录第4章集成验证

参考基因组与注释：版本、格式、索引都要验

做生信 pipeline，第一步就是拉参考基因组和注释文件。这一步要是出问题，后面的分析全是废的。我见过太多人随便找个版本就跑，结果比对率低、变异假阳性高，最后查半天发现是参考基因组版本不匹配。

今天咱们就把这三件事聊透：参考基因组版本怎么选、GTF/GFF 注释文件的格式怎么验、基因组索引完整性怎么查。

参考基因组版本选择：别小看这个决定

选参考基因组，不是简单的「GRCh38 vs hg19」这种选择题。版本差异会影响多少结果？嗯，我算过，不同 patch 版本之间，拼接区域、基因号都可能变。你想想看，如果 RNA-seq 定量用的注释来自 GRCh38.p13，而你比对的参考是 GRCh38.p10，可能几千个转录本对不上。

我的个人习惯：先确认数据来源和下游工具要求。比如 1000 Genomes 项目释放的 VCF 用的是 hg19（其实是 GRCh37），那你要么也用 GRCh37，要么做 liftover。别自作聪明用 GRCh38，否则变异坐标全偏。

下面这张表，是我平时做版本决策的 checklist：

考量维度	推荐做法	我的避坑经验
数据来源	优先与原始研究保持同一版本	我曾遇到 TCGA 和 GTEx 用了不同 patch，跨项目整合时花了三天做 liftover，后悔没提前对齐
工具兼容性	查工具文档中列出的推荐版本	有些 aligner 对较新 patch 支持不好，比如 BWA-MEM2 在某个版本对 GRCh38.p14 的 ALT 区域有 bug
注释质量	选择有高覆盖率注释的版本	GRCh38.p13 的 GENCODE v40 注释质量远高于早期的 patch 版
可重复性	固定版本号，不要用「latest」	我常用的做法：在 pipeline 配置里写死 GRCh38.p13，并记录 MD5

小技巧：用 wget 下载参考基因组时，顺便拉 MD5 校验文件。以后谁问你「你用的哪一版」，直接把 MD5 抛过去，谁也赖不掉。

注释文件格式验证：GTF vs GFF，一个符号都别错

GTF 和 GFF，长得像，但结构差异不小。我在项目中见过最坑的一次：有人把 GFF3 文件当成 GTF 喂给 featureCounts，结果数出来的 reads 全成了 0。为什么？因为 GFF 的字段顺序和 GTF 不一样，第 9 列之后 GFF 用「Tag=Value」而 GTF 用「key "value"」，featureCounts 解析失败就全部跳过。

先明确两种格式的核心差异：

特性	GTF (Gene Transfer Format)	GFF3 (General Feature Format)
第9列格式	`key "value"` 双引号值，分号分隔	`Tag=value` 等号值，分号分隔
必需字段	`gene_id`, `transcript_id`	`ID`, `Parent`（可选但有最佳实践）
多值属性	用空格分隔的多个值（不推荐）	用逗号分隔列表：`Dbxref=PMID:123,NCBI:456`
常见工具偏好	STAR, featureCounts (默认)	GATK Funcotator, bedtools 新版本

验证格式，我一般三步走：

看文件头：第一行 ##gff-version 3 就是 GFF3，没有 ## 头但有 #! 可能是 GTF。
数列数：用 awk 检查每行是不是 9 列。
解析第9列：确认 tag 命名规范。

# 快速检查列数
awk '$0 !~ /^#/ {if(NF != 9) print "行 " NR " 有 " NF " 列"}' your_file.gtf

# 检查第9列是否包含双引号（GTF特征）
grep -v '^#' your_file.gtf | head -5 | sed 's/\t/\n/g' | tail -1 | grep -q '"' && echo "看起来是GTF" || echo "不是GTF"

# GFF3验证：等号分隔
grep -v '^#' your_file.gff3 | awk -F'\t' '$9 ~ /=/' | head -1 | cut -f9

我的硬性要求：在 pipeline 里加一个格式校验 step。用 gtftk 或 gffread 都行。

# 用 gffread 验证 GTF 完整性
gffread your_file.gtf -T -o /dev/null 2>gffread.log
if [ -s gffread.log ]; then
  echo "GTF 文件有问题，请检查日志"
  cat gffread.log
fi

注意：有些工具（如 cellranger）要求 GTF 必须包含 transcript_id 和 gene_id，缺一个直接报错。我遇到过用 NCBI 的 GFF 没有 transcript_id，手动 awk 补了一天才填上。

基因组索引的完整性检查：比对前的最后一道关卡

你下载了参考基因组 fasta，然后用 bwa index 跑了一晚上，第二天跑比对，发现报错「index file not found」。怎么回事？啊，可能是磁盘写满了，索引只生成了一半。

我个人习惯：每次比对前，写一个自动检查脚本，确认所有索引文件都存在、大小不为0、且文件头符合预期。

不同比对工具生成的索引文件都不一样，我列个常用表：

比对工具	索引文件列表	检查要点
BWA (mem, aln)	`.bwt, .pac, .ann, .amb, .sa`	五个文件缺一不可，且 `.bwt` 文件的头部字符串包含版本号
Bowtie2	`.bt2 或 .bt2l (for large)`，共4-6个	文件名前缀必须一致，可以用 `bowtie2-inspect` 验证
STAR	`SA, SAindex, chrLength.txt, chrStart.txt, genomeParameters.txt` 等约10个文件	用 `STAR --genomeDir` 时，看 `genomeParameters.txt` 是否存在、大小是否合理
HISAT2	`.ht2` 系列，通常8个文件	检查文件列表是否存在，以及 `hisat2-build` 日志最后一行是否成功

我常用这个脚本来快速验证索引完整性：

#!/bin/bash
# check_index_integrity.sh
REF=$1
TOOL=$2

case $TOOL in
  bwa)
    SUFFIXES=("bwt" "pac" "ann" "amb" "sa")
    for s in "${SUFFIXES[@]}"; do
      f="${REF}.${s}"
      if [ ! -f "$f" ] || [ ! -s "$f" ]; then
        echo "ERROR: $f 缺失或空文件"
        exit 1
      fi
    done
    echo "BWA 索引完整"
    ;;
  star)
    if [ ! -f "${REF}/genomeParameters.txt" ] || [ ! -s "${REF}/genomeParameters.txt" ]; then
      echo "ERROR: STAR 索引不完整（没有 genomeParameters.txt）"
      exit 1
    fi
    echo "STAR 索引完整"
    ;;
  hisat2)
    ls "${REF}".*.ht2 &> /dev/null || { echo "ERROR: HISAT2 索引缺失"; exit 1; }
    echo "HISAT2 索引完整"
    ;;
  *)
    echo "未知工具"; exit 1
    ;;
esac

一个小细节：BWA 的 .bwt 文件有魔数字符串，可以用 xxd 看前几个字节是否包含 "BOWTIE"。我曾在集群上发现 NFS 挂载导致文件头损坏，但文件大小正常，用 xxd 一眼就抓出来了。

最后，关于索引的完整性问题，其实更重要的是保证索引与参考基因组版本严格对应。我曾经犯过一个错：手快用 GRCh38.p10 的参考做了索引，但下载的 GTF 是 GRCh38.p13 的，比对时直接忽略了注释不匹配的区域，导致几十个基因没捡到。后来我就学乖了，把版本信息写进文件描述行的注释里：

# 在参考 fasta 的头部注释中加入版本
>1 GRCh38.p13 Primary Assembly chr1
ACGTACGT...

然后索引脚本里自动提取这个第一行的版本号，跟你 pipeline 配置文件里的版本比对。不一致直接报错退出。

简单来说，验证过程就是三部曲：

版本：对齐参考、注释、工具三者版本
格式：检查列数、分隔符、必需属性
索引：确认文件齐全、大小非零、头部信息一致

每一步都不能省，省了就是在给自己埋坑。

至此，参考基因组和注释文件的验证环节就全部铺开了。每一步都有坑，每一步都能用自动化脚本封住。下次你建新 pipeline 的时候，不妨把这几步写进去——你会感谢自己的。