R code: Figure 3

Load the following packages:

library(limma)
library(qvalue)
library(RColorBrewer)

Set parameters:

n <- 1394 # number of proteins
d0 <- 4.428296 # prior degrees of freedom
s02 <- 0.03197716 # prior variance

nit <- 1000 # number of iterations
hm <- 100 # number of significant proteins
hm.fc <- 1.25 # fold change for significant proteins

Code case and control groups according to limma-package framework:

design <- model.matrix(~factor(c(rep(1,4), rep(2,4))))
colnames(design) <- c("Intercept", "Diff")

Perform simulations for detecting significant changes in protein abundance, Figure 3 (left):

set.seed(1)

fdr <- seq(0.01,0.1,0.01)
nf <- length(fdr)

mtp.mod <- mfp.mod <- mtp.ord <- mfp.ord <- matrix(nrow=nf,ncol=nit)

for(i in 1:nit){
  v <- rchisq(n,d0)
  v <- d0*s02/v
  x <- matrix(NA,nrow=n,ncol=8)
  for(k in 1:n) x[k,] <- round(rnorm(8,0,sqrt(v[k])),3)
  x[1:hm,1:4] <- x[1:hm,1:4] - log2(1.5)
  fit <- lmFit(x,design)
  fit.eb <- eBayes(fit)
  t.ord <- fit.eb$coefficients[,2]/fit.eb$sigma/fit.eb$stdev.unscaled[,2]
  t.mod <- fit.eb$t[, 2]
  p.ord <- 2*pt(-abs(t.ord),fit.eb$df.residual)
  p.mod <- fit.eb$p.value[,2]
  z <- data.frame(p.ord,p.mod)
  z$q.ord <- qvalue(p.ord)$q
  z$q.mod <- qvalue(p.mod)$q
  for(k in 1:nf){
    mtp.ord[k,i] <- sum(z$q.ord[1:hm]<fdr[k])
    mtp.mod[k,i] <- sum(z$q.mod[1:hm]<fdr[k])
    mfp.ord[k,i] <- sum(z$q.ord[-(1:hm)]<fdr[k])
    mfp.mod[k,i] <- sum(z$q.mod[-(1:hm)]<fdr[k])
  }
}

vtp.ord <- apply(mtp.ord,1,mean)
vtp.mod <- apply(mtp.mod,1,mean)
vfp.ord <- apply(mfp.ord,1,mean)
vfp.mod <- apply(mfp.mod,1,mean)

mcol <- c(brewer.pal(9,"Blues")[c(5,5,8)],brewer.pal(9,"Reds")[c(4,4,7)])
par(las=1,mar=c(7,5,5,0),font.lab=2,cex.lab=1.7,font.axis=2,cex.axis=1.7,xaxs="i",yaxs="i")
plot(c(0.005,0.105),c(0,100),type="n",
     xlab="false discovery rate",
     ylab="declared differentially expressed", xaxt="n")
abline(v=c(0,fdr),h=seq(0,100,10),col="lightgray",lty="dotted")
lines(fdr,vtp.ord,col=mcol[1],type="b",cex=0.6,lwd=0.5,pch=15)
lines(fdr,vtp.mod,col=mcol[4],type="b",cex=0.6,lwd=0.5,pch=15)
lines(fdr,vfp.ord,col=mcol[2],type="b",cex=0.9,lwd=0.5,pch=18)
lines(fdr,vfp.mod,col=mcol[5],type="b",cex=0.9,lwd=0.5,pch=18)
lines(fdr,vfp.ord+vtp.ord,col=mcol[3],type="b",cex=1.1,lwd=2,pch=19)
lines(fdr,vfp.mod+vtp.mod,col=mcol[6],type="b",cex=1.1,lwd=2,pch=19)
axis(1,seq(0.01,0.1,0.01),paste0(seq(1,10,1),"%"))
axis(2,seq(0,100,10),rep("",length(seq(0,100,10))),tck=-0.01)
legend("topleft",col=mcol,lty=1,pch=rep(c(15,18,19),2),bg="white",cex=1.5,
       paste0(rep(c("ordinary","moderated"),rep(3,2))," ",
              rep(c("true positives","false positives","total"),2)))

Perform simulations for ROC curves, Figure 3 (right):

set.seed(1)

p.mod <- p.ord <- matrix(nrow=n,ncol=nit)

for(i in 1:nit){
  v <- rchisq(n,d0)
  v <- d0*s02/v
  x <- matrix(NA,nrow=n,ncol=8)
  for(k in 1:n) x[k,] <- round(rnorm(8,0,sqrt(v[k])),3)
  x[1:hm,1:4] <- x[1:hm,1:4] - log2(1.5)
  fit <- lmFit(x,design)
  fit.eb <- eBayes(fit)
  t.ord <- fit.eb$coefficients[,2]/fit.eb$sigma/fit.eb$stdev.unscaled[,2]
  p.ord[,i] <- 2*pt(-abs(t.ord),fit.eb$df.residual)
  p.mod[,i] <- fit.eb$p.value[,2]
}

nfp <- 21
mfp <- 0.1
x.ord <- x.mod <- matrix(nrow=nfp,ncol=nit)
fp <- seq(0,mfp,length=nfp)

for(i in 1:nit){
  p <- p.ord[-(1:hm),i]
  cut <- quantile(p,fp)
  p <- p.ord[1:hm,i]
  x.ord[,i] <- unlist(lapply(lapply(cut,">",p),mean))
  p <- p.mod[-(1:hm),i]
  cut <- quantile(p,fp)
  p <- p.mod[1:hm,i]
  x.mod[,i] <- unlist(lapply(lapply(cut,">",p),mean))
}

tp.ord <- apply(x.ord,1,mean)
tp.mod <- apply(x.mod,1,mean)

mcol <- c(brewer.pal(9,"Blues")[8],brewer.pal(9,"Reds")[7],brewer.pal(9,"Greens")[7])
par(las=1,mar=c(5,7,5,7),xaxs="i",yaxs="i",font.lab=2,cex.lab=1.7,font.axis=2,cex.axis=1.7)
plot(c(-0.005,mfp+0.005),c(0.5,1),xlab="false positive rate",ylab="",type="n",xaxt="n",yaxt="n")
abline(v=seq(0,0.1,0.01),h=seq(0.5,1,0.05),col="lightgray",lty="dotted")
lines(fp[-1],tp.ord[-1],col=mcol[1],type="b",cex=1.1,lwd=2,pch=19)
lines(fp[-1],tp.mod[-1],col=mcol[2],type="b",cex=1.1,lwd=2,pch=19)
lines(fp[-1],2.0*(tp.mod-tp.ord)[-1]+0.5,col=mcol[3],type="b",cex=0.7,lwd=0.5,pch=19)
axis(1,seq(0,0.1,0.005),rep("",length(seq(0,0.1,0.005))),tck=-0.01)
axis(1,seq(0,0.1,0.01),rep("",length(seq(0,0.1,0.01))),tck=-0.02)
axis(1,seq(0,0.1,0.05),paste0(seq(0,10,5),"%"))
axis(2,seq(0.5,1,0.1),paste0(seq(50,100,10),"%"))
axis(4,seq(0.5,1,length=6),paste0(seq(0,25,5),"%"),col.axis=mcol[3],col.ticks=mcol[3])
mtext("improvement in sensitivity",side=4,line=5,col=mcol[3],las=0,font=2,cex=1.7)
mtext("true positive rate",side=2,line=5,las=0,font=2,cex=1.7)
legend("bottomright",col=mcol[1:2],lty=1,pch=rep(19,2),bg="white",cex=1.5,c("ordinary","moderated"))