FN.fdr <- function(fn, fp, target) {
  out <- c()
  for (j in 1:length(target)) {
    temp <- c()
    for (k in 1:100) {
      LL <- match(target[j], fp[k,])
      if (is.na(LL)) {
        LL <- max(which(fp[k,]<target[j]), 1)
        UL <- min(which(fp[k,]>target[j]), 100)
      }  else  UL <- max(which(target[j]==fp[k,]))
      temp <- c(temp, mean(fn[k, LL:UL]))
    }
    out <- c(out, mean(temp))
  }
}


FN <- function(fn, fp, target) {
  out <- c()
  for (j in 1:length(target)) {
    temp <- c()
    for (k in 1:100) {
      LL <- match(target[j], fp[k,])
      if (is.na(LL)) {
        LL <- max(which(fp[k,]>target[j]), 1)
        UL <- min(which(fp[k,]<target[j]), 100)
      }  else  UL <- max(which(target[j]==fp[k,]))
      temp <- c(temp, mean(fn[k, LL:UL]))
    }
    out <- c(out, mean(temp))
  }
}

pp <- c(0.01, 0.02, 0.05, 0.1, 0.2)  # true proportion of genes differentially expressed (5, 10, 25, 50, 100)
par(mfrow=c(2,2))

plot(-1,-1, main=paste("p =", pp[4]), xlab="Desired FDR", ylab="Actual FDR", xlim=c(0.05,0.6), ylim=c(0.05,0.6))
points(((1:100)*0.01)[(0:19)*5+1], (apply(gg.fdr$fp[4,,], 2, mean) / (apply(gg.fdr$fp[4,,], 2, mean) + apply(gg.fdr$tp[4,,], 2, mean)))[(0:19)*5+1], type="b", lwd=0.3, pch=4, cex=0.6)
points(((1:100)*0.01)[(0:19)*5+2], (apply(lnn.fdr$fp[4,,], 2, mean) / (apply(lnn.fdr$fp[4,,], 2, mean) + apply(lnn.fdr$tp[4,,], 2, mean)))[(0:19)*5+2], type="b", lwd=0.3, pch=1, cex=0.6)
points(((1:100)*0.01)[(0:19)*5+3], (apply(gg.pr.fdr$fp[4,,], 2, mean) / (apply(gg.pr.fdr$fp[4,,], 2, mean) + apply(gg.pr.fdr$tp[4,,], 2, mean)))[(0:19)*5+3], type="b", lwd=0.3, pch=17, cex=0.6)
points(((1:100)*0.01)[(0:19)*5+4], (apply(elnn.fdr$fp[4,,], 2, mean) / (apply(elnn.fdr$fp[4,,], 2, mean) + apply(elnn.fdr$tp[4,,], 2, mean)))[(0:19)*5+4], type="b", lwd=0.3, pch=16, cex=0.6)
points(((1:100)*0.01)[(0:20)*5+0], (apply(limma.fdr$fp[4,,], 2, mean) / (apply(limma.fdr$fp[4,,], 2, mean) + apply(limma.fdr$tp[4,,], 2, mean)))[(0:20)*5+0], type="b", lwd=0.3, pch=6, cex=0.6)
points(((1:100)*0.01)[(0:19)*5+1], (apply(samr.fdr$fp[4,,], 2, mean) / (apply(samr.fdr$fp[4,,], 2, mean) + apply(samr.fdr$tp[4,,], 2, mean)))[(0:19)*5+1], type="l", lwd=0.3)
points(((1:100)*0.01)[(0:19)*5+2], (apply(bridge.fdr$fp[4,,], 2, mean) / (apply(bridge.fdr$fp[4,,], 2, mean) + apply(bridge.fdr$tp[4,,], 2, mean)))[(0:19)*5+2], type="b", lwd=0.3, pch=15, cex=0.6)
abline(0,1,col="gray")

plot(-10,-10, main=paste("p =", pp[4]), xlab="FN", ylab="FP", xlim=c(10,40), ylim=c(0,80))
target <- (0:16)*5
points(FN(gg.map$fn[4,,], gg.map$fp[4,,], target), target, type="b", lwd=0.3, pch=4, cex=0.6)
points(FN(lnn.map$fn[4,,], lnn.map$fp[4,,], target), target, type="b", lwd=0.3, pch=1, cex=0.6)
points(FN(gg.pr.map$fn[4,,], gg.pr.map$fp[4,,], target), target, type="b", lwd=0.3, pch=17, cex=0.6)
points(FN(elnn.map$fn[4,,], elnn.map$fp[4,,], target), target, type="b", lwd=0.3, pch=16, cex=0.6)
points(FN.fdr(limma.fdr$fn[4,,], limma.fdr$fp[4,,], target), target, type="b", lwd=0.3, pch=6, cex=0.6)
points(FN.fdr(samr.fdr$fn[4,,], samr.fdr$fp[4,,], target), target, type="l", lwd=0.3)
points(FN(bridge.map$fn[4,,], bridge.map$fp[4,,], target), target, type="b", lwd=0.3, pch=15, cex=0.6)


