用agat把embl格式转换成gff格式

我最近在用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提交给官方了。

https://github.com/NBISweden/AGAT/issues/371

那么有什么办法可以解决这个问题吗?

毕竟当我有多个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/,d'{id}:使用sed命令编辑每个文件(由变量id表示)。

-i选项表示就地编辑文件。

'/##FASTA/,$d'是一个sed表达式,表示删除从包含"##FASTA"的行开始到文件末尾的所有行。

done:结束循环。

总的来说,这条命令搜索当前目录中的所有GFF文件,并从每个文件中删除"##FASTA"行及其后面的所有行。

Update:

修改配置文件以丢弃fasta

当我写完上述内容去吃个了饭回来后发现我提交的issue被回复了...

震惊,效率也太高了。

https://github.com/NBISweden/AGAT/issues/371

作者回复说, 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序列给丢掉了。

萌哥碎碎念

  1. 今天是美国国庆,也是女朋友生日。
    祝万万老师十八岁生日快乐!~(别问,问就是每年都是十八岁∠( ᐛ 」∠)_)
  2. 今天放假在家。以后争取每次放假都发布一个更新。当然,flag立起来就是用来打的,我们将会看。
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 193,812评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,626评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,144评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,052评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 60,925评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,035评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,461评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,150评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,413评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,501评论 2 307
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,277评论 1 325
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,159评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,528评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,868评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,143评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,407评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,615评论 2 335

推荐阅读更多精彩内容