agat：基于perl的GTF/GFF基因注释文件处理工具
写在前面的
以前在做微生物功能基因组学研究时，因为微生物基因组小，基因结构相对简单。所以许多关于微生物功能基因组分析过程中我对基因注释文件（GTF/GFF）的处理都是自己很快写写shell或perl的命令完成分析。后来开始做人类基因组分析后，逐渐发现以前针对微生物基因注释文件处理的命令在面对基因元件及结构多样化的人类基因组时心有余而力不足。如果重新针对人类基因组的GTF和GFF文件重构相关代码，即使有AI机器人帮助，也会显得效率低下，因为不熟悉文件信息是最大限制因素。后来借助Manus AI我找到了一个高效处理高级物种基因注释文件的神器。
agat介绍
agat软件全称是Another GFF Analysis Toolkit (AGAT)，该软件由瑞典国家生物信息中心（National Bioinformatics Infrastructure Sweden NBIS）团队开发。软件主体基于perl语言编写，已经过开发者团队多年对超过42种不同类型的GTF/GFF文件的测试，能够做到排查、订正并补充各种相关文件中缺失的特征及属性等信息。从而生成完整的或经过排序的标准的格式。agat与其说是一个软件，不如说是整合了许多处理注释文件的工具包。例如它足完成与GTF/GFF和对应的fasta文件相关的数据清洗、转换格式、合并文件、调整内容、筛选信息、提取FASTA序列、添加新数据等。并且开发者团队在长期保持对软件的更新、补充和维护。目前最新版本是1.5.1。
安装软件
agat的安装建议通过conda完成，特别建议单独为该软件创建新的独立虚拟环境，因为我在已有的以perl为依赖的环境中安装agat时碰见了Solving environment: - Killed的报错。经过调查发现应该是agat依赖的perl程序包和我已经安装的其他perl程序包有版本冲突。当然有已有的环境中安装的perl依赖很多，我也没有查出到底是哪个冲突了。
#创建环境conda create -p ~/software/miniconda3/envs/agat perl=5.32.1conda activate ~/software/miniconda3/envs/agatconda install -c bioconda agat#查看agat软件中的功能或工具ls ~/software/miniconda3/envs/agat/bin/*pl
准备分析数据
在UCSC官网http://genome.ucsc.edu/cgi-bin/hgTables下载人类基因组T2T版本的GTF注释文件，下载后解压缩。人类基因组端到端注释文件大约753Mb.
wget https://hgdownload.soe.ucsc.edu/goldenPath/hs1/bigZips/genes/hs1.ncbiRefSeq.gtf.gzzcat hs1.ncbiRefSeq.gtf.gz > hgT2T.ncbiRefSeq.gtfdu -ha hgT2T.ncbiRefSeq.gtf753Mb hgT2T.ncbiRefSeq.gtf#挑选chr1做测试awk -v OFS='\t' '{if($1=="chr1"){print $0}}' hgT2T.ncbiRefSeq.gtf > hgT2T.ncbiRefSeq.chr1.gtf
统计gtf/gff文件
统计gtf文件一般是统计文件内转录本、基因、外显子等元件的数量、长度分布信息。可以使用如下命令完成。在该命令中建议带上-p参数可以自动生成基因元件长度分布信息的pdf图片。
~/software/miniconda3/envs/agat/bin/agat_sp_statistics.pl -i hgT2T.ncbiRefSeq.gtf -o hgT2T.ncbiRefSeq.stat -p
输出的.stat文件内容超过了200行，详细统计了转录本或非转录本的基因元件数量等信息。在输出的pdf文件夹内包含了多个基因元件的长度分布信息。
为gtf文件增加带有内含子信息
常见的gtf文件中不包含对内含子的注释信息，因为内含子默认是处于两个相邻外显子之间，传统的提取内含子信息的手段依赖自己手写命令或编程，在agat的帮助下这个过程可以简单借用工具完成。下面的代码中我借助shell命令提出第三列排序去重后观察gene feature内是否增加了intron这个特征，以判断是否成功添加内含子信息。
cut -f 3 hgT2T.ncbiRefSeq.chr1.gtf | sort | uniq | grep -v \#3UTR5UTRCDSexonstart_codonstop_codontranscript~/software/miniconda3/envs/agat/bin/agat_sp_add_introns.pl --ref hgT2T.ncbiRefSeq.gtf --out hgT2T.ncbiRefSeq.introns.gtfcut -f 3 hgT2T.ncbiRefSeq.introns.chr1.gtf | sort | uniq | grep -v \#3UTR5UTRCDSexongeneintronstart_codonstop_codontranscript
生成内含子bed文件
如果要只生成某一特定feature的bed文件，事先要借助shell命令从总gff/gtf文件中提取出只含有该特征的子文件作为subgtf。然后再利用agat转换为bed。我用1号染色体做测试时发现我提取的原始intron有18万行，但是最后转成bed的intron只有1万6千行。在排除问题的过程中我发现这是因为1个基因有多个转录本，不同的转录本可能共享了多个内含子，对于这些内含子而言它们的信息在gtf文件中是冗余出现的，agat在转换bed的过程中对于冗余出现的内含子都去除了重复，并且对于去重后的内含子如果发现有任意内含子之间存在overlap，那么就合并带有overlap的内含子为一行bed信息。同时agat可能也会过滤掉一些不完整的内含子信息（这或许解释了为什么除去数量冗余和合并位置存在overlap的内含子后还比转完bed后的内含子多了3000个的原因，当然读者可以更详细补充这里可能的其他原因，总之agat保留的1万6千个内含子是非常确定的内含子信息，当然也欢迎读者解答是否还存在第二种途径是将全部的gtf文件转bed后在bed中保留各个feature，最后只挑选内含子，但是这个途径如何处理那种在A转录本里面是外显子但是在B转录本里面的内含子的区域呢）。
#只提取intron的gtf部分awk -v OFS='\t' '{if($3=="intron"){print $0}}' hgT2T.ncbiRefSeq.introns.chr1.gtf > hgT2T.ncbiRefSeq.introns.only.chr1.gtf#转换intron的信息为~/software/miniconda3/envs/agat/bin/agat_convert_sp_gff2bed.pl --gff hgT2T.ncbiRefSeq.introns.only.chr1.gtf --sub intron -o chr1.sub.introns.only.bed#计数bed文件函数wc -l chr1.sub.introns.only.bed16,235cut -f 1,4,5 hgT2T.ncbiRefSeq.introns.only.chr1.gtf | bedtools sort -i | uniq | bedtools merge -i - | wc -l19,361
提取蛋白序列
相比于只能提取核苷酸序列的bedtools而言，agat可以直接提取蛋白序列。用下面的命令可以提取cds或者gene的蛋白序列，如果有特殊需求，如果不提取蛋白序列，可以取消--protein选项，并使用--type参数指定要提取的特征成为核苷酸序列。
~/software/miniconda3/envs/agat/bin/agat_sp_extract_sequences.pl -g hgT2T.ncbiRefSeq.introns.chr1.gtf -f hgT2T.fa --protein --output cds.faaagat其他功能
因为该软件包含的功能较多，详细可参考官网https://agat.readthedocs.io/en/latest/
注意
在我测试过程中，我发现agat对于计算机内存有一定需求。建议用户在处理人类基因组GTF/GFF文件时确保有足够资源的内存，此外在提取序列的过程中，似乎对长度还有一定限制。