# Exercise 1 - False discovery rate
# Author Roland Krause/Utz Pape
ngenes <- 1000
ndiff <- 50

fgen.data <- function(nrepl)
{
  m.norm <- matrix(rnorm((ngenes-ndiff)*nrepl*2,mean=10,sd=5),ncol=2*nrepl,nrow=ngenes-ndiff)
  m.diff <- matrix(rnorm(ndiff*nrepl,mean=10,sd=5),ncol=nrepl,nrow=ndiff)
  m.add <- matrix(rnorm(ndiff*nrepl,mean=5,sd=1),ncol=nrepl,nrow=ndiff)
  m.diff <- cbind(m.diff,m.diff+m.add)
  m <- rbind(m.norm,m.diff)

  vt <- vector()
  for (i in 1:nrow(m))
  {
    vt[i] <- t.test(m[i,1:nrepl],m[i,(nrepl+1):ncol(m)],paired=T)$p.value
  }
  return(vt)
}

# two replicates
d2 <- sort(fgen.data(2))
d4 <- sort(fgen.data(3))
d8 <- sort(fgen.data(4))

par(mfcol=c(1,2))
lw = 2
plot(sort(d2), type="l", xlab="Genes", ylab="p-value", lwd=lw)
points(sort(d4), col="red", type = "l", lwd=lw)
points(sort(d8), col="green", type = "l", lwd=lw)
# fdr p(kP <= k/m * alpha (m =1000)  - 0.05 
abline(0,1/1000, col= "blue", lwd =lw)
abline(0,0.05/200, col= "orange", lwd = lw)
abline(0,0.05/1000, col= "purple", lwd = lw)

legend(1,1, c(2,3,4, "FDR", "Bonferroni"), col=c("black","red", "green", "orange","purple"), lwd = lw)

xlimit=100
plot(sort(d2), type="l", xlab="Genes", ylab="p-value", lwd=lw, xlim=c(0,xlimit), ylim=c(0,xlimit/ngenes))
points(sort(d4), col="red", type = "l", lwd=lw)
points(sort(d8), col="green", type = "l", lwd=lw)
# fdr p(kP <= k/m * alpha (m =1000)  - 0.05 
abline(0,1/1000, col= "blue", lwd =lw)
abline(0,0.05/200, col= "orange", lwd = lw)
abline(0,0.05/1000, col= "purple", lwd = lw)