我试图从下面的sbml/xml文件解析信息
https://dl.dropboxusercontent.com/u/10712588/file.xml从这个代码
http://search.bioconductor.jp/codes/11172似乎可以通过
正常导入文件doc <- xmlTreeParse(filename,ignoreBlanks = TRUE)
但是我不能通过
恢复节点属性atrr <- xpathApply(doc, "//species[@id]", xmlGetAttr, "id")
或
xpathApply(doc, "//species", function(n) xmlValue(n[[2]]))
文件的一个节点在…
<species id="M_10fthf_m" initialConcentration="1" constant="false" hasOnly
SubstanceUnits="false" name="10-formyltetrahydrofolate(2-)" metaid="_metaM_10fth
f_m" boundaryCondition="false" sboTerm="SBO:0000247" compartment="m">
<notes>
<body xmlns="http://www.w3.org/1999/xhtml">
<p>FORMULA: C20H21N7O7</p>
<p>CHARGE: -2</p>
<p>INCHI: InChI=1S/C20H23N7O7/c21-20-25-16-15(18(32)26-20)23-11(7-22
-16)8-27(9-28)12-3-1-10(2-4-12)17(31)24-13(19(33)34)5-6-14(29)30/h1-4,9,11,13,23
H,5-8H2,(H,24,31)(H,29,30)(H,33,34)(H4,21,22,25,26,32)/p-2/t11-,13+/m1/s1</p>
<p>HEPATONET_1.0_ABBREVIATION: HC00212</p>
<p>EHMN_ABBREVIATION: C00234</p>
</body>
</notes>
<annotation>
...
我想检索物种节点内的所有信息,有人知道怎么做吗?
存在一个SBML解析库libSBML (http://sbml.org/Software/libSBML)。
这包括一个到R的绑定,允许在R中使用类似于
的代码直接访问SBML对象。document = readSBML(filename);
errors = SBMLErrorLog_getNumFailsWithSeverity(
SBMLDocument_getErrorLog(document),
enumToInteger("LIBSBML_SEV_ERROR", "_XMLErrorSeverity_t")
);
if (errors > 0) {
cat("Encountered the following SBML errors:n");
SBMLDocument_printErrors(document);
q(status=1);
}
model = SBMLDocument_getModel(document);
if (is.null(model)) {
cat("No model present.n");
q(status=1);
}
species = Model_getSpecies(model, index_of_species);
id = Species_getId(species);
conc = Species_getInitialConcentration(species)
每个可能的属性都有一个Species_get(NameOfAttribute)函数;与Species_isSet(NameOfAttribute)一起;Species_set(NameOfAttribute)和Species_unset(NameOfAttribute).
与任何SBML元素交互的API都是类似的。
libSBML版本包括R安装程序,可从
http://sourceforge.net/projects/sbml/files/libsbml/5.8.0/stable导航到您选择的操作系统和体系结构的R_interface子目录
libSBML的源代码发行版包含一个examples/r目录,其中有许多在r环境中使用libSBML与SBML交互的示例。
我想这取决于当您说要"检索"物种节点中的所有信息时的意思,因为检索到的数据可以强制转换为任意数量的不同格式。下面的示例假设您希望所有内容都放在一个数据框架中,其中每行是XML文件中的一个物种节点,列表示不同的信息片段。
当试图提取信息时,我通常发现使用列表比使用XML更容易。
doc <- xmlTreeParse(xml_file, ignoreBlanks = TRUE)
doc_list <- xmlToList(doc)
一旦它在列表中,你就可以找出物种数据存储的位置:
sapply(x, function(x)unique(names(x)))
[[1]]
NULL
[[2]]
NULL
[[3]]
NULL
[[4]]
[1] "species"
[[5]]
[1] "reaction"
[[6]]
[1] "metaid"
$.attrs
[1] "level" "version"
所以你只需要doc_list[[4]]
中的信息。看看doc_list[[4]]
的第一个组件:
str(doc_list[[4]][[1]])
List of 9
$ : chr "FORMULA: C20H21N7O7"
$ : chr "CHARGE: -2"
$ : chr "HEPATONET_1.0_ABBREVIATION: HC00212"
$ : chr "EHMN_ABBREVIATION: C00234"
$ : chr "http://identifiers.org/obo.chebi/CHEBI:57454"
$ : chr "http://identifiers.org/pubchem.compound/C00234"
$ : chr "http://identifiers.org/hmdb/HMDB00972"
$ : Named chr "#_metaM_10fthf_c"
..- attr(*, "names")= chr "about"
$ .attrs: Named chr [1:9] "M_10fthf_c" "1" "false" "false" ...
..- attr(*, "names")= chr [1:9] "id" "initialConcentration" "constant" "hasOnlySubstanceUnits" ...
这样就得到了包含在前8个列表中的信息,以及包含在属性中的信息。
获取属性信息很容易,因为它已经被命名了。下面将属性信息格式化为每个节点的数据帧:
doc_attrs <- lapply(doc_list[[4]], function(x) {
x <- unlist(x[names(x) == ".attrs"])
col_names <- gsub(".attrs.", "", names(x))
x <- data.frame(matrix(x, nrow = 1), stringsAsFactors = FALSE)
colnames(x) <- col_names
x
})
一些节点似乎没有属性信息,因此返回空数据帧。这引起了后来的问题,所以我创建了NAs的数据帧:
doc_attrs_cols <- unique(unlist(sapply(doc_attrs, colnames)))
doc_attrs[sapply(doc_attrs, length) == 0] <-
lapply(doc_attrs[sapply(doc_attrs, length) == 0], function(x) {
df <- data.frame(matrix(rep(NA, length(doc_attrs_cols)), nrow = 1))
colnames(df) <- doc_attrs_cols
df
})
当涉及到提取非属性数据时,变量的名称和值通常包含在同一个字符串中。我最初试图用一个正则表达式来提取这些名称,但它们的格式都不同,所以我放弃了,只是识别了这个特定数据集中所有的可能性:
flags <- c("FORMULA:", "CHARGE:", "HEPATONET_1.0_ABBREVIATION:",
"EHMN_ABBREVIATION:", "obo.chebi/CHEBI:", "pubchem.compound/", "hmdb/HMDB",
"INCHI: ", "kegg.compound/", "kegg.genes/", "uniprot/", "drugbank/")
此外,有时非属性信息仅作为值列表保存,如我上面所示的节点,而其他时候它包含在"notes"one_answers"annotation"子列表中,因此我必须包含if else
语句以使事情更加一致。
doc_info <- lapply(doc_list[[4]], function(x) {
if(any(names(x) != ".attrs" & names(x) != "")) {
names(x)[names(x) != ".attrs"] <- ""
x <- unlist(do.call("c", as.list(x[names(x) != ".attrs"])))
} else {
x <- unlist(x[names(x) != ".attrs"])
}
x <- gsub("http://identifiers.org/", "", x)
need_names <- names(x) == ""
names(x)[need_names] <- gsub(paste0("(", paste0(flags, collapse = "|"), ").+"), "\1", x[need_names], perl = TRUE)
#names(x) <- gsub("\s+", "", names(x))
x[need_names] <- gsub(paste0("(", paste0(flags, collapse = "|"), ")(.+)"), "\2", x[need_names], perl = TRUE)
col_names <- names(x)
x <- data.frame(matrix(x, nrow = 1), stringsAsFactors = FALSE)
colnames(x) <- col_names
x
})
为了将所有内容整合到一个数据帧中,我建议使用plyr
包的rbind.fill
。
require(plyr)
doc_info <- do.call("rbind.fill", doc_info)
doc_attrs <- do.call("rbind.fill", doc_attrs)
doc_all <- cbind(doc_info, doc_attrs)
dim(doc_all)
[1] 3972 22
colnames(doc_all)
[1] "FORMULA:" "CHARGE:" "HEPATONET_1.0_ABBREVIATION:" "EHMN_ABBREVIATION:"
[5] "obo.chebi/CHEBI:" "pubchem.compound/" "hmdb/HMDB" "about"
[9] "INCHI: " "kegg.compound/" "kegg.genes/" "uniprot/"
[13] "drugbank/" "id" "initialConcentration" "constant"
[17] "hasOnlySubstanceUnits" "name" "metaid" "boundaryCondition"
[21] "sboTerm" "compartment"
作为部分答案,文档使用了名称空间,而'species'是'id'名称空间的一部分。所以
> xpathSApply(doc, "//id:species", xmlGetAttr, "id", namespaces="id")
[1] "M_10fthf_c" "M_10fthf_m" "M_13dampp_c" "M_h2o_c" "M_o2_c"
[6] "M_bamppald_c" "M_h2o2_c" "M_nh4_c" "M_h_m" "M_nadph_m"
...
与id:species
和namespaces="id"
不同,