Re: Reading xml pdb molecular data: Expert advice for xml/R needed
Bert Gunter <[email protected]> Sun, 5 Apr 2026 14:09:25 -0700
| Newsgroups | gmane.comp.lang.r.general |
|---|---|
| Message-ID | <CAGxFJbT7-oGSYNyNVjuHbqMnH+1ACfY9dfoUPaYgcq=CGcs0gA@mail.gmail.com> |
I think you would do better posting this on the appropriate Bioconductor forum, no? Cheers, Bert On Sun, Apr 5, 2026 at 12:44 PM Leo Mada via R-help <[email protected]> wrote: > Dear R-Users, > > I want to extract the molecular structure from PDB XML files using R. I > came up with 3 variants of R code, but none is ideal. I will describe the 3 > ways below and hope to find some expert advice how to do it best. > > Functional versions of the 3 variants of code are posted on GitHub: > https://github.com/discoleo/Rpdb/issues/4 > > Various example files can be downloaded from: > https://www.rcsb.org/ > - after searching for the PDB code (e.g. 4wrg), see button: > "Download File" => "PDB XML.gz format"; > > Examples of various sizes are: > > # Atoms: 2485 (Smallish) > doc = xmlParse("4wrg.xml") > > # Atoms: 12500 (Medium) > doc = xmlParse("3fxi.xml") > > # Atoms: 25450 (Large) > # - contains 10 models of 2545 atoms each; > doc = xmlParse("7k3g.xml") > > # Atoms: 1,359,288 (Very Large) > # - contains 378 "chains" of 3596 atoms each; > doc = xmlParse("6x63.xml") > # Note: slow to open & very slow to convert to data.frame (~30 minutes); > > My 3 approaches follow now: > > 1. Using package XML > Unfortunately, I do not have much experience with this package. I used > mostly the xml2 package. However, this version works best - though still > NOT perfect. > > The effort to convert to a data.frame is small and it is pretty fast, too. > However, I did not find a way to extract a data.frame which includes one of > the xml attributes. The atom id is stored as an attribute in the > "atom_site" parent. I therefore parse first the "atom_site" subtree to add > the value of the id attribute as a proper child-node to each of the > "atom_site" nodes: > > ### Atoms: > at = getNodeSet(doc, > "/PDBx:datablock/PDBx:atom_siteCategory/PDBx:atom_site[@id]") > > tmp = lapply(seq_along(at), function(id) { > idN = xmlAttrs(at[[id]])[["id"]]; > addChildren(at[[id]], kids = idN); > }) > > # Extract now the data.frame: > tmp = xmlToDataFrame(nodes = tmp, > homogeneous = TRUE, collectNames = F) > > It works "fastesh", except for the last example, where it takes ~30 > minutes (but ultimately did complete). > > > 2. Using package xml2 > Involves explicitly reading all sub-nodes of the "atom_site" node. Much > more code to write (but is already done on GitHub). > > However, it has its own quirks: once I iterate over each "atom_site" node, > I need to convert it to text and read the small subtree as xml again. > Otherwise, it would take forever to process the whole tree. > > I tried as_list as well (on the "whole" xml-branch), but it is slower and > converting it to a data.frame is even more cumbersome (and slow). > > 3. Using package xslt > Seemed a reasonable approach - but I still do NOT know how to convert it > to a data.frame. I generate a quasi-csv, which needs some post-processing > to become a correct csv and then read the csv. > > The xslt stylesheet is included on the GitHub issue page - although I am a > little bit rusty with xslt. It may be possible to improve it. > > # - TODO: convert efficiently to data.frame; > style = read_xml(xs) # see GitHub; > tmp = xml_xslt(doc, style) > tmp = as.character(tmp) > tmp = substr(tmp, 39, nchar(tmp)) # Remove "<xml ...>" > tmp = gsub("[\n \r\t]+", "", tmp) > tmp = gsub(">", "\n", tmp) > tmp = read.csv(text = tmp, header = FALSE) > head(tmp) > > > FINAL REMARKS > > I prefer to extract the actual "id" and not infer anything about it. These > ids were not continuous in the legacy pdb format. Molecular chains had > chain terminators with id, but those were not atoms. The next atom would > therefore skip 1 id. > > The legacy pdb format had fixed width fields and it ran out of pdb codes > some 2 years ago. The xml pdb format seems a reasonably accessible way to > read newer pdb files. But I struggle to find the ideal approach to read > these files. Some expert advice is well appreciated. > > It may be warranted to improve some of the existing xml packages as well. > The Protein Data Bank contains over 250,000 protein structures determined > experimentally and even more models. It is a large resource to explore and > I find the R language more versatile than traditional tools used by > biochemists. > > There is a newer text-based pdb format. I started to read the > documentation some 2 years ago - but it was like reading a horror story - > and this is not only a joke! E.g. it is described how to loop over the > atoms-rows and duplicate the header for each row: in the statistics > community this is commonly known as csv-style data with headers. If anyone > is interested, I will search for the documentation and post a link to it. > > Sincerely, > > Leonard > > > [[alternative HTML version deleted]] > > ______________________________________________ > [email protected] mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > https://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > [[alternative HTML version deleted]]