第一章 GFF3 文件概述与格式基础

好，咱们正式开始。GFF3 这个东西，说白了就是生物学家用来描述基因组上各种“特征”的一种标准格式。基因在哪？外显子怎么拼？这个基因有几种剪接体？这些小信息，全靠 GFF3 来记录。

我最初接触它的时候，第一次看那个文件，简直头皮发麻——全是 tab 分隔的字段，还经常混着注释和 FASTA。但后来用得多了就发现，只要搞清楚那 9 列到底在说什么，这格式其实非常简洁。

一句话总结： GFF3 是“序列特征注释”的格子本，每一行代表一个特征（feature），由 9 个 tab 分隔的列组成。

1.1 核心格式规则

GFF3 文件是纯文本，每行一个特征。注释行以 # 开头。真正的数据行必须包含 9 列，列之间用 Tab 分隔，不能用空格。我把这 9 列的定义列在下表，你对照着看会清晰很多。

列号 列名 说明 个人踩坑点
1 seqid 参考序列 ID（染色体、scaffold 等） 不能有大写以外的字符？其实大小写都可以，但要跟参考基因组 fasta 的 ID 完全一致。我遇到过因为多了一个空格导致后续标注全部错位的案例
2 source 注释来源（软件、人工、数据库） 这里常用“.”代表未知。但很多下游工具会过滤掉 source 为空的条目，所以最好写上工具名称，比如 BRAKER2
3 type 特征类型（gene, mRNA, CDS, exon 等） 大小写敏感！SO 术语库中 exon 是小写，CDS 是大写。我曾经把 CDS 写成 cds，结果 IGV 直接不显示
4 start 特征起始位置（1-based） 一定是正整数。如果是负链特征，start 是左边界？不，start 永远取较小的坐标值（5' 端）
5 end 特征终止位置（包含） 必须 ≥ start。区间长度 = end - start + 1。有一次我写的脚本算错了，导致一个外显子长度为 0，直接被 GFF3 验证器报错
6 score 可信度分数（常用 0–1000） 如果工具没给分，写“.”。注意：很多质控工具会把这个字段当作锚点，若它是小数一定要保留原样，别取整
7 strand 链方向：+ 正链，- 负链，. 缺失写“.”，但建议都标上。如果正负链混着写，后续做剪接体统计时会被误解
8 phase CDS 相位（0, 1, 2, .） 只对 CDS 特征有意义。我曾见过有人把外显子的 phase 也填了数字，导致翻译框架计算全乱套
9 attributes 属性键值对（key=value; 分号分隔） 这里坑最多。ID 必须是唯一的，Parent 必须指向已存在的 ID。我习惯用一个 check_gff3_attributes.py 的小脚本提前跑一遍
一个小技巧： 写代码解析 GFF3 时，永远用 split('\t') 分割行，而不是按照空格。因为第九列的属性里也可能有空格，用 split() 默认按空白字符会搞乱列计数。

1.2 层级结构与 ID 规则

GFF3 最牛的地方，就是用 ID 和 Parent 把特征串成树。比如：

chr1	BRAKER2	gene	1000	9000	.	+	.	ID=gene001;Name=ABC1
chr1	BRAKER2	mRNA	1000	8500	.	+	.	ID=transcript_g001.1;Parent=gene001
chr1	BRAKER2	exon	1000	1500	.	+	.	ID=exon_g001.1.1;Parent=transcript_g001.1
chr1	BRAKER2	CDS	1200	1500	.	+	0	ID=cds_g001.1.1;Parent=transcript_g001.1

看到没？mRNA 的 Parent 指向了 gene001，exon 和 CDS 的 Parent 指向了那个 mRNA。这样，一个基因下面可以挂多个转录本，每个转录本下面再挂外显子和 CDS。这种层级结构也是质控的检查重点——任何孤儿特征（有 Parent 但 Parent 不存在）都是不合格的。

下面这张图，我画出了最典型的 GFF3 层级结构，你一看就明白特征之间的父子关系了。

GFF3 特征层级关系（典型真核基因模型） 染色体 / scaffold (seqid) gene (ID=gene001) mRNA-1 (Parent=gene001) mRNA-2 (Parent=gene001) mRNA-3 (Parent=gene001) exon-1 CDS-1 ⚠ 必须保证所有 Parent 引用都能在文件中找到 对应的 ID，否则报错
我曾经犯过的错： 在写批量转换脚本时，把所有特征的 Parent 都写成了 gene001，但每个 mRNA 的 ID 是自动生成的（比如 gene001.t1），结果外显子和 CDS 全部指向了基因，而不是 mRNA。这样的文件到了 IGV 里，基因结构显示就是一条直线，没有外显子间隔。后来我花了一整天用 gff3_validate 才找出问题——所以，一定要验证 GFF3 的层级完整性。

1.3 第九列 attributes 的常见陷阱

第九列是坑最多的地方。我总结几个常见的：

  1. ID 重复： 同一个文件里，ID 必须全局唯一。如果你有两个不同的基因都叫 gene001，后面的会覆盖前面的。检查时可以用 grep -oP '(?<=ID=)[^;]+' file.gff3 | sort | uniq -d 看看。
  2. 值中包含特殊字符： 分号 ; 是属性分隔符，等号 = 是键值分隔符。如果基因名字本身就带分号？别慌，用 URL 编码 %3B 代替。我见过一个物种的基因名叫 ABC1;DEF2，直接写进去导致解析器以为有两个属性。
  3. 属性名拼写： 官方推荐的小写属性名包括 ID, Parent, Name, Alias, Note, Dbxref, Ontology_term。大小写敏感，Name 和 name 是两个不同的属性。很多工具只认 Name。
  4. 缺少必须属性： 按照规范，每个特征必须有 ID（除非是假想特征如 region）。但有些工具生成时漏掉，这种文件在大部分下游流程中直接报废。
我自己的避坑心得： 每次处理新 GFF3 文件前，先用 gff3-validator 跑一遍（NCBI 提供的小工具），它会提示缺失 Parent、ID 重复、不合理相位等问题。如果文件太大，至少用 head -1000 取样验证，别等到分析到一半才报错。

1.4 格式验证小工具的使用

这里给一个我常用的命令行检查组合（假设你安装了 gff3-validator）：

# 基本语法验证
gff3-validator my_annotation.gff3

# 检查是否存在孤儿特征（有 Parent 但该 Parent 不存在）
awk -F'\t' '{if($9 ~ /Parent=/) print}' my_annotation.gff3 | \
  grep -oP '(?<=Parent=)[^;]+' | \
  sort -u > /tmp/parents.txt
awk -F'\t' '{if($9 ~ /ID=/) print}' my_annotation.gff3 | \
  grep -oP '(?<=ID=)[^;]+' | \
  sort -u > /tmp/ids.txt
comm -23 /tmp/parents.txt /tmp/ids.txt   # 输出孤儿特征对应的 Parent

# 统计特征类型分布（快速质控）
awk -F'\t' '{print $3}' my_annotation.gff3 | sort | uniq -c | sort -rn

嗯，这些都是我早期踩坑之后慢慢攒下来的脚本。你想想看，如果一个 50MB 的 GFF3 文件，读到一半报“Parent 不匹配”，你得花多少时间定位？所以，每次上游交付文件，我都先跑一遍这个组合拳。

← 上一章 📖 返回目录 下一章 →