Create a gist now

Instantly share code, notes, and snippets.

@rajarshi /rp-dude.R
Last active Feb 16, 2016

Embed
What would you like to do?
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())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment