rmsd-helper.R

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 = ',')))
}