刘小泽写于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
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
网友评论