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("&gt;", "\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]]