报错
#1
ERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
The specified gene identifier attribute is 'gene_id'
An example of attributes included in your GTF annotation is 'transcript_id "rna-NC_001941.1:1..68";'
The program has to terminate.
ERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file
#2
使用featucounts,程序崩溃报错GTF文件。
使用的文件是NCBI RefSeq中未经任何修改的参考基因组和GTF注释
使用独立的featurecots可执行文件(2.0.1)或在Rsubread包(2.2.2)中的featurecots函数,得到了相同的崩溃和错误消息:
ERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
The specified gene identifier attribute is 'gene_id'
An example of attributes included in your GTF annotation is 'gene_id ""; transcript_id "unknown_transcript_1"; gbkey "tRNA"; product "tRNA-Phe"; exon_number "1"; '
The program has to terminate.
原因
有些GTF文件如此:NC_000001.11 BestRefSeq exon 11874 12227 . + . transcript_id "rna0"; gene_id "gene0"; gene_name "DDX11L1";
而featureCounts 需要输入的GTF格式为:NC_012920.1 hg38.GenePred transcript 15956 16023 . - . gene_id "gene60958"; transcript_id "rna171196"; gene_name "gene60958";
差异在于gene_id的顺序问题。
此外,应注意RNAseq使用的参考基因组.fa与参考注释文件.gtf文件染色体ID名称应该一致。
解决 1
#似乎使用ucsc开发的工具即可;
gff3ToGenePred GCF_000001405.38_GRCh38.p12_genomic.gff hg38.GenePred
genePredToGtf database hg38.GenePred hg38.gtf
# 去其他数据库找找,或者在已有基础上自己写脚本更换顺序
对于GTF的实际使用而言,只需要
exon
区间信息就可以区分不同的转录本了,而且在定量的过程中,也是只需要参考exon的位置信息。所以,自己编写转换脚本,只需要保留eoxn
信息。
另外一个问题就是,第九列提供哪些属性。根据我的经验,只需要以下6个属性:1. gene_id 2. gene_name 3. transcript_id 4. transcript_name 5. gene_type 6. transcript_type
gene_id
可以用来保存不同数据库中的基因ID,比如NCBI Entrez Id, Ensembl gene Id, 当然也可以和gene_name属性保持一致;gene_name
属性用来保存gene symbol, 相比id, symbol在文章中的使用频率更高。
transcript_id
和transcript_name
表征转录本的id和名称,可以是RefSeq ID,也可以是Ensembl transcript id, 用于区分不同的转录本。
gene_type
和transcript_type
表示基因和转录本的类型,比如是protein_coding
, 还是lncRNA
,rRNA
等。在分析时,我们通常会根据类型选择其中的部分转录本来分析,比如只分析蛋白编码的转录本。
以上6种属性就能够满足几乎100%的场景,对于不同数据库中的文件,只需要自己写脚本提取这些信息,就可以了。
基因组注释文件(GFF,GTF)下载数据库
详解GFF转换为GTF文件
解决 2
如屏幕输出所示,NCBI RefSeq的GTF注释包含空的gene_id值。这在featucounts中是不允许的。
您可以考虑删除带有空gene_id值的行。作者表示会在下一个版本中找到此类问题的通用解决方案。
zcat raw.gtf.gz | grep gene_id | wc -w #检查原始gtf注释文件gene_id个数
zcat raw.gtf.gz | grep gene_id\ \"\"\; | wc -w #检查空gene_id值的行数
sed -i -e '/gene_id\ \"\"\;/d' raw.gtf #删除包括空gene_id值的行 最好重命名一下 如:raw.no.na.gtf
cat raw.no.na.gtf | grep gene_id\ \"\"\; | wc -w #修改后的检查一下 检查包括空gene_id值的行 0为已删除掉
featureCounts -T 5 -t exon -g gene_id -a raw.no.na.gtf -o raw.counts.txt raw.sorted.bam