Archive for December, 2007

RMSD calculation in R

Sunday, December 23rd, 2007

RMSD stands for root mean square deviation. In structural biology, it is a good way to evaluation the similarity (closeness) of structures. Here is my script in R to calculate RMSD.

rmsd <- function(refX, expX) { # refX and expX should by n*3 matrices. # refX is the coordinates for the reference structure. # expX is the coordinates for the structure needs to be aligned to the reference structure. library(MASS) refo <- apply(refX, 2, mean) refDX <- refX-matrix(rep(refo,dim(refX)[1]), dim(refX)[1], dim(refX)[2], byrow=T) expO <- apply(expX, 2, mean) expDX <- expX-matrix(rep(expO, dim(expX)[1]), dim(expX)[1], dim(expX)[2], byrow=T) A <- matrix(0, 3, 3) for (i in 1:3) { for (j in 1:3) { for (k in 1:dim(refDX)[1]){ A[i,j] <- A[i,j] + refDX[k,i]*expDX[k,j] } } } TR <- sqrt(t(A)*A)*ginv(A) TS <- -TR * expO + refo DX <- refX - expX*TR - matrix(rep(TS,dim(expX)[1]), dim(expX)[1], dim(expX)[2], byrow=T) rmsdv <- sqrt(sum(apply(DX^2,1,sum))/dim(DX)[1]) return(rmsdv) }