rp-dude.R

library(parallel)
library(dplyr)
library(ggplot2)

## Fig 6A from http://www.jcheminf.com/content/8/1/1

dudefiles <- list.files(path='dude', pat=glob2rx("*txt"), full=TRUE)
dat <- mclapply(dudefiles, function(f) {
  cat(f, '\n')
  d <- read.table(f, header=TRUE, comment='')
  names(d) <- c('chemblid', 'klass', 'Glide', 'Surflex', 'Flexx', 'GOLD')
  d$GOLD <- -1*d$GOLD ## Since higher scores are better
  d$Surflex <- -1*d$Surflex ## Since higher scores are better (pers. comm.)
  ranks <- apply(d[,-c(1,2)], 2, rank)
  rp <- rowSums(log(ranks))
  tmp <- data.frame(d[,1:2], ranks, lrp=rp)
  tmp <- tmp[order(tmp$lrp),]
  pos.active <- which(tmp$klass == 'ACTIF')

  ## generate pos.active for individual docking scores
  dock.pos.active <- list()
  for (i in 1:4) {
    tmp <- data.frame(d[,1:2], r=ranks[,i])
    tmp <- tmp[order(tmp$r),]
    dpos.active <- which(tmp$klass == 'ACTIF')
    dock.pos.active[[i]] <- dpos.active
  }
  names(dock.pos.active) <- c('Glide', 'Surflex', 'Flexx', 'GOLD')
  
  list(fname=f, target=gsub("dude/", "", strsplit(f,"_")[[1]][1]),
       nr=nrow(d),
       pos.active=pos.active,
       dock.pos.active=dock.pos.active,
       nact=length(which(d$klass=='ACTIF')))
}, mc.cores=5)

ntests <- seq(1, 50, by=1)
f1 <- function(pos, ntest) {
  length(pos[pos<=ntest]) > 0
}

tmp1 <- do.call(rbind, lapply(ntests, function(ntest) {
  tmp <- sapply(dat, function(x) f1(x$pos.active, ntest))
  data.frame(method='RP', ntest=ntest, ts=length(which(tmp))/102)
}))

tmp2 <- lapply(names(dat[[1]]$dock.pos.active), function(dname) {
  do.call(rbind, lapply(ntests, function(ntest) {
    tmp <- sapply(dat, function(x) f1(x$dock.pos.active[[dname]], ntest))
    data.frame(method=dname, ntest=ntest, ts=length(which(tmp))/102)
  }))
})

tmp <- rbind(tmp1, do.call(rbind, tmp2))

ggplot(tmp,aes(ntest, ts*100, group=method, colour=method)) +
  geom_line(size=1)+
  xlab(expression(n[test]))+
  ylab(expression(T[h>0],"(%)"))+
  scale_y_continuous(breaks=seq(0,100,10), limits=c(0,100))+  
  theme(axis.title = element_text(size=14, face='bold'),
        axis.text = element_text(size=10),
        legend.title = element_blank())