今天有一个需求,就是要将gtf中的转录本转成fasta序列,一开始是想着用bedtools getfasta实现,awk取出来坐标做成bed文件输入bedtools,但是结果发现bedtools是单纯按照坐标取出来的,也懒得自己写脚本取了,搜一下发现cufflinks中有个程序可以实现。
如上图所示,“ENSMUST00000082908.1”转录本是这两个exons,取出这个转录本的fasta序列其实就是这两个exons对应的序列位置,需要把两个序列连起来。比如说,exon1是 chr 1 10-50 : AAAAAAAAAA; exon2是chr 1 70-80: TTTTTTTTT(仅做例子)。那么取出来的序列为:AAAAAAAAAATTTTTTTTT。(exons中间空出的部分并没有真的转录出来)。
gffread可以直接实现这个功能,这来自于cufflinks(一直不知道这个老软件竟然还有这个功能),直接conda install cufflinks之后即可使用gffread。使用如下代码即可转换:
gffread transcripts.gtf -g reference.fasta -w transcripts.fasta
转出来效果:
使用:
gffread -h
即可查看所有参数。
网友评论