How to calculate spyrmsd in mdanalysis

[ldx022]

Hello!
At present, I align and calculate rmsd in the following ways in mdanalysis just as the tutorial says

R = rms.RMSD(md,  # coordinates to align
             ref,  # reference coordinates
             select = '(protein and not name H*) ', #group to superimpose and calculate RMSD
             groupselections= ['(resname LIG)']) #groups only for RMSD
R.run(start=0, stop=1000, step=1)

My aim is to align the protein first and then calculate the rmsd of the ligand. I was curious if spyrmsd could be added to calculate the rmsd of the ligand.
Thank you!

[davidoskky]

Hello @ldx022,
I'm afraid spyrmsd cannot be applied directly to mdanalysis.
You can use the mdanalysis rmsd function to calculate the rmsd directly over the whole trajectory, though I do not remember whether that one does consider symmetry or not.
I didn't understand if you wish to calculate the rmsd between the unaligned trajectory and the aligned trajectory or between the frames of the aligned trajectory, in any case you can do bot: doc.
If you wish to use spyrmsd, then the most straightforward approach would be to select the ligand in the trajectory, export it to some format readable by openbabel or rdkit — namely PDB — and use that to import the molecule in spyrmsd. I would not recommend this option because by passing through the PDB format, you may introduce some errors.
The last option would be to manually build spyrmsd Molecule objects using the mdanalysis data. To do this, you will need to select the molecule in mdanalysis and extract the atom types as well as the coordinates of all atoms which are part of the ligand. Moreover, you will need to build the connectivity matrix; I do not believe mdanalysis provides an internal representation of the matrix itself nor that it has any built in method to generate it. If you want to go this way, all you need should be explained here, the procedure may be easier if you then convert your selection (making sure you only selected a fully connected ligand) to a TopologyGroup which provides some helper functions such as bonds() which could make your life easier.

I hope this may help you.

[RMeli]

Thank you very much for taking the time for an extensive reply, it's really appreciated!

The last option would be to manually build spyrmsd Molecule objects using the mdanalysis data.

I just want to point out that there is no need to go through a spyrmsd.molecule.Molecule object. One can pass information extracted from MDAnalysis to a lower-level spyrmsd function such as symmrmsd.

The spyrmsd.molecule.Molecule object is mainly there for I/O convenience and for easily interacting with RDKit/OpenBabel. Since MDAnalysis deals directly with NumPy arrays (as does spyrmsd), they can be used directly.

But you are totally right that the relevant information needs to be manually extracted from MDAnalysis data.

[RMeli]

Hello! spyrmsd can certainly interoperate with MDAnalysis and compute symmetry-corrected RMSD over a trajectory. MDAnalysis does not use symmetry correction.

spyrmsd only needs the following information from your trajectory, as NumPy arrays: adjacency matrix, coordinates, and atomic properties (usually elements or atom types). This information can be extracted from an MDAnalysis trajectory and fed directly into spyrmsd low-level functions.

As @davidoskky already pointed out, the adjacency matrix can be constructed from bonding information but the procedure is a bit convoluted because MDAnalysis does not have a function to build the connectivity matrix.

I wanted to integrate spyrmsd into MDAnalysis, but I went for a MDAKit instead: you can have a look at the spyrmsdkit. Unfortunately I did not have enough time to finish off that project, but it should be usable (please use with care, test it extensively, and report bugs!) or at least give you a solid idea of how to extract the information you need from MDAnalysis.

---

Relevant code from the MDAKit (please be aware it is still under GPL-v2 but will be re-licensed when MDAnslysis re-licensing will be finished):

• Compute adjacency matrix (vectorized implementation)
• Compute symmetry-corrected RMSD

With spyrmsdkit, computing symmetry-corrected RMSDs on an MD trajectory should be as easy as

# This import will likely change in the future to something simpler...
from spyrmsdkit.analysis.spyrmsd import SPyRMSD

