vlambda博客
学习文章列表

R语言提取GTF文件信息

GFF和GTF是两种最常用的数据库注释格式,在信息分析中建库时除了需要fasta文件一般还会需要这两种文件,提取需要的信息进行注释。

GFF全称为general feature format,这种格式主要是用来注释基因组。

GTF全称为gene transfer format,主要是用来对基因进行注释。

GTF文件由9列数据组成,这两种文件的前8列都是相同的(一些小的差别),gtf文件是以tab键分割的9列组成,以下为每一列的对应信息:

1)  seq_id:序列的编号,一般为chr或者scanfold编号;

2)  source: 注释的来源,一般为数据库或者注释的机构,如果未知,则用点"."代替;

3)  type: 注释信息的类型,比如Gene、cDNA、mRNA、CDS等

4)  start: 该基因或转录本在参考序列上的起始位置;

5)  end: 该基因或转录本在参考序列上的终止位置;

6)  score: 得分,数字,是注释信息可能性的说明,可以是序列相似性比对时的E-values值或者基因预测是的P-values值,"."表示为空;

7)  strand: 该基因或转录本位于参考序列的正链(+)或负链(-)上;

8)  phase: 仅对注释类型为"CDS"有效,表示起始编码的位置,有效值为0、1、2 (对于编码蛋白质的CDS来说,本列指定下一个密码子开始的位置。每3个核苷酸翻译一个氨基酸,从0开始,CDS的起始位置,除以3,余数就是这个值,,表示到达下一个密码子需要跳过的碱基个数。该编码区第一个密码子的位置,取值0, 1, 2。0表示该编码框的第一个密码子第一个碱基位于其5'末端;1表示该编码框的第一个密码子的第一个碱基位于该编码区外;2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外;如果Feature为CDS时,必须指明具体值。);

9)  attributes:一个包含众多属性的列表,格式为"标签=值"(tag=value),标签与值之间以空格分开,且每个特征之后都要有分号;(包括最后一个特征),其内容必须包括gene_id和transcript_id。以多个键值对组成的注释信息描述,键与值之间用"=",不同的键值用“


gtf文件介绍主要来自:

GTF基因注释文件详解_人工智能_宁生信-CSDN博客:

https://blog.csdn.net/sinat_38163598/article/details/72851239



如果说要从gtf文件中提取有用的东西呢?比如gene_id,gene_name (symbol),gene_type......,所以我来介绍一个超好用的GTF加载器:rtracklayer。


The rtracklayer package is an interface (or layer ) between R and genome browsers. Its main purpose is the visualization of genomic annotation tracks, whether generated through experimental data analysis performed in R or loaded from an external data source. The features of rtracklayer may be divided into two categories: 1) the import/export of track data and 2) the control and querying of external genome browser sessions and views.

该包也属于Bioconductor (以R语言为平台的一个高通量基因组数据的分析工具),关于Bioconductor,徐州更有一篇简书介绍的超好:https://www.jianshu.com/p/f528f1fad97e,里面有很多使用的包,可以在之后使用。所以他的安装就需要Bioconductor类包的特定安装方式:

BiocManager::install("rtracklayer")

来看 rtracklayer,它其实不仅可用于gtf文件的导入导出,还可以操作其他多种数据格式(BED,WIG,bigWig......)。至于gtf文件,请去ensemble下载。

gtf <- rtracklayer::import('Homo_sapiens.GRCh38.99.gtf.gz.gtf')  #导入gtf_df=as.data.frame(gtf) #转化为矩阵colnames(gtf)  #查看一下包含的内容 [1] "seqnames" "start"  [3] "end" "width"  [5] "strand" "source"  [7] "type" "score"  [9] "phase" "gene_id" [11] "gene_version" "gene_name" [13] "gene_source" "gene_biotype" [15] "transcript_id" "transcript_version" [17] "transcript_name" "transcript_source" [19] "transcript_biotype" "tag" [21] "transcript_support_level" "exon_number" [23] "exon_id" "exon_version" [25] "protein_id" "protein_version" [27"ccds_id"                 

看到这我觉得可以开心了,至少我是很开心了。。。

直接用"$"提取你想要的列就可以了。其中gene_id是ensemble_id,gene_name是symbol。但是注意,提取出的信息包含重复,所以可以用duplicated函数进行过滤,比如:

test[!duplicated(test),] #删掉所有列上都重复的


好了,祝你结果顺利咯~