从Gencode解析GTF文件

Question

我使用该data.table包编写了一个脚本来解析GENCODE gtf文件的最后一列.对于那些不知道的人来说,该列包含一些键值项,每行由一个分号分隔.我正在使用的特定文件包含大约250万行.我索引前100个,然后前1000行只是为了测试脚本,输出正是我需要的.但是,尽管使用了该set功能,但运行时间并不像我预期的那么快.它是前100行的瞬间,但前1000行需要大约一两分钟.这是脚本.

#LOAD DATA.TABLE LIBRARY
require(data.table)
#READ GTF ANNOTATION FILE
info <- fread("gencodeAnnotation.gtf")
colnames(info)[9] <- "AdditionalInfo"
info <- info[1:1000]
#CREATE LIST OF 'KEYS' TO PARSE OUT
pars <- as.character(list("gene_id", "gene_type", "gene_status", "gene_name", " level ", "transcript_name", "transcript_id", "transcript_type", "transcript_support_level", "havana_gene"))
#NESTED FOR LOOP TO PARSE KEY-VALUE PAIR
for (i in 1:length(pars)) {
    for (j in 1:nrow(info)) {
        infoRow <- info[,tstrsplit(AdditionalInfo, ';', fixed = T)][j]
        headerCheck <- like(infoRow, pars[i])
        if (any(headerCheck) == TRUE) {
          keyVal <- length(tstrsplit(infoRow[[which(headerCheck == T)]], " ", fixed = T))
          set(info, i = j, j = toupper(pars[i]), value = tstrsplit(infoRow[[which(headerCheck == T)]], " ", fixed = T)[[keyVal]])   
        } else {
          set(info, i = j, j = toupper(pars[i]), value = NA)    
        }
    }
}

正如我之前所说,在前100行,1000行测试时,输出是完美的.根据代码,它必须遍历所有行乘以要添加的列数或其中的项目pars.我的问题是,我的脚本中缺少什么或者我可以进行哪些编辑以减少运行时间？以下是正在使用的gtf文件的链接:http://www.gencodegenes.org/releases/current.html.它是第一个标记为"综合基因注释"的链接.提前致谢.

每个行的样本如下:

gene_id ENSG00000223972.5; gene_type transcribed_unprocessed_pseudogene; gene_status KNOWN; gene_name DDX11L1; level 2; havana_gene OTTHUMG00000000961.2; remap_status full_contig; remap_num_mappings 1; remap_target_status overlap;

Answer 1

我发现readGFFbioconductor包rtracklayer的功能非常合适.

require(rtracklayer)
system.time(gtf <- readGFF("~/Downloads/gencode.v24.annotation.gtf", version=2L))
#    user  system elapsed 
#  34.558   1.541  36.737 
dim(gtf)
# [1] 2572840      26

您也可以只选择您喜欢的标签.

system.time(gtf_tags <- readGFF("~/Downloads/gencode.v24.annotation.gtf", version=2L, 
      tags = c("gene_id", "transcript_id")))
#    user  system elapsed 
#  16.830   0.733  17.780 
dim(gtf_tags)
# [1] 2572840      10