RMSD between re-docked complex and co-crystallized complex

Question

[ksraji251]

Hello,

I would like to find the rmsd between re-docked complex and co-crystallized complex for docking validation? What are the tools available for this purpose? I tried with pymol using align command, but I got the rmsd between two proteins not for complexes.

Answer 1

[jgreener]

You could try using MDAnalysis, see for example https://userguide.mdanalysis.org/stable/examples/analysis/alignment_and_rms/rmsd.html.

Answer 2

[Michael]

I remember doing this using the rmsd function in UCSF Chimera in this paper. To achieve that, I loaded the docking result and the original structure from PDB (that from which the ligand was deleted prior to docking). Then I deleted all protein residues for simplification. Then, RMSD was computed only taking C, N, and O into account, but to do this, all atoms needed to be aligned (as to appear as proper pairs) between the docking and the ligand. I think I did this manually all the time by superimposing the ligand structures and entering the atom identifiers. I feel there must be a better way but it worked for the few structures we had.

Edit: So, I found the helper script in R. To use that, one has to create a new PDB file for the purpose of aligning the atoms, don't overwrite your docking file. There, the compounds are superimposed such that all atoms are nearest neighbor of their respective image in the other ligand. The PDB should exactly contain two models, the docking ligand and the PDB ligand.

Then function align.nn.pdb gives the respective atom ids. The function name is a bit misleading because you would have done the alignment yourself already before.

Then, formatAs.chimeraAtomIds can convert the output into a string of correctly paired atom identifiers to use with RMSD.

If you think this is crude, I totally agree and please provide a better solution.

RMSD was sensitive for even slight rotation of the docking pose, so I think distance between centroids is also a good measure to provide.

