R语者小case之——从GTF文件生成注释表格做基因ID转换

 基因的注释表格是经常需要用到的，可以从GTF文件中获得。用R可以简单地实现这个功能。

简易的GTF文件实际上可以认为是用制表符分隔为9列的TSV。
第一列是seqid, 通常是染色体编号；
第二列是source，表示信息的来源；
第三列是feature,表示类型；
第四和第五列分别是区间的起始位置与终止位置；
第六列是score, 软件提供的统计值；
第七列是strand, 代表正负链的信息, +表示正链，-表示负链，?表示不清楚正负链的信息；
第八列是phase，当第三列是CDS时，需要指定翻译时开始的位置，取值范围有0,1,2；
第九列是attributes, 表示属性，必须有gene_id和transcript_id这两个属性， 多个属性用分号分隔

通常基因的feature值为gene，而gene的注释信息集中在第九列。
OK，了解了以上信息，stop talking, show me the code.

library(stringr)
library(readr)
gtf_in <- read_delim("Homo_sapiens.gtf",
"\t", escape_double = FALSE, col_names = FALSE,
comment = "#", trim_ws = TRUE)
gtf_in <- gtf_in[gtf_in$X3=="gene",]#第三列feature是gene的保留下来
gene_info <- data.frame(str_split_fixed(unique(gtf_in$X9),";",5))
# 多个属性用分号分隔
head(gene_info)#瞅一眼结果
                         X1                X2                       X3                     X4
1 gene_id "ENSG00000223972"  gene_version "5"      gene_name "DDX11L1"   gene_source "havana"
2 gene_id "ENSG00000227232"  gene_version "5"       gene_name "WASH7P"   gene_source "havana"
3 gene_id "ENSG00000278267"  gene_version "1"    gene_name "MIR6859-1"  gene_source "mirbase"
4 gene_id "ENSG00000243485"  gene_version "5"  gene_name "MIR1302-2HG"   gene_source "havana"
5 gene_id "ENSG00000284332"  gene_version "1"    gene_name "MIR1302-2"  gene_source "mirbase"
6 gene_id "ENSG00000237613"  gene_version "2"      gene_name "FAM138A"   gene_source "havana"
                                                   X5
1  gene_biotype "transcribed_unprocessed_pseudogene";
2              gene_biotype "unprocessed_pseudogene";
3                               gene_biotype "miRNA";
4                             gene_biotype "lincRNA";
5                               gene_biotype "miRNA";
6                             gene_biotype "lincRNA";
gene_info2 <- unique(data.frame("gene_id"=str_split_fixed(gene_info$X1,'"',3)[,2],
                         "gene_name"=str_split_fixed(gene_info$X3,'"',3)[,2],
                         "gene_biotype"=str_split_fixed(gene_info$X5,'"',3)[,2]))
#再用str_split_fixed获取引号内的内容，这样就获得了gene_id、gene_name、gene_biotype构成的基因注释表格

通过这个注释表格，我们就可以做基因ID转换等很多工作了~
示例的GTF文件第九列比较规整，每个基因都有5个相同属性，如果gene的GTF文件第九列不规整，那又如何来做这个事情呢？这里获取的是gene的注释信息，那么对于transcript或者CDS又如何处理呢？TB搜索店铺：R语者，获取指导。