115-GTF与GFF的小趣事
刘小泽写于19.5.14 最近在分析单细胞数据时,发现了有的数据需要将两个genome、两个GTF组合起来分析,但其中一个又是不常见的物种,只能自己去下载gff,然后转换,当然转换过程也有一些小坑
GFF与GTF的概念
GFF:General Feature Format,目前常用3版本,即GFF3。主要用来注释基因组,包含tab分隔的9列:
1. seq_id:序列编号(可以是chr或者scaffold) 2. source:注释来源(一般是数据库或机构),未知用.表示 3. type: 注释类型,如Gene、cDNA、mRNA、CDS等 4. start:基因/转录本在参考序列起始位置(1-based) 5. end:(同上)终止位置 6. score:(数字表示)得分,如序列相似性的E值或者基因预测的P值,未知用.表示 7. strand:基因/转录本在参考序列的+/-链 8. phase:当第3列type为编码蛋白的CDS时,这一列才有意义。表示下一个密码子开始的位置,或者说到下一个密码子需要跳过的碱基数,用0、1、2表示 # 其中,前8列在gff的三个版本(gff1、gff2、gff3)中记录的信息都是一样的,只是名称有差异:比如第一列虽然都记录序列编号,但gff1叫seqname,gff2叫reference sequence;type在gff1、gff2中叫feature;phase在gff1、gff2中叫frame 9. attributes:键值对表示的属性,格式是:key=value,不同的键值对用分号分隔;一个键内的多个值用逗号隔开。其中有一些已经定义好的键名:如 ID表示注释类型(也就是type的名称); Name:注释名称,可以重复; Parent:type属于的上一级名称,如exon上一级是transcript,而transcript的上一级是gene
GTF:General Transfer Format,主要对基因进行注释,目前主要使用gtf2。与GFF有两个主要的地方存在不同之处:
- 第三列的feature中一定包含CDS、start_condon、stop_codon
- 第九列的attribute中必须以gene_id和transcript_id开头,并且键值对之间用空格隔开,而不是用等号;另外每个键值对之间用分号分开,并且最后一个键值对的末尾也要有分号
NC_020166.2 Gnomon exon 2398 3338 . - . gene_id "1"; transcript_id "1.1"; NC_020166.2 Gnomon exon 3781 3884 . - . gene_id "1"; transcript_id "1.1";
GTF与GFF的转换
做项目时发现,自己研究的物种没有gtf,只有gff。不要管,只管下载下来,然后进行转换即可,而转换工具也有不少
工具一:gffread
# 它是由cufflinks开发的,使用如下
gffread -T xxx.gff -o xxx.gtf
# 生成的结果
NC_000001.11 BestRefSeq exon 11874 12227 . + . transcript_id "rna0"; gene_id "gene0"; gene_name "DDX11L1";
生成的GTF文件中只有exon
与CDS
信息,第九列只包含transcript_id
、gene_id
、gene_name
工具二:genometools中的gff3_to_gtf
# 工具下载安装
http://genometools.org/pub/binary_distributions/gt-1.5.10-Linux_x86_64-64bit-complete.tar.gz
# 它会将一些非exon、CDS的gene_type作为warning信息,可以直接不输出这些无用信息
gt gff3_to_gtf xxx.gff >xxx.gtf 2>/dev/null
不知道其他应用场景,反正在分析10X的时候,利用cellranger mkgtf过程中,使用genometools生成的gtf是不报错的,而gffread会报错
如果连gff也没有怎么办?
做cellranger发现,有时会用到组合两个物种的注释信息,比如其中一个是人,另一个是病毒,人的信息好找,但是病毒的只能找到基因组序列,没有注释信息,对于一些小型的序列可以选择手动去自己做。
另外网上还有一些脚本可以提供参考:例如https://github.com/jorvis/biocode/blob/master/gff/convert_genbank_to_gff3.py
但是自己构建GTF时,一般是先在文本编辑器编辑好,然后再cat >xxx.gtf
到命令行中。这个过程一定要注意使用tab分隔,即使数据表明看上去像也不可信
检查的方法有两个:
一个是直接运行cellranger mkgtf
看是否报错,可能会提示:GTF的行数不对;
另一个是直接检查:awk -F '\t' '{print NF}' mcv.gtf
,如果显示1,那么说明没有tab分隔,而是用的多个空格,利用sed可以把这些空格替换成tab:
sed 's/ \+ /\t/g' inputfile > outputfile