水稻基因组注释整合实战：从零散GFF到结构化注释文件
水稻作为重要的模式生物和粮食作物，其基因组注释质量直接影响着分子生物学研究的效率。然而许多研究者第一次从RAP-DB或RGAP下载注释文件时，往往会遇到一个令人头疼的问题——官方提供的GFF/GTF文件信息分散，缺少完整的基因结构层级。本文将分享一套完整的解决方案，帮助您将零散的注释文件整合为可直接用于下游分析的标准化格式。
1. 水稻基因组注释现状解析
目前主流的水稻基因组注释主要来自两个权威数据库：RAP-DB和RGAP。虽然两者基于相同的日本晴（Nipponbare）参考基因组，但在注释策略和ID命名体系上存在显著差异：
- RAP-DB：采用形如
Os01g0100100的编号系统，与KEGG等代谢通路数据库兼容性好 - RGAP：使用
LOC_Os01g01010格式的标识符，在部分功能预测工具中更常见
这两个项目提供的GFF文件都存在信息碎片化的问题。以RAP-DB为例，其下载页面提供以下文件：
| 文件类型 | 包含信息 | 缺失信息 |
|---|---|---|
| locus.gff | 基因位置 | mRNA、外显子等结构 |
| transcripts.gff | mRNA、UTR、CDS | 基因级注释 |
| transcripts_exon.gff | mRNA、外显子 | UTR、CDS信息 |
这种分散的存储方式导致研究人员无法直接获得完整的基因模型，进而影响后续的变异注释、差异表达分析等工作流程。
2. 数据准备与环境配置
2.1 原始数据获取
首先从RAP-DB官网下载所需文件：
# 下载GFF压缩包 wget https://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_representative_annotation_2021-12-15.tar.gz # 解压获取关键文件 tar -xzvf IRGSP-1.0_representative_annotation_2021-12-15.tar.gz解压后将得到三个关键GFF文件：
locus.gff：包含基因级位置信息transcripts.gff：包含转录本、UTR和CDS注释transcripts_exon.gff：包含转录本和外显子信息
2.2 工具准备
我们需要以下工具链来完成注释整合：
# 安装基础工具 conda create -n rice_annotation python=3.8 conda activate rice_annotation conda install -c bioconda gffutils biopython pandas提示：建议使用Conda管理环境以避免依赖冲突
3. 注释文件整合流程
3.1 染色体命名统一化
RAP-DB使用chr01-chr12的命名方式，而许多分析工具更适应Chr1-Chr12的简写格式。我们先进行染色体名称转换：
import gffutils def standardize_chr_names(input_gff, output_gff): with open(input_gff) as infile, open(output_gff, 'w') as outfile: for line in infile: if line.startswith('chr'): line = line.replace('chr0', 'Chr').replace('chr', 'Chr') outfile.write(line) standardize_chr_names('locus.gff', 'locus_std.gff') standardize_chr_names('transcripts.gff', 'transcripts_std.gff')3.2 构建完整基因模型
关键步骤是将基因级注释与转录本级注释合并：
def merge_annotations(gene_gff, transcript_gff, output_gff): # 创建临时数据库 db = gffutils.create_db( ':memory:', dbs=[gene_gff, transcript_gff], merge_strategy='merge', id_spec=['ID', 'Parent'] ) # 输出完整GFF with open(output_gff, 'w') as f: for gene in db.features_of_type('gene'): f.write(str(gene) + '\n') for mRNA in db.children(gene, featuretype='mRNA'): f.write(str(mRNA) + '\n') for exon in db.children(mRNA, featuretype='exon'): f.write(str(exon) + '\n') for cds in db.children(mRNA, featuretype='CDS'): f.write(str(cds) + '\n') for utr in db.children(mRNA, featuretype=('five_prime_UTR', 'three_prime_UTR')): f.write(str(utr) + '\n') merge_annotations('locus_std.gff', 'transcripts_std.gff', 'rice_complete.gff')这个流程确保了最终GFF文件包含完整的层级结构：
- gene → mRNA → exon/CDS/UTR
- 所有特征具有正确的Parent-Child关系
- 染色体命名格式统一
4. 下游应用适配
4.1 生成TxDb对象
对于R用户，可以将整合后的GFF转换为TxDb对象：
library(GenomicFeatures) # 创建TxDb并保存 txdb <- makeTxDbFromGFF("rice_complete.gff", format="gff3") saveDb(txdb, file="Osativa_MSU_TxDb.sqlite") # 使用示例 genes <- genes(txdb) transcripts <- transcriptsBy(txdb, by="gene")4.2 构建BSgenome包
如需创建可共享的BSgenome数据包：
library(BSgenome) # 准备种子文件 writeLines( c( "BSgenomeName: OsativaMSU", "Provider: RAP-DB", "Species: Oryza sativa", "Genome: MSU7", "Source: https://rapdb.dna.affrc.go.jp", "chromosomes: Chr1,Chr2,Chr3,Chr4,Chr5,Chr6,Chr7,Chr8,Chr9,Chr10,Chr11,Chr12", "circular: NA", "masked: NA" ), "seed.txt" ) # 构建包 forgeBSgenomeDataPkg("seed.txt", seqs_srcdir="path/to/fasta", destdir=".")5. 质量验证与问题排查
为确保整合后的注释文件质量，建议进行以下检查：
结构完整性验证：
db = gffutils.FeatureDB("rice_complete.db") problematic_genes = [] for gene in db.features_of_type('gene'): mRNAs = list(db.children(gene, featuretype='mRNA')) if not mRNAs: problematic_genes.append(gene.id)ID转换表生成（RAP↔RGAP）：
import pandas as pd # 从RAP-DB下载ID对应表 rap_rgap_map = pd.read_csv("RAP_RGAP_mapping.txt", sep='\t') rap_rgap_map.to_csv("id_mapping.csv", index=False)
常见问题解决方案：
- Parent关系缺失：检查原始文件是否包含完整的ID/Parent属性
- 坐标冲突：使用
gffutils.merge策略处理重叠特征 - 版本不一致：确保所有输入文件来自同一注释版本
这套流程不仅适用于水稻基因组，经过适当调整也可应用于其他作物的注释文件整合。在实际项目中，建议将完整流程封装为Snakemake或Nextflow工作流，实现可重复的自动化处理。
