在R语言中读取GTF文件的最好方法

下载的TCGA数据是没有注释的,需要从ensemble上面下载GTF文件,现在需要把GTF文件读入R

第一种方法rtracklayer::import

1
2
3
4
5
source("https://bioconductor.org/biocLite.R")
biocLite("rtracklayer")
biocLite("SummarizedExperiment")
gtf1 <- rtracklayer::import('Homo_sapiens.GRCh38.90.chr.gtf')
gtf_df <- as.data.frame(gtf1)

最终读入27个变量,2612129个观测,测试一下显示的不错

1
2
test <- gtf_df[1:5,]
View(test)

取出我需要的gene_id,gene_biotype,gene_name,成为新的数据框

1
geneid_df <- dplyr::select(gtf_df,c(gene_name,gene_id,gene_biotype))

第二种方法read.table

1
gtf2 <- read.table('Homo_sapiens.GRCh38.90.chr.gtf', header = FALSE, sep = '\t')

最后发现,速度奇慢无比,取消参考了Y叔的帖子增加参数,读入成功,需要1分钟,读入9个变量
根据GTF画基因的多个转录本结构

1
2
gtf2 <- read.table('Homo_sapiens.GRCh38.90.chr.gtf', stringsAsFactors = F, header = FALSE, sep = '\t',comment = "#")
test2 <- gtf2[1:5,]

第三种方法readr包里面的read_table2

速度很快,有进度条,读入了18个变量,但是最gene_biotype显示不对

1
2
3
4
gtf3 <- readr::read_table2("Homo_sapiens.GRCh38.90.chr.gtf",comment = "#")
gtf_df3 <- as.data.frame(gtf3) #转成data.frame
test3 <- gtf_df3[1:5,]
View(test3)

第四种方法read.table,fread

读入数据跟read.table一样,

1
2
3
4
library(data.table)
genes <- fread("Homo_sapiens.GRCh38.90.chr.gtf")
test4<- genes[1:5,]
View(test4)

上一步用时12秒,速度极快,只是读入后没有名字,需要自己添加
读入9个变量,跟read.table一样

1
setnames(genes, names(genes), c("chr","source","type","start","end","score","strand","phase","attributes") )

可选步骤,type那一项有转录本和gene,选取gene,其他有很多转录本

1
genes <- genes[type == "gene"]

我要的信息粗存在attributes中,构建函数提取

1
2
3
4
5
6
7
8
extract_attributes <- function(gtf_attributes, att_of_interest){
att <- strsplit(gtf_attributes,"; ")
att <- gsub("\"","",unlist(att))
if(!is.null(unlist(strsplit(att[grep(att_of_interest, att)], " ")))){
return( unlist(strsplit(att[grep(att_of_interest, att)], " "))[2])
}else{
return(NA)}
}

举例子解释,加入我们需要gene_id

1
2
gtf_attributes <- genes[1,9]
att <- strsplit(gtf_attributes,"; ")

这一步提示non-character argument,后面我查询才知道, 对于第九列的第一行的的获取两种方式返回的格式不一样

1
2
genes[1,9] #data.frame
genes$attributes[1] #character

而strsplit的对象是character,重来

1
2
3
gtf_attributes <- genes$attributes[1]
att <- strsplit(gtf_attributes,"; ")
att <- gsub("\"","",unlist(att))

grep获取位置,获取到内容后分隔,得到的第二个元素是我们需要的

1
unlist(strsplit(att[grep("gene_id", att)], " "))[2]

利用函数获得基因列表

1
genes$gene_id <- unlist(lapply(genes$attributes, extract_attributes, "gene_id"))

如果我们选择使用第一种方法rtracklayer::import,这些事情都是不需要做的!!!

#总结:rtracklayer::import最简单,fread最快,read.table可选,read_table没戏!

------ 本文结束------