How to calculate RMSD and binding energy from the already docked ligand-receptor

Question

I have the following PDB file that looks like this. It is downloadable here. This PDB already contains ligand and receptor together.

The red molecule is the ligand and the yellow is the receptor. They're already docked with ClusPro. My question is how to calculate the RMSD and binding energy from it? Any existing tools (with or without command line) that can help?

I don't have any data to try it out at hand but have a look at the documentation I linked. There you can find ways of constructing an Ensemble which can then be superposed and RMSDs calculated, even pairwise.

Reference is the red one (receptor). I have hundreds of other complexed ligands. I want to calculate how good are those complexes. What's the reasonable measure to do that?

Example code using Prody? please :-)

Answer 1

NB. This is a long comment as opposed to an answer which requires a colossal working chunk of code.

Control

Science is all about controls failing. A great control for docking is the barnase-barnstar complex (PDB:1BRS): both small, but form a strong interactions. If this is done, the RMSD against the original structure is very meaningful and useful.

RMSF

A RMSD is a metric of how much two structures differ. In an ensemble, the metric of how much a structure and the ensemble average is called a RMSF (root-mean square fluctuation). Here the RMSF is sought as the RMSD from first/top pose or worse the manually thrown together pre-docking pose does not mean much.

Clusters

The decoys (PBDs trails) need to be clustered into different ensembles first, before a RMSF can be calculated.

∆∆G and water

For the calculation of difference in Gibbs free energy between the docked minus the docked poses, the structures need to be energy minimised with the forcefield scorefunction that is to be used for scoring. Water is very important for protein-protein interactions, but adds a lot of complication. Say Gromacs is used, gmx energy will return the single state (conformer) Gibbs free energy, but it requires a water box. For the difference, the water box of undocked needs to be the same size and preferably with the same amount of molecules. A much simpler solution is using an implicit solvent model, e.g. Rosetta. As a bonus, the energy minimisation (relax) is relative fast. And specific explicit waters can be added. However, calculating ∆∆G for every decoy will take hundreds of CPU-hours and is a bit nonsensical.

A better approach is to look at residue conservation (consurf), cross-linking mass spec data if available, and if human post-translational modifications (phosphosite) and healthy variants (gnomAD). But ultimately wet lab validation by ITC of different mutants/labelling based upon the docking results etc. is required.