PDBx/mmCIF Reader/Topology Reader

[JoaoRodrigues]

Hi all,

I've been using OpenMM to run simulations for the past couple of years. Unlike other tools, OpenMM does not write a 'topology' file (e.g. .gro, .prmtop) but instead creates a Topology object on the fly when loading structures. Because simulation systems tend to be quite large, I've take to use PDBx/mmCIF files as my default file format for writing structures that I use as topologies.

I'd like to start using MDAnalysis but right now this involves jumping through a bunch of hoops to get my topologies in a format that is parseable. It'd be much easier if I could just load a PDBx/mmCIF file as a topology, specially since it's now the default file format for structures in the PDB.

To this end, I've started working on writing a simple PDBxParser class. I wouldn't mind extending it to a PDBxReader/Writer class but that is not necessarily a priority (specially the writer). Would this be an interesting feature to add in your opinion? See the code here, I modeled the class after PDBParser.

Thanks for the great work with the library so far!

EDIT: Since I cannot label issues, I'm just editing the title of the issue for now to make it clear!

[orbeckst]

Hi @JoaoRodrigues , new formats are a great addition. Can you just create a PR and then we can comment directly on the code, run tests, etc?

[JoaoRodrigues]

Sounds good. Will do. Thanks!

[orbeckst]

Pure python implementation: https://github.com/Electrostatics/mmcif_pdbx

[orbeckst]

Also, chemfiles reads mmCIF and can be used inside MDAnalysis, see Reading trajectories with chemfiles.

[orbeckst]

Note on trying to use chemfiles:

With the 2.0.0b,

u = mda.Universe("1ake.cif.gz", topology_format="chemfiles")

fails with

'CHEMFILES' isn't a valid topology format, nor a coordinate format

I can't just read it with

u = mda.Universe("1ake.cif.gz", format="chemfiles")

either as it gives the same ValueError.

Just trying something stupid with the topology as a PDB also fails

u = mda.Universe("1ake.pdb.gz", "1ake.cif.gz", format="chemfiles")

with

ChemfilesError: can not read file '1ake.cif.gz' at step 0, it does not contain any step.

Conclusion: I couldn't get it to work with chemfiles. Maybe @Luthaf has some ideas but a key problem seems to be that our (EDIT) chemfiles converter do not work for topologies —— ??

[Luthaf]

I think that topology reading is not registered with MDA, since the chemfiles adapter is implemented as a coordinate reader/writer. Although the conversion from chemfiles to MDA topologies is already implemented, so it should be mostly a question of adding a new ChemfilesParser class. I'll have a look at this next week!

ChemfilesError: can not read file '1ake.cif.gz' at step 0, it does not contain any step.

This is a bit strange and probably a bug, is this 1ake from the wwwPDB?

[Luthaf]

This is a bit strange and probably a bug, is this 1ake from the wwwPDB?

Ok, I understand this part. The core of the issue is that to relatively different formats want to use the same extension: mmCIF and crystallography CIF. While they both use the same STAR format, they specify data in different ways. Unfortunately, chemfiles is associating the .cif extension with crystallography CIF files, and uses .mmcif for mmCIF.

There is a simple workaround though, since you can specify the format to use manually, with something like

u = mda.Universe("1ake.pdb.gz", "1ake.cif.gz", format="chemfiles", chemfiles_format="mmCIF / GZ")

Unfortunately, this still fails in the case of 1AKE (but should work for other files). I'll fix the 1AKE issue, it should be working in the next patch release.

I would also like to introduce a better format guessing functionality, to decide between mmCIF and crystallography CIF on the fly instead of having the use specify it manually.