我最近在用YGAP (Yeast Genome Annotation Pipeline) 做yeast基因组的注释。
YGAP生成的注释结果默认是bed格式不是gff格式,并且同时提供了embl和genbank的格式,因此如果想要得到gff格式的注释文件,就需要手动把格式转换一下。
YGAP (
http://wolfe.ucd.ie/annotation/
) 的使用方法很简单,提交序列上去,邮箱填一下,根据自己的需求设置一下选项就好了。傻瓜式操作,这里就不展开讲了。
接下来介绍如何将embl格式转换成gff3格式,使用的工具是agat
。
首先,agat是什么?
AGAT(Another GTF/GFF Analysis Toolkit)是一个GTF/GFF工具包,允许您执行几乎所有您可能想要完成的任务。AGAT是一套工具,能够处理任何GTF/GFF格式,并执行您可能需要的大多数任务。这使得AGAT在处理和分析基因组注释数据时具有极大的灵活性和实用性。
为什么这里说的是任何的GTF/GFF格式呢?难道不应该就俩格式吗?
因为:
GTF/GFF格式是用于描述和表示基因组特征的9列文本格式。自1997年以来,这些格式已经发展许多,并且尽管现在有明确定义的规范,但它们具有很大的灵活性,允许容纳各种信息。
然而,这种灵活性有一个缺点,即格式的种类繁多,包括GFF、GFF1、GFF2、GFF2.5、GFF3、GTF、GTF2、GTF2.1、GTF2.2、GTF2.5和GTF3。由于特定的期望,很多使用GTF/GFF格式的工具很难理解和区分所有的GFF/GTF格式。
说人话就是看起来都是gff
格式,但是每个gff
互相之间都有一些不同的地方。每次都要具体情况具体分析才能得到想要的结果。非常恼人。
(具体不同的点可以参考这个链接:https://agat.readthedocs.io/en/latest/gxf.html
)
之后我也有计划全部重新总结一下gxf
格式的变迁。(flag已立好)
无他,深受其扰尔。
agat
包含了一系列以agat
开头的工具,用于gxf
与各种格式之间的转换。
具体的工具的列表可以查看:https://agat.readthedocs.io/en/latest/index.html
agat安装
conda一键完成:
conda install agat
目前最新版为1.1.0
数据处理
注释得到的文件是embl.zip文件,先解压开
unzip embl.zip
进入到文件夹之后就可以看到单条染色体的embl格式的文件。
agat中用于embl和gff格式之间转换的工具是 agat_convert_embl2gff.pl
下面是它的帮助文档简化版:
Using standard /path/to/envs/daily/lib/perl5/site_perl/auto/share/dist/AGAT/config.yaml file
Name:
agat_converter_embl2gff.pl
Description:
The script takes an EMBL file as input, and will translate it in gff format.
Usage:
agat_converter_embl2gff.pl --embl infile.embl [ -o outfile ]
Options:
--embl Input EMBL file that will be read
--emblmygff3
Bolean - Means that the EMBL flat file comes from the EMBLmyGFF3
software. This is an EMBL format dedicated for submission and
contains particularity to deal with. This parameter is needed to
get a proper sequence id in the GFF3 from an embl made with
EMBLmyGFF3.
--primary_tag, --pt, -t
List of "primary tag". Useful to discard or keep specific
features. Multiple tags must be coma-separated.
-d Bolean - Means that primary tags provided by the option
"primary_tag" will be discarded.
-k Bolean - Means that only primary tags provided by the option "primary_tag" will be kept.
--throw_fasta Bolean - Means that you do not want to keep the fasta sequence at the end of the gff output.
-o, --output, --out, --outfile or --gff Output GFF file. If no output file is specified, the output will be written to STDOUT.
-h or --help Display this helpful text.
最基础的使用方法就是:
agat_convert_embl2gff.pl --embl my.embl -o my.gff3
如果不出意外的话,运行完这条命令你就得到了你想要的gff文件了。
如果你有很多个embl文件,那么一个循环就完事儿了:
ls *embl | while read id ; do agat_convert_embl2gff.pl --embl ${id} -o ${id%.*}.gff3 ;done
接下来只要把多个gff文件cat成一个总的gff3文件就可以啦。
我在使用的过程中遇到了个问题:
生成的每个gff3文件下面都会自带对应的fasta序列。
当然,这个理论上讲不是bug,因为gff3格式底部带序列是正常的做法,但是这不是我想要的,并且这个脚本提供的--throw_fasta
选项失效了。
bug修复
虽然 agat_convert_embl2gff.pl
提供了一个 --throw_fasta
的参数,但是实际上这个参数是无法使用的。如果你使用如下的命令:
agat_convert_embl2gff.pl --embl test.embl --throw_fasta -o test.gff
是会直接得到一个报错的:
Unknown option: throw_fasta
Failed to parse command line
我使用的agat的版本是 1.1.0
,是目前的最新版。
因此我已经把这个bug提交给官方了。
那么有什么办法可以解决这个问题吗?
毕竟当我有多个embl文件,转换成gff文件之后我是想要合并到一起的。我不想要中间穿插着fasta文件。
用sed删除掉底部fasta序列
我想到的第一个方法是用sed来完成这个去除掉fasta的操作:
ls *gff | while read id; do sed -i '/##FASTA/,$d' ${id} ; done
让GPT给你解释一下这条命令:
这条Linux命令是一个组合命令,用于批量处理目录中的GFF文件。它的功能是删除每个GFF文件中的"##FASTA"行及其后面的所有行。
下面是对这条命令的简要解释:
ls *gff:列出当前目录中所有以"gff"结尾的文件。
|:管道符号,用于将前一个命令的输出作为下一个命令的输入。
while read id:对于通过管道传递的每个文件名,将其赋值给变量id。
do:开始一个循环,对每个文件执行以下操作。
sed -i '/##FASTA/,{id}:使用sed命令编辑每个文件(由变量id表示)。
-i选项表示就地编辑文件。
'/##FASTA/,$d'是一个sed表达式,表示删除从包含"##FASTA"的行开始到文件末尾的所有行。
done:结束循环。
总的来说,这条命令搜索当前目录中的所有GFF文件,并从每个文件中删除"##FASTA"行及其后面的所有行。
Update:
修改配置文件以丢弃fasta
当我写完上述内容去吃个了饭回来后发现我提交的issue被回复了...
震惊,效率也太高了。
作者回复说, throw_fasta
这个参数已经不再生效了,如果想要丢掉自带的fasta,可以修改agat运行时使用的 config.yaml
文件。
运行:
agat config --expose
就会在当前文件夹生成一个config.yaml
文件,只要手动地把throw_fasta
这个参数给改成true就可以了:
# should fasta embedded in GFF/GTF be thrown away from input?
throw_fasta: true
再次运行的时候就会使用本地的config.yaml
文件,然后就可以把fasta序列给丢掉了。
萌哥碎碎念
- 今天是美国国庆,也是女朋友生日。
祝万万老师十八岁生日快乐!~(别问,问就是每年都是十八岁∠( ᐛ 」∠)_) - 今天放假在家。以后争取每次放假都发布一个更新。当然,flag立起来就是用来打的,我们将会看。
网友评论