【发布时间】:2016-05-05 22:55:53
【问题描述】:
我使用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;
【问题讨论】:
-
您的最小可重现示例在哪里?
-
i = j和j = i-- 你真的想搞混这样的事情吗?此外,"NA"不是您通常想在 R 中使用的东西;也许是NA。 -
我无法真正重现数据集,我提供了它的链接。您必须从网页上查看它的外观。它不像我可以复制的一堆值,它是一个像“GENE_ID”这样的键,后跟一个字符串或某些情况下的数字,用空格分隔。每个“键值”(可能每行 8 个或 9 个)都由分号分隔。如果你只是下载文件并读入,你可以只索引出一小部分并测试代码。 @eddi。
-
@abbas786 你给出了 0 个为什么你不能创建一个简单、小、可重复的例子的理由。
-
我编辑了问题并复制了该列的第一行。这有帮助吗? @eddi
标签: r performance data.table bioinformatics gtfs