Here is the Gist (pls, don't ask for details on how this even works :) )

	library(Rpdb)
	library(pdist)
	#mypdb <- read.pdb("SCARB_S361384_results_itasser/1dve-rmsd.pdb", HETATM = T, ATOM = T, MODEL = 1:2)
	my.pdist <- function(x,y) {
	ret <- as.matrix(pdist(x,y), nrow = nrow(x), ncol=nrow(y))
	dimnames(ret) <- list(row.names(x), row.names(y))
	return(list(ret))
	}
	dist.centres <- function(mypdb) {
	stopifnot( inherits(mypdb$MODEL.1, 'pdb'), inherits(mypdb$MODEL.2, 'pdb'))
	c1 <- as.numeric(centres(mypdb$MODEL.1))
	c2 <- as.numeric(centres(mypdb$MODEL.2))
	dist(rbind(c1,c2))
	}
	align.nn.pdb <- function(mypdb, include.atoms=c('C','N','O')) {
	stopifnot( inherits(mypdb$MODEL.1, 'pdb'), inherits(mypdb$MODEL.2, 'pdb'),
	nrow(mypdb$MODEL.1$atoms) == nrow(mypdb$MODEL.2$atoms))
	m1.elem <- substr(mypdb$MODEL.1$atoms$elename,1,1) # bit crude, missing two letter elements, but ok for now
	m2.elem <- substr(mypdb$MODEL.2$atoms$elename,1,1)
	## split the coordinate matrix by elements
	m1.coords <- as.matrix(coords(mypdb$MODEL.1))
	row.names(m1.coords) <- 1:nrow(m1.coords)
	m2.coords <- as.matrix(coords(mypdb$MODEL.2))
	row.names(m2.coords) <- 1:nrow(m2.coords)
	m1.cbe <- split.data.frame(m1.coords, m1.elem)
	m2.cbe <- split.data.frame(m2.coords, m2.elem)
	atom.dists <- mapply(m1.cbe[include.atoms], m2.cbe[include.atoms], FUN=my.pdist)
	duplist <- unlist(lapply(lapply(atom.dists, apply, 1, which.min), duplicated))
	if (any(duplist)) {
	warning(paste("Some next-neighbors were duplicated:\n", which(duplist)))
	}
	atom.nn.id <- (sapply(atom.dists, apply, 1, function(x) {as.numeric(names(x[which.min(x)]))}))
	atom.nn.name <- sapply(atom.nn.id, function(x) {cbind(mypdb$MODEL.1$atoms$elename[as.numeric(names(x))],mypdb$MODEL.2$atoms$elename[x])})
	#print (atom.nn.name)
	sqd <- unlist(sapply(atom.dists, apply, 1, min))^2
	rmsd <- sqrt(
	sum(sqd/length(sqd))
	)
	ret <- list(atom.dists=atom.dists, atom.nn.id=atom.nn.id, atom.nn.name=atom.nn.name, RMSD=rmsd, distance.centres=dist.centres(mypdb))
	class (ret) <- 'atom.nn.structure'
	return (ret)
	}
	formatAs.chimeraAtomIds <- function(o) {
	stopifnot(inherits(o, 'atom.nn.structure'))
	s1 <- unlist(sapply(o$atom.nn.name, '[',,1))
	s2 <- unlist(sapply(o$atom.nn.name, '[',,2))
	return(paste0('#0.1@',paste(s1, sep=',', collapse = ','),' #0.2@' ,paste(s2, sep=',', collapse = ',')))
	}

view raw rmsd-helper.R hosted with ❤ by GitHub

Here is an example PDB file to try this on. One model contains ATOM's (the docked ligand) and the other (the original PDB model) HETATOM definitions, but I think that doesn't matter I think.

	MODEL 1
	ATOM 1 C1 UNK 1 46.878 21.312-121.026 1.00 0.00 C
	ATOM 2 C2 UNK 1 46.008 21.716-122.189 1.00 0.00 C
	ATOM 3 C3 UNK 1 46.381 22.574-123.243 1.00 0.00 C
	ATOM 4 C4 UNK 1 47.629 23.170-123.351 1.00 0.00 C
	ATOM 5 C5 UNK 1 47.990 24.006-124.415 1.00 0.00 C
	ATOM 6 N1 UNK 1 47.231 24.376-125.447 1.00 0.00 N
	ATOM 7 C6 UNK 1 47.918 25.174-126.264 1.00 0.00 C
	ATOM 8 C7 UNK 1 47.497 25.776-127.450 1.00 0.00 C
	ATOM 9 C8 UNK 1 46.216 25.625-127.970 1.00 0.00 C
	ATOM 10 C9 UNK 1 45.777 26.245-129.151 1.00 0.00 C
	ATOM 11 C10 UNK 1 44.459 25.876-129.359 1.00 0.00 C
	ATOM 12 C11 UNK 1 44.086 25.014-128.291 1.00 0.00 C
	ATOM 13 N2 UNK 1 45.175 24.882-127.460 1.00 0.00 N
	ATOM 14 H1 UNK 1 45.203 24.352-126.649 1.00 0.00 H
	ATOM 15 C12 UNK 1 42.833 24.441-128.161 1.00 0.00 C
	ATOM 16 C13 UNK 1 42.487 23.582-127.110 1.00 0.00 C
	ATOM 17 N3 UNK 1 43.249 23.188-126.090 1.00 0.00 N
	ATOM 18 C14 UNK 1 42.565 22.366-125.294 1.00 0.00 C
	ATOM 19 C15 UNK 1 43.006 21.743-124.120 1.00 0.00 C
	ATOM 20 C16 UNK 1 44.275 21.925-123.594 1.00 0.00 C
	ATOM 21 C17 UNK 1 44.729 21.314-122.408 1.00 0.00 C
	ATOM 22 N4 UNK 1 45.306 22.690-124.095 1.00 0.00 N
	ATOM 23 H2 UNK 1 45.277 23.220-124.906 1.00 0.00 H
	ATOM 24 C21 UNK 1 41.214 22.204-125.848 1.00 0.00 C
	ATOM 25 C22 UNK 1 41.165 22.956-126.969 1.00 0.00 C
	ATOM 26 C23 UNK 1 39.987 23.119-127.894 1.00 0.00 C
	ATOM 27 C29 UNK 1 46.595 27.153-130.033 1.00 0.00 C
	ATOM 28 C30 UNK 1 49.273 25.336-125.685 1.00 0.00 C
	ATOM 29 C31 UNK 1 49.277 24.597-124.537 1.00 0.00 C
	ATOM 30 C32 UNK 1 50.430 24.437-123.579 1.00 0.00 C
	ATOM 31 H33 UNK 1 44.170 23.459-125.948 1.00 0.00 H
	ATOM 32 H34 UNK 1 46.310 24.103-125.584 1.00 0.00 H
	ATOM 33 C18 UNK 1 43.934 20.379-121.532 1.00 0.00 C
	ATOM 34 C19 UNK 1 42.874 21.176-120.767 1.00 0.00 C
	ATOM 35 C20 UNK 1 41.875 20.225-120.158 1.00 0.00 C
	ATOM 36 O1 UNK 1 42.064 19.774-119.043 1.00 0.00 O
	ATOM 37 O2 UNK 1 40.878 19.908-120.784 1.00 0.00 O
	ATOM 38 C24 UNK 1 40.101 21.363-125.277 1.00 0.00 C
	ATOM 39 C25 UNK 1 40.288 19.907-125.711 1.00 0.00 C
	ATOM 40 C26 UNK 1 40.427 19.034-124.491 1.00 0.00 C
	ATOM 41 O3 UNK 1 39.432 18.632-123.913 1.00 0.00 O
	ATOM 42 O4 UNK 1 41.535 18.730-124.081 1.00 0.00 O
	ATOM 43 C33 UNK 1 50.385 26.128-126.236 1.00 0.00 C
	ATOM 44 C34 UNK 1 51.297 25.548-127.008 1.00 0.00 C
	ATOM 45 C27 UNK 1 43.601 26.298-130.480 1.00 0.00 C
	ATOM 46 C28 UNK 1 43.602 25.608-131.615 1.00 0.00 C
	ENDMDL
	MODEL 2
	TER 1743 GLU A 223
	HETATM 1744 FE HEM A 300 45.050 24.022-125.819 1.00 31.66 Fe
	HETATM 1745 CHA HEM A 300 42.970 21.736-124.203 1.00 31.51 C
	HETATM 1746 CHB HEM A 300 42.628 24.457-128.166 1.00 35.58 C
	HETATM 1747 CHC HEM A 300 47.134 26.103-127.388 1.00 31.88 C
	HETATM 1748 CHD HEM A 300 47.449 23.291-123.503 1.00 29.47 C
	HETATM 1749 NA HEM A 300 43.170 23.219-126.118 1.00 33.16 N
	HETATM 1750 C1A HEM A 300 42.460 22.337-125.340 1.00 33.89 C
	HETATM 1751 C2A HEM A 300 41.099 22.118-125.890 1.00 35.06 C
	HETATM 1752 C3A HEM A 300 41.037 22.897-127.021 1.00 34.56 C
	HETATM 1753 C4A HEM A 300 42.302 23.569-127.153 1.00 34.84 C
	HETATM 1754 CMA HEM A 300 39.859 23.050-127.974 1.00 35.13 C
	HETATM 1755 CAA HEM A 300 40.047 21.227-125.291 1.00 36.31 C
	HETATM 1756 CBA HEM A 300 40.172 19.704-125.676 1.00 41.11 C
	HETATM 1757 CGA HEM A 300 39.425 18.791-124.703 1.00 42.62 C
	HETATM 1758 O1A HEM A 300 38.414 19.188-124.088 1.00 42.42 O
	HETATM 1759 O2A HEM A 300 39.870 17.588-124.546 1.00 41.49 O
	HETATM 1760 NB HEM A 300 44.935 25.052-127.437 1.00 34.67 N
	HETATM 1761 C1B HEM A 300 43.839 25.146-128.307 1.00 35.08 C
	HETATM 1762 C2B HEM A 300 44.164 26.101-129.396 1.00 35.59 C
	HETATM 1763 C3B HEM A 300 45.408 26.559-129.170 1.00 35.26 C
	HETATM 1764 C4B HEM A 300 45.868 25.898-127.957 1.00 33.69 C
	HETATM 1765 CMB HEM A 300 43.210 26.445-130.515 1.00 36.93 C
	HETATM 1766 CAB HEM A 300 46.214 27.560-129.981 1.00 35.37 C
	HETATM 1767 CBB HEM A 300 47.150 27.224-130.844 1.00 38.02 C
	HETATM 1768 NC HEM A 300 46.951 24.582-125.492 1.00 28.84 N
	HETATM 1769 C1C HEM A 300 47.653 25.500-126.250 1.00 30.81 C
	HETATM 1770 C2C HEM A 300 48.979 25.721-125.703 1.00 30.20 C
	HETATM 1771 C3C HEM A 300 49.058 24.911-124.613 1.00 29.71 C
	HETATM 1772 C4C HEM A 300 47.783 24.216-124.499 1.00 30.07 C
	HETATM 1773 CMC HEM A 300 50.000 26.667-126.293 1.00 31.46 C
	HETATM 1774 CAC HEM A 300 50.168 24.652-123.612 1.00 30.78 C
	HETATM 1775 CBC HEM A 300 51.372 25.169-123.561 1.00 30.47 C
	HETATM 1776 ND HEM A 300 45.194 22.732-124.152 1.00 29.99 N
	HETATM 1777 C1D HEM A 300 46.269 22.600-123.324 1.00 30.30 C
	HETATM 1778 C2D HEM A 300 46.004 21.676-122.276 1.00 30.93 C
	HETATM 1779 C3D HEM A 300 44.716 21.236-122.479 1.00 30.24 C
	HETATM 1780 C4D HEM A 300 44.217 21.903-123.658 1.00 30.46 C
	HETATM 1781 CMD HEM A 300 46.966 21.277-121.171 1.00 30.78 C
	HETATM 1782 CAD HEM A 300 43.987 20.242-121.624 1.00 32.66 C
	HETATM 1783 CBD HEM A 300 42.922 20.951-120.727 1.00 33.65 C
	HETATM 1784 CGD HEM A 300 42.199 19.973-119.888 1.00 33.06 C
	HETATM 1785 O1D HEM A 300 42.371 19.906-118.651 1.00 34.91 O
	HETATM 1786 O2D HEM A 300 41.409 19.194-120.489 1.00 32.10 O
	ENDMDL
	CONECT 124 1744
	CONECT 1750 1745 1749 1751
	CONECT 1761 1746 1760 1762
	CONECT 1769 1747 1768 1770
	CONECT 1777 1748 1776 1778
	CONECT 1751 1750 1752 1755
	CONECT 1762 1761 1763 1765
	CONECT 1770 1769 1771 1773
	CONECT 1778 1777 1779 1781
	CONECT 1752 1751 1753 1754
	CONECT 1763 1762 1764 1766
	CONECT 1771 1770 1772 1774
	CONECT 1779 1778 1780 1782
	CONECT 1753 1746 1749 1752
	CONECT 1764 1747 1760 1763
	CONECT 1772 1748 1768 1771
	CONECT 1780 1745 1776 1779
	CONECT 1755 1751 1756
	CONECT 1766 1763 1767
	CONECT 1774 1771 1775
	CONECT 1782 1779 1783
	CONECT 1756 1755 1757
	CONECT 1767 1766
	CONECT 1775 1774
	CONECT 1783 1782 1784
	CONECT 1757 1756 1758 1759
	CONECT 1784 1783 1785 1786
	CONECT 1745 1750 1780
	CONECT 1746 1753 1761
	CONECT 1747 1764 1769
	CONECT 1748 1772 1777
	CONECT 1754 1752
	CONECT 1765 1762
	CONECT 1773 1770
	CONECT 1781 1778
	CONECT 1744 1846 124 1749 1760
	CONECT 1744 1768 1776
	CONECT 1749 1750 1753
	CONECT 1749 1744
	CONECT 1760 1761 1764
	CONECT 1760 1744
	CONECT 1768 1769 1772
	CONECT 1768 1744
	CONECT 1776 1777 1780
	CONECT 1776 1744
	CONECT 1758 1757
	CONECT 1785 1784
	CONECT 1759 1757
	CONECT 1786 1784
	CONECT 1846 1744
	END

view raw 1DVE-rmsd-example.pdb hosted with ❤ by GitHub