R code: Figure 3

Load the following packages:

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

Set parameters:

d0 <- 4 # prior degrees of freedom
s02 <- 0.035 # prior variance
n <- 5000 # number of proteins

hm <- 100 # number of significant proteins
hm.fc <- 1.50 # fold change for significant proteins

Coding of treatment (tr) and control (ct) groups according to limma-package:

design <- model.matrix(~factor(c(rep(1,5), rep(2,5))))
colnames(design) <- c("Intercept", "tr-ct")

Simulate data and calculate ordinary and moderated statistics:

set.seed(1)

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

v <- rchisq(n,d0)
v <- d0*s02/v

x <- matrix(NA, nrow=n, ncol=10)
for(k in 1:n) x[k,] <- round(rnorm(10,0,sqrt(v[k])),5)
x[1:hm,1:5] <- x[1:hm,1:5] - log2(hm.fc)
fit <- lmFit(x, design)
fit.eb <- eBayes(fit)
log2FC <- fit.eb$coefficients[, 2]
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]
q.ord <- qvalue(p.ord)$q
q.mod <- qvalue(p.mod)$q

Create volcano plots:

# save significant/non-significant IDs
sig.ord <- which(q.ord<0.05)
sig.mod <- which(q.mod<0.05)
nsig.ord <- which(q.ord>0.05)
nsig.mod <- which(q.mod>0.05)

# split according to index < / > 100 into tp, fp, tn, and fn:
sig.ord.tp <- sig.ord[sig.ord <=100]
sig.ord.fp <- sig.ord[sig.ord >=100]
sig.mod.tp <- sig.mod[sig.mod <=100]
sig.mod.fp <- sig.mod[sig.mod >=100]
sig.ord.tn <- nsig.ord[nsig.ord >=100]
sig.ord.fn <- nsig.ord[nsig.ord <=100]
sig.mod.tn <- nsig.mod[nsig.mod >=100]
sig.mod.fn <- nsig.mod[nsig.mod <=100]

# volcano plot
par(mfrow=c(1,2), font.lab=2, cex.lab=1.2, font.axis=2, cex.axis=1.2)
par(las=1, xaxs="i", yaxs="i")
rx <- c(-1, 1)*1.6
ry <- c(0, 6.5)

plot(log2FC, -log10(p.ord), pch=21, bg="lightgrey", cex=0.5, 
     xlim=rx, ylim=ry, xaxt="n",
     xlab="fold change", ylab="-log10  p-value")
abline(v=seq(-2,2,1), col="lightgray", lty="dotted")
abline(h=seq(0,ry[2],1), col="lightgray", lty="dotted")
axis(1, seq(-2,2,1), paste(c("1/4","1/2","1/1","2/1","4/1")))
title("volcano plot of ordinary p-values")

points(log2FC, -log10(p.ord),
       pch=21, bg="lightgrey", cex=0.5)

# mark true positives
points(log2FC[sig.ord.tp], -log10(p.ord[sig.ord.tp]),
       pch=21, bg="blue", cex=1.3)

# mark false positives
points(log2FC[sig.ord.fp], -log10(p.ord[sig.ord.fp]),
       pch=21, bg="red", cex=1.3)

# mark false negatives
points(log2FC[sig.ord.fn], -log10(p.ord[sig.ord.fn]),
       pch=21, bg="yellow", cex=1.3)

plot(log2FC, -log10(p.mod), pch=21, bg="lightgrey", cex=0.5,
     xlim=rx, ylim=ry, xaxt="n",
     xlab="fold change", ylab="-log10  p-value")
abline(v=seq(-2,2,1), col="lightgray", lty="dotted")
abline(h=seq(0,6,1), col="lightgray", lty="dotted")
axis(1, seq(-2,2,1), paste(c("1/4","1/2","1/1","2/1","4/1")))
title("volcano plot of moderated p-values")

points(log2FC, -log10(p.mod),
       pch=21, bg="lightgrey", cex=0.5)

# mark true positives
points(log2FC[sig.mod.tp], -log10(p.mod[sig.mod.tp]),
       pch=21, bg="blue", cex=1.3)

# mark false positives
points(log2FC[sig.mod.fp], -log10(p.mod[sig.mod.fp]),
       pch=21, bg="red", cex=1.3)

# mark false negatives
points(log2FC[sig.mod.fn], -log10(p.mod[sig.mod.fn]),
       pch=21, bg="yellow", cex=1.3)