I found a function in RDKit (rdMolAlign.AlignMol), but I can't import PDBQT file to RDKit...
I docked a ligand with AutoDock Vina and was generated a PDBQT file with the results of docking and I want to compare this results with the same molecule, present in the PDB file (original ligand).
http://www.rdkit.org/Python_Docs/rdkit.Chem.rdMolAlign-module.html#AlignMol
If you are OK with matrices (I am not so you actually do not have to be either) and some programming (e.g., Python/Numpy), you can calculate RMSD between ligand conformers A and B without any GUI at all.
---
If the protein structures have identical coordinates, you can calculate RMSD between ligand positions A and B directly from the distances between x, y, z coordinates of equivalent atoms in A and B. Distance between two points defined by x, y, z coordinates is calculated as a square root of sum of squares of dx, dy, dz.
If the protein structures have different coordinates, align the protein-ligand complexes using the protein coordinates, and then calculate RMSD between the ligand positions A and B.
---
If you are interested in the difference in conformation of A and B, rather than in the relative orientation of A and B in the binding site of the protein, you will need to align B to A, which can be done as follows:
1. Calculate centroids for A and for B.
2. Calculate rotation matrix.
3. Rotate the matrix.
4. Translate (move B onto A).
5. Calculate RMSD.
Let M1 and M2 be the matrices of x, y, z coordinates of A and B centroids, respectively. Using Python . . .
import numpy
Covariance matrix can be calculated . .
C = numpy.dot(numpy.transpose(M2, M1)
and decomposed . . .
V, S, W = numpy.linalg.svd(C)
d = (numpy.linalg.det(V) * numpy.linalg.det(W)) > 0.0
if d:
S[-1] = -S[-1]
V[:, -1] = -V[:, -1]
and rotation matrix can be calculated . . .
U = numpy.dot(V, W)
and B can be rotated . . .
M2 = numpy.dot(M2, U)
and translated (by removing the average dx, dy, dz between A and the rotated B from x, y, z of the rotated B).
And RMSD can be calculated using x, y, z coordinates of A and the rotated/translated B.
You can use VMD to do this. First, add the non-polar hydrogens to your PDBQT file then use VMD to align the crystallographic complex structure to your docked complex structure and calculate the RMSD.
I never used RDKit before, but if its only issue is that it doesn't accept PDBQT files, then you should add the non-polar hydrogens to your docked complex and delete the info in the PDB after the x,y,z coordinates of each atom.
Hope this helps!
I suggest pymol, where you can directly read pdbqt file of ligand conformers generated through Vina docking run. Choose the best conformer. Also import the original pdb of protein with cocrystallized ligand. You can manually pick the atoms from docked & cocrystallized ligand and through pymol wizard you can measure RMSD directly.
I want to know from Martin who suggested calculation of RMSD using python that is it possible to calculate it by using R or Rstudio and which package is required for it? Thanks Martin any way I am going to try your suggested method using R manually.
Regards
You can use chimera or fconv probably. Both are free.
http://pc1664.pharmazie.uni-marburg.de/drugscore/fconv_download.php
Gromacs has a tool to calculate RMSD.
gmx rms -s cocrystalized_ligand.pdb -f docked_ligand.pdb -fit none -o output.xvg
Or obrms d'OpenBabel
"obrms
Usage: obrms [-firstonly]
Computes the heavy-atom RMSD of identical compound structures.
Structures in multi-structure files are compared one-by-one unless -firstonly
is passed, in which case only the first structure in the reference file is used."
- Badges
- Science topic
- Similar topics
- Filing