Add option for `gemmi` backend for PDB reading

[marinegor]

Is your feature request related to a problem?

It's more of a suggestion.

Describe the solution you'd like

Pull #4712 is about to introduce gemmi-based backend for mmCIF reading. It wasn't intended but, MMCIFReader/Parser can also read pdb files:

import MDAnalysis as mda
mda.Universe('testsuite/MDAnalysisTests/data/4E43.pdb', format='mmcif')
# <Universe with 1877 atoms>

since gemmi, being de-facto standard for RCSB stuff, can read both into the same gemmi.Model interface. Perhaps it's worth including an option to read pdb with gemmi backend then? I'm not sure how to technically do that though -- for instance, we could add format = 'pdb_gemmi' to MMCIFParser, and document this behavior in PDBParser.

Describe alternatives you've considered

Another option would be to just add a paragraph in documentation about this in MMCIFParser and/or PDBParser without changing any code -- I'd be probably fine with that, but just don't want this information to stay undocumented.

@MDAnalysis/coredevs do you have an opinion on that?

[orbeckst]

Can you do a performance comparison between reading with our PDBParser/Reader vs gemmi?

If gemmi is comparable or faster then we just have to do a feature comparison to see if we should make gemmi the default (and possibly retire our version).

[orbeckst]

For format: perhaps expand its syntax along the lines of including an optional library for parsing like gemmi|pdb or chemfiles|xtc .

 format="<library>|<fileformat>"

We could introduce "*" for the file format and then chemfiles|* means "use chemfiles for anything it can read". Similarly, gemmi|* then says, "use gemmi for CIF and PDB files".

---

I originally wanted to use : as separator but this will collide with URIs such as the ones we use for streaming such as "localhost://8888", which could be written explicitly with the imdclient as the library as "imdclient|localhost://8888".

Or does "[imdclient]localhost://8888" or "[gemmi]*" look better??? Has the disadvantage that [ and ] may occur directly in URIs as reserved chars (not %-encoded) (RFC 3986, Section 2.

[orbeckst]

How does gemmi deal with "PDB trajectories", i.e., multi-MODEL files, often with hundreds of frames? This is still the lowest common denominator for a MD trajectory, I thik...

[marinegor]

@orbeckst

How does gemmi deal with "PDB trajectories", i.e., multi-MODEL files, often with hundreds of frames?

I think it does, since it handles NMR depositions as well, which are multi-model. I haven't implemented that in my PR, but it can be added in the future, and this moment can also be a moment to retire a legacy PDB reader.

Or does "[imdclient]localhost://8888" or "[gemmi]*" look better???

I would not be in favor of neither [] nor | based syntax, and rather explicitly specify reader = ... or format = ... for the Reader. You can have a look at e.g. pandas syntax for read_parquet, where they explicitly specify engine=....

[marinegor]

Can you do a performance comparison between reading with our PDBParser/Reader vs gemmi?

our PDBParser/Reader is faster, if you just read a file:

In [13]: %timeit u = mda.Universe('/Users/egormarin/Downloads/pdb-mmcif-bench/1FFK.cif')
467 ms ± 9.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [14]: %timeit u = mda.Universe('/Users/egormarin/Downloads/pdb-mmcif-bench/1FFK.pdb')
140 ms ± 1.11 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [15]: %timeit u = mda.Universe('/Users/egormarin/Downloads/pdb-mmcif-bench/1FFK.cif.gz')
480 ms ± 850 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [16]: %timeit u = mda.Universe('/Users/egormarin/Downloads/pdb-mmcif-bench/1FFK.pdb.gz')
168 ms ± 1.57 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

But I feel like there's some constant overhead that I don't know about yet. Perhaps we should add feature for performance optimization of it, what's your opinion?

[orbeckst]

As part of the EOSS4 performance enhancements @richardjgowers looked into fast text file readers. The omnireader in C++ with Cython bindings was a very promising proof of principle (potential for 5-10x speed-up) but we didn't have time to fully package it and use it for the core library. I'll share the final report with you on discord.

[orbeckst]

would not be in favor of neither [] nor | based syntax, and rather explicitly specify reader = ... or format = ... for the Reader. You can have a look at e.g. pandas syntax for read_parquet, where they explicitly specify engine=....

That's fair.

Using an explicit engine=... has the advantage that it's programmatically clearer. I feel that the word "engine" is not very semantically clear although I'll just use it for this discussion. In the context of a Universe we'd then have to use engine for coordinates (format) and topology_engine (for topology_format).

I am pretty sure that explicitly specifying a ReaderClass or a ParserClass is too complicated for most folks as you have to know from where to import it in the first place.

[orbeckst]

Based on your performance comparison, I would not want to default to gemmi for PDB reading, as long as using the PDBReader produces correct results.

[BradyAJohnston]

our PDBParser/Reader is faster, if you just read a file

Which results are which for this? Currently I am just seeing comparisons between PDB / CIF which would have different speeds, and the comparisons with the .gz versions will have the additional overhead of decompressing before reading as well?

[marinegor]

@BradyAJohnston haha you're right, I didn't actually do the comparison😁

here it is:

In [1]: import MDAnalysis as mda

In [2]: %timeit u = mda.Universe('/Users/egormarin/Downloads/pdb-mmcif-bench/1FFK.pdb.gz')
166 ms ± 2.13 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [3]: %timeit u = mda.Universe('/Users/egormarin/Downloads/pdb-mmcif-bench/1FFK.pdb.gz', format='pdb_gemmi')
423 ms ± 9.62 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

it's still slower though.

I'll look into the reasons for it, and after it's speed is on par, we might switch to it fully, but for now I suggest just having this option explicitly for opt-in.

As for this @orbeckst

I am pretty sure that explicitly specifying a ReaderClass or a ParserClass is too complicated for most folks as you have to know from where to import it in the first place.

I think we had this pattern earlier with backend=... for AnalysisBase -- it's type is Literal['dask', 'multiprocessing'] | BackendBase, which covers both "I want to use built-in" and "I want to customize the behavior".

[orbeckst]

If it's a string then it's not a problem; I was thinking about passing a class (which you can do).