R = SPyRMSD(
  u.select_atoms("resname LIG"),
)
R.run()

SPyRMSD inherit from mdanalysis.analysis.base.AnalysisBase and has a simpler API than the standard MDAnalysis RMSD class.

As mentioned above, there are test in the MDAKit, but it's not really fully battle-tested (i.e. I haven't used it on a real project yet), so there might be bugs and sharp edges.

[ldx022]

Thank you for your prompt response. I appreciate your suggestions and assistance. As I am new to MDAnalysis, I'm still getting familiar with its code logic. I will take time to thoroughly understand your provided advice when possible.

However, I'm currently in a bit of a hurry. I have an existing MDAnalysis conda environment where I've installed spyrmsdkit for testing. Following your advice, I tried the following script module:

R = SPyRMSD(md,  # coordinates to align
            ref,  # reference coordinates
            select='(protein and not name H*) or (resname XMA)', # group to superimpose and calculate RMSD
            groupselections=['(resname LIG)']) # groups only for RMSD
R.run(start=0, stop=1000, step=1)

Unfortunately, I encountered this error:

R = SPyRMSD(md,  # coordinates to align
   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/export/home/ldx022/software/spyrmsdkit-main/spyrmsdkit/analysis/spyrmsd.py", line 87, in __init__
   if len(self.ref_atoms) != len(self.mobile_atoms):
      ^^^^^^^^^^^^^^^^^^^
TypeError: object of type 'Universe' has no len()

Interestingly, using rms.RMSD instead of SPyRMSD doesn't cause any issues. Could you please help me understand why this might be happening?

By the way, for aligning protein residues, I believe whether to use spyrmsd or not doesn't significantly impact the outcome. My main concern is about calculating the RMSD of the ligand in the trajectory against the reference ligand using spyrmsd. Specifically, I've already aligned the trajectory, so how do I directly compute the SPyRMSD for the ligands in this trajectory? Any further advice or suggestions would be greatly appreciated.

Thank you once again for your support.

[RMeli]

Great you are trying out spyrmsdkit. As I mentioned above, please make sure the results are sensible, since the project is a bit unfinished and not extensively battle-tested (unlike spyrmsd itself).

The SPyRMSD analysis class, which is unfortunately not well documented (see the SPyRMSD class docstring), works differently from MDAnslysis RMSD class (because the use case is narrower). You need to pass atom groups directly, not universes.

Once you have your trajectory pre-processed (wrapped, aligned, ...) the way you want, you could do something like the following:

lig_sel = "resname LIG"

R = SPyRMSD(md.select_atoms(lig_sel),
            ref.select_atoms(lig_sel)
)
R.run(start=0, stop=1000, step=1)

The SPyRMSD class does not do superimposition of the protein structure.

---

Interestingly, using rms.RMSD instead of SPyRMSD doesn't cause any issues. Could you please help me understand why this might be happening?

As mentioned above, the SPyRMSD analysis class is different from the RMSD analysis class, so they should be used differently. Unfortunately the former still lacks proper documentation.

By the way, for aligning protein residues, I believe whether to use spyrmsd or not doesn't significantly impact the outcome. My main concern is about calculating the RMSD of the ligand in the trajectory against the reference ligand using spyrmsd. Specifically, I've already aligned the trajectory, so how do I directly compute the SPyRMSD for the ligands in this trajectory?

Aligning the protein is usually done using backbone atoms only anyways, so symmetry correction would not have an impact (and symmetry correction is expensive, so it would take a long time for a whole protein). Your approach is what I would do: align the protein with standard RMSD, and then compute the symmetry-corrected RMSD for the ligand.

[ldx022]

Hello, thank you for your patient response!

Following your suggestions, I implemented the process again and realized that bonding information is needed. Therefore, the topology file must be a .mol2 file. This seems to make direct selection and calculation on proteins difficult, so I had to first align the protein and then extract the ligand's trajectory part:

md = mda.Universe("./lig_H.mol2", 'lig.xtc') 
ref = mda.Universe("./lig_H.mol2","./LIG.pdb")
LIG = 'resname LIG'
R = SPyRMSD(md.select_atoms('resname LIG and not name H*'), 
             ref.select_atoms('resname LIG and not name H*'),  
)
R.run(start=0, stop=1000, step=1)

However, I found that using '(resname LIG) and (not name H*)' to exclude H atoms causes an error:

    R.run(start=0, stop=1000, step=1)
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/export/home/ldx022/software/miniconda3/envs/spyrmsdkitb/python3.12/site-packages/MDAnalysis/analysis/base.py", line 436, in run
    self._prepare()
  File "/export/home/ldx022/software/spyrmsdkit-main/spyrmsdkit/analysis/spyrmsd.py", line 112, in _prepare
    self.ref_adj = utils.adjacency_matrix(self.ref_atoms)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/export/home/ldx022/software/spyrmsdkit-main/spyrmsdkit/analysis/utils.py", line 36, in adjacency_matrix
    A[indices[:, 0], indices[:, 1]] = 1
    ~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IndexError: index 28 is out of bounds for axis 1 with size 28

Index 28 is precisely the first H atom I excluded. But when I changed it to 'resname LIG', it ran successfully. I suspect this issue is caused by the hydrogen atoms. Should I create a separate .mol2 file for the ligand without H atoms, or is there a more convenient alternative?

Thank you for your suggestions!

[ldx022]

Alright, I had a look at the reproducer of the error you reported. The issue is that when creating an AtomGroup in MDAnalysis with select_atoms, the bonds attribute also contains dangling bonds. That is, bonds from atoms of the AtomGroup to atoms that are not part of the AtomGroup. This happens when stripping hydrogen atoms from the ligand (since they are bonded to your selection).
I added a filtering step to only consider bonds within atoms of the AtomGroup when building the adjacency matrix, and the error went away on my side.

Yes! Just like my previous approach, I found that in the ligand PDB files which include hydrogen atoms, the CONNECT part also contains some hydrogen atoms. Therefore, it's necessary to regenerate a PDB file for the ligand without hydrogen atoms, to exclude the hydrogen atom serial numbers from the CONNECT field.

Would you be able to test the strip-hydrogens-issue branch in the spyrmsdkit repo, which contains this filtering? Your feedback would be extremely valuable, before merging it.

Sure. May I ask if I need to create a new environment and reinstall SPyRMSDkit? Or is there a simpler way to retest in my original SPyRMSDkit environment?

[RMeli]

May I ask if I need to create a new environment and reinstall SPyRMSDkit? Or is there a simpler way to retest in my original SPyRMSDkit environment?

Assuming you cloned my repo and installed spyrmsdkit with pip, the easiest way is to do the following within your current Python environment:

Go to the cloned spyrmsdkit repo and switch to the strip-hydrogens-issue branch

git fetch
git checkout strip-hydrogens-issue

If you installed spyrmsdkit with pip install -e . (-e is for "editable" mode), there should be nothing else to do. If you installed spyrmsdkit with pip install . then you need to uninstall it (pip uninstall spyrmsdkit) and re-install it (pip install ./pip install -e .).

[ldx022]

Hello! Since I don't have root access on the server and also lack the git command, I usually download the code from GitHub in zip format. Then, I create a conda environment to install it. However, when I tried to perform the same installation operation on the "strip-hydrogens-issue" branch, it failed, resulting in the following error:

[ldx022@login02 spyrmsdkit-strip-hydrogens-issue]$ conda create --name spyrmsdkit1
(spyrmsdkit1) [ldx022@login02 spyrmsdkit-strip-hydrogens-issue]$ conda env update --name spyrmsdkit --file devtools/conda-envs/test_env.yaml
(spyrmsdkit1) [ldx022@login02 spyrmsdkit-strip-hydrogens-issue]$ conda env update --name spyrmsdkit --file docs/requirements.yaml
(spyrmsdkit1) [ldx022@login02 spyrmsdkit-strip-hydrogens-issue]$ /export/home/ldx022/software/miniconda3/envs/python39/bin/pip install -e .
Obtaining file:///export/home/ldx022/software/spyrmsdkit-strip-hydrogens-issue
  Installing build dependencies ... done
  Checking if build backend supports build_editable ... done
  Getting requirements to build editable ... error
  error: subprocess-exited-with-error
  
  × Getting requirements to build editable did not run successfully.
  │ exit code: 1
  ╰─> [27 lines of output]
      Traceback (most recent call last):
        File "/export/home/ldx022/software/miniconda3/envs/python39/lib/python3.9/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 353, in <module>
          main()
        File "/export/home/ldx022/software/miniconda3/envs/python39/lib/python3.9/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 335, in main
          json_out['return_val'] = hook(**hook_input['kwargs'])
        File "/export/home/ldx022/software/miniconda3/envs/python39/lib/python3.9/site-packages/pip/_vendor/pyproject_hooks/_in_process/_in_process.py", line 132, in get_requires_for_build_editable
          return hook(config_settings)
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/build_meta.py", line 441, in get_requires_for_build_editable
          return self.get_requires_for_build_wheel(config_settings)
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/build_meta.py", line 325, in get_requires_for_build_wheel
          return self._get_build_requires(config_settings, requirements=['wheel'])
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/build_meta.py", line 295, in _get_build_requires
          self.run_setup()
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/build_meta.py", line 311, in run_setup
          exec(code, locals())
        File "<string>", line 26, in <module>
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/__init__.py", line 103, in setup
          return distutils.core.setup(**attrs)
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/_distutils/core.py", line 147, in setup
          _setup_distribution = dist = klass(attrs)
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/dist.py", line 314, in __init__
          self.metadata.version = self._normalize_version(self.metadata.version)
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/dist.py", line 350, in _normalize_version
          normalized = str(Version(version))
        File "/tmp/pip-build-env-zaotp68v/overlay/lib/python3.9/site-packages/setuptools/_vendor/packaging/version.py", line 198, in __init__
          raise InvalidVersion(f"Invalid version: '{version}'")
      setuptools.extern.packaging.version.InvalidVersion: Invalid version: 'refs/pull/1/head'
      [end of output]
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
error: subprocess-exited-with-error

× Getting requirements to build editable did not run successfully.
│ exit code: 1
╰─> See above for output.

note: This error originates from a subprocess, and is likely not a problem with pip.

[RMeli]

Mmmh... I've never seen this error before.

setuptools.extern.packaging.version.InvalidVersion: Invalid version: 'refs/pull/1/head'

If I git clone the repository and then switch to the strip-hydrogens-issue branch and then python -m pip install ., everything works as expected. However, I can confirm that trying to python -m pip install . from the archived downloaded with git results in the same error you have.

I'll look into it (I never download archives from GitHub), but in the meantime a workaround could be to git clone and git checkout strip-hydrogens-issue locally and then upload the repository to your server?

[RMeli]

@ldx022 I looked into this (sorry for the delay) and the reason for the error is that the version is determined via git. I'll try to figure out if this is a wider issue with MDAKits in general and open an issue there if it's the case.

[RMeli]

@ldx022, thank you very much for testing out spyrmsdkit and providing feedback!

I just merged RMeli/spyrmsdkit#1, which should solve the issues you were seeing, and make spyrmsdkit usable to compute small-molecule RMSDs over a (pre-aligned) MDAnalysis trajectory.

If you end up using it and encounter other problems, don't hesitate to ask questions, or open an issue in the spyrmsdkit repository. I will try to find some time to add documentation and add it to the MDAKit registry.

[ldx022]

OK,Thank you for your help!