plot(-1,-1, main=paste("p =", pp[5]), xlab="Desired FDR", ylab="Actual FDR", xlim=c(0.05,0.6), ylim=c(0.05,0.6))
points(((1:100)*0.01)[(0:19)*5+1], (apply(gg.fdr$fp[5,,], 2, mean) / (apply(gg.fdr$fp[5,,], 2, mean) + apply(gg.fdr$tp[5,,], 2, mean)))[(0:19)*5+1], type="b", lwd=0.3, pch=4, cex=0.6)
points(((1:100)*0.01)[(0:19)*5+2], (apply(lnn.fdr$fp[5,,], 2, mean) / (apply(lnn.fdr$fp[5,,], 2, mean) + apply(lnn.fdr$tp[5,,], 2, mean)))[(0:19)*5+2], type="b", lwd=0.3, pch=1, cex=0.6)
points(((1:100)*0.01)[(0:19)*5+3], (apply(gg.pr.fdr$fp[5,,], 2, mean) / (apply(gg.pr.fdr$fp[5,,], 2, mean) + apply(gg.pr.fdr$tp[5,,], 2, mean)))[(0:19)*5+3], type="b", lwd=0.3, pch=17, cex=0.6)
points(((1:100)*0.01)[(0:19)*5+4], (apply(elnn.fdr$fp[5,,], 2, mean) / (apply(elnn.fdr$fp[5,,], 2, mean) + apply(elnn.fdr$tp[5,,], 2, mean)))[(0:19)*5+4], type="b", lwd=0.3, pch=16, cex=0.6)
points(((1:100)*0.01)[(0:20)*5+0], (apply(limma.fdr$fp[5,,], 2, mean) / (apply(limma.fdr$fp[5,,], 2, mean) + apply(limma.fdr$tp[5,,], 2, mean)))[(0:20)*5+0], type="b", lwd=0.3, pch=6, cex=0.6)
points(((1:100)*0.01)[(0:19)*5+1], (apply(samr.fdr$fp[5,,], 2, mean) / (apply(samr.fdr$fp[5,,], 2, mean) + apply(samr.fdr$tp[5,,], 2, mean)))[(0:19)*5+1], type="l", lwd=0.3)
points(((1:100)*0.01)[(0:19)*5+2], (apply(bridge.fdr$fp[5,,], 2, mean) / (apply(bridge.fdr$fp[5,,], 2, mean) + apply(bridge.fdr$tp[5,,], 2, mean)))[(0:19)*5+2], type="b", lwd=0.3, pch=15, cex=0.6)
abline(0,1,col="gray")

plot(-10,-10, main=paste("p =", pp[5]), xlab="FN", ylab="FP", xlim=c(15,80), ylim=c(0,80))
target <- (0:16)*5
points(FN(gg.map$fn[5,,], gg.map$fp[5,,], target), target, type="b", lwd=0.3, pch=4, cex=0.6)
points(FN(lnn.map$fn[5,,], lnn.map$fp[5,,], target), target, type="b", lwd=0.3, pch=1, cex=0.6)
points(FN(gg.pr.map$fn[5,,], gg.pr.map$fp[5,,], target), target, type="b", lwd=0.3, pch=17, cex=0.6)
points(FN(elnn.map$fn[5,,], elnn.map$fp[5,,], target), target, type="b", lwd=0.3, pch=16, cex=0.6)
points(FN.fdr(limma.fdr$fn[5,,], limma.fdr$fp[5,,], target), target, type="b", lwd=0.3, pch=6, cex=0.6)
points(FN.fdr(samr.fdr$fn[5,,], samr.fdr$fp[5,,], target), target, type="l", lwd=0.3)
points(FN(bridge.map$fn[5,,], bridge.map$fp[5,,], target), target, type="b", lwd=0.3, pch=15, cex=0.6)

dev.print(device=postscript, "res_lnn_variable_variance3.eps", onefile=FALSE, horizontal=FALSE)