How to compute the rmsd between the docked pose and crystal pose?

[leaves520]

Hi, I encounter some error on spyrmsd.rmsd.hrmsd when between the docked pose and crystal pose. I suspect the docked pose from vina didnt have the same atom number as crystal pose. how to deal with this problem? Looking for your reply. thx.

[RMeli]

Hi @leaves520. Which error message are you getting exactly? IIRC Vina removes non-polar hydrogen atoms from the output. Therefore, the crystal pose and the docked pose usually have a different number of atoms. You probably want to compute the heavy-atom RMSD anyways.

If you are working with spyrmsd.molecule.Molecules, then you can use the strip() method to remove hydrogen atoms, and compute the heavy-atom RMSD. See the tutorial for an example. Otherwise, you can simply remove the hydrogen atom from both molecules using your method of choice, before giving them to spyrmsd.

This stripping is performed automatically by spyrmsd.rmsd.rmsdwrapper (see the strip=True argument).

However, the spyrmsd.rmsd.rmsdwrapper only exposes functionality from spyrmsd.rmsd.rmsd (standard RMSD) and from spyrmsd.rmsd.symmrmsd (symmetry-corrected RMSD using graph isomorphism). spyrmsd.rmsd.hrmsd is not exposed, because as we discussed in the accompanying paper the Hungarian algorithm does not respect the molecular graph and can lead to artificially low RMSD values. The spyrmsd.rmsd.hrmsd is only provided for comparison, and should not be used IMHO.

[leaves520]

Thx for your reply. the error message is here : assert atomicn1.shape == atomicn2.shape. The reason is just as you said, the crystal pose and the docked pose usually have a different number of atoms. However, I use the strip() and io.loadmol() refer to tutorial but error still arise.

However, I use open babel to remove the hydrogen atoms of molecules before giving them to spyrmsd，it will be succeed.

[RMeli]

I use open babel to remove the hydrogen atoms of molecules before giving them to spyrmsd，it will be succeed.

Great, I'm glad we figured out the issue and that you found a workaround.

However, I use the strip() and io.loadmol() refer to tutorial but error still arise.

Can you please provide a code snippet? It's odd that this problem persist when using strip(), and some code reproducing the issue might be helpful to track down potential bugs.

[leaves520]

Can you please provide a code snippet? It's odd that this problem persist when using strip(), and some code reproducing the issue might be helpful to track down potential bugs.

hi, RMeli.
I am sorry for my mistake. I carefully check my code and debug. Actually, the main cause of the problem is in rdkit.Chem.MolFromPDBFile rather than strip(). The output is None when use the rdkit.Chem.MolFromPDBFile to load the pdb or mol2 file, the error message is here: Explicit valence for atom # 29 C, 6, is greater than permitted. This output will cause the error in strip(). However, pdb file will successfully loaded by openbabel: next(pybel.readfile(fmt, fname)). Anyway, I appreciate your response.
By the way, based on your experience, would it be better to use openbabel rather than rdkit to read molecules ? Looking forward to your reply, thx.

[RMeli]

Hi @leaves520, thank you for coming back to me with more details. I'm glad to hear the issue was not within spyrmsd. Please do not hesitate to open other issues/discussions if you encounter further problems.

By the way, based on your experience, would it be better to use openbabel rather than rdkit to read molecules?

In my experience, rdkit sanitisation is much more strict than openbabel's by default. Therefore, I started preferring rdkit over openbabel for this reason. (openbabel remains nonetheless very useful at times.)

Looking at you error message, there might be something very wrong in your input file since rdkit detects an explicit valence of 6 on a carbon atom. Instead of switching to openbabel to hide the issue, I would look carefully into the input file to see where the issue actually comes from (and if your structure is sensible).

[leaves520]

Hi, @RMeli
Thank you for your sincere response in this problem. I got a lot of helpful information from you.