library(EBarrays)
library(limma)
library(samr)
library(bridge)
source("summary.pp.R")
source("gg_pr.R")
source("elnn.R")


n <- 500  # number of genes
pp <- c(0.01, 0.02, 0.05, 0.1, 0.2)  # true proportion of genes differentially expressed (5, 10, 25, 50, 100)
data <- matrix(0,n,6)


temp <- array(0, c(length(pp), 100, 100))     # (proportion index, trial index, cutpoint index) 
tlist <- list(tp=temp, fp=temp, tn=temp, fn=temp, sens=temp, spec=temp, fdr=temp)
gg.map <- gg.fdr <- lnn.map <- lnn.fdr <- gg.pr.map <- gg.pr.fdr <- elnn.map <- elnn.fdr <- limma.fdr <- samr.fdr <- bridge.map <- bridge.fdr <- tlist


for (m in 1:length(pp)) {

  for (k in 1:100) { cat(m,k)

    hmu <- 2
    hs2 <- 1
    aa0 <- 1
    vv <- 20

    aa <- rlnorm(n, hmu, sqrt(hs2))
    theta <- rgamma(n*pp[m], shape=aa0, scale=1/vv)
    for (i in 1:(n*pp[m])) data[i,1:3] <- rgamma(3, shape=aa[i], scale=1/theta[i])
    theta <- rgamma(n*pp[m], shape=aa0, scale=1/vv)
    for (i in 1:(n*pp[m])) data[i,4:6] <- rgamma(3, shape=aa[i], scale=1/theta[i])

    theta <- rgamma(n-n*pp[m], shape=aa0, scale=1/vv)
    for (i in (n*pp[m]+1):n) data[i,] <- rgamma(6, shape=aa[i], scale=1/theta[i-n*pp[m]])



    ##### EBarrays(GG)

    pattern <- ebPatterns(c("1,1,1,1,1,1", "1,1,1,2,2,2"))
    gg.fit <- emfit(data, family="GG", hypotheses=pattern)
    gg.post <- postprob(gg.fit, data)
    Z <- gg.post[,2]

    for (i in 1:100) {
      res <- summary.pp(Z, method="MAP", cutoff=i*0.01)
      gg.map$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      gg.map$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      gg.map$tn[m,k,i] <- n*(1-pp[m]) - gg.map$fp[m,k,i]
      gg.map$fn[m,k,i] <- n*pp[m] - gg.map$tp[m,k,i]
      gg.map$sens[m,k,i] <- gg.map$tp[m,k,i] / (gg.map$tp[m,k,i] + gg.map$fn[m,k,i])
      gg.map$spec[m,k,i] <- gg.map$tn[m,k,i] / (gg.map$tn[m,k,i] + gg.map$fp[m,k,i])
      gg.map$fdr[m,k,i] <- gg.map$fp[m,k,i] / (gg.map$tp[m,k,i] + gg.map$fp[m,k,i])
    }

    for (i in 1:100) {
      res <- summary.pp(Z, method="FDR", cutoff=i*0.01)
      gg.fdr$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      gg.fdr$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      gg.fdr$tn[m,k,i] <- n*(1-pp[m]) - gg.fdr$fp[m,k,i]
      gg.fdr$fn[m,k,i] <- n*pp[m] - gg.fdr$tp[m,k,i]
      gg.fdr$sens[m,k,i] <- gg.fdr$tp[m,k,i] / (gg.fdr$tp[m,k,i] + gg.fdr$fn[m,k,i])
      gg.fdr$spec[m,k,i] <- gg.fdr$tn[m,k,i] / (gg.fdr$tn[m,k,i] + gg.fdr$fp[m,k,i])
      gg.fdr$fdr[m,k,i] <- i*0.01
    }





    ##### EBarrays(LNN)

    pattern <- ebPatterns(c("1,1,1,1,1,1", "1,1,1,2,2,2"))
    lnn.fit <- emfit(data, family="LNN", hypotheses=pattern)
    lnn.post <- postprob(lnn.fit, data)
    Z <- lnn.post[,2]

    for (i in 1:100) {
      res <- summary.pp(Z, method="MAP", cutoff=i*0.01)
      lnn.map$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      lnn.map$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      lnn.map$tn[m,k,i] <- n*(1-pp[m]) - lnn.map$fp[m,k,i]
      lnn.map$fn[m,k,i] <- n*pp[m] - lnn.map$tp[m,k,i]
      lnn.map$sens[m,k,i] <- lnn.map$tp[m,k,i] / (lnn.map$tp[m,k,i] + lnn.map$fn[m,k,i])
      lnn.map$spec[m,k,i] <- lnn.map$tn[m,k,i] / (lnn.map$tn[m,k,i] + lnn.map$fp[m,k,i])
      lnn.map$fdr[m,k,i] <- lnn.map$fp[m,k,i] / (lnn.map$tp[m,k,i] + lnn.map$fp[m,k,i])
    }

    for (i in 1:100) {
      res <- summary.pp(Z, method="FDR", cutoff=i*0.01)
      lnn.fdr$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      lnn.fdr$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      lnn.fdr$tn[m,k,i] <- n*(1-pp[m]) - lnn.fdr$fp[m,k,i]
      lnn.fdr$fn[m,k,i] <- n*pp[m] - lnn.fdr$tp[m,k,i]
      lnn.fdr$sens[m,k,i] <- lnn.fdr$tp[m,k,i] / (lnn.fdr$tp[m,k,i] + lnn.fdr$fn[m,k,i])
      lnn.fdr$spec[m,k,i] <- lnn.fdr$tn[m,k,i] / (lnn.fdr$tn[m,k,i] + lnn.fdr$fp[m,k,i])
      lnn.fdr$fdr[m,k,i] <- i*0.01
    }





    ##### GG_prior_ln_MM


    gg.pr.fit <- gg.pr(data, pattern="1 1 1 2 2 2")
    Z <- gg.pr.fit$z

    for (i in 1:100) {
      res <- summary.pp(Z, method="MAP", cutoff=i*0.01)
      gg.pr.map$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      gg.pr.map$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      gg.pr.map$tn[m,k,i] <- n*(1-pp[m]) - gg.pr.map$fp[m,k,i]
      gg.pr.map$fn[m,k,i] <- n*pp[m] - gg.pr.map$tp[m,k,i]
      gg.pr.map$sens[m,k,i] <- gg.pr.map$tp[m,k,i] / (gg.pr.map$tp[m,k,i] + gg.pr.map$fn[m,k,i])
      gg.pr.map$spec[m,k,i] <- gg.pr.map$tn[m,k,i] / (gg.pr.map$tn[m,k,i] + gg.pr.map$fp[m,k,i])
      gg.pr.map$fdr[m,k,i] <- gg.pr.map$fp[m,k,i] / (gg.pr.map$tp[m,k,i] + gg.pr.map$fp[m,k,i])
    }

    for (i in 1:100) {
      res <- summary.pp(Z, method="FDR", cutoff=i*0.01)
      gg.pr.fdr$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      gg.pr.fdr$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      gg.pr.fdr$tn[m,k,i] <- n*(1-pp[m]) - gg.pr.fdr$fp[m,k,i]
      gg.pr.fdr$fn[m,k,i] <- n*pp[m] - gg.pr.fdr$tp[m,k,i]
      gg.pr.fdr$sens[m,k,i] <- gg.pr.fdr$tp[m,k,i] / (gg.pr.fdr$tp[m,k,i] + gg.pr.fdr$fn[m,k,i])
      gg.pr.fdr$spec[m,k,i] <- gg.pr.fdr$tn[m,k,i] / (gg.pr.fdr$tn[m,k,i] + gg.pr.fdr$fp[m,k,i])
      gg.pr.fdr$fdr[m,k,i] <- i*0.01
    }





    ##### eLNN

    elnn.fit <- elnn(data, pattern="1 1 1 2 2 2")
    Z <- elnn.fit$z

    for (i in 1:100) {
      res <- summary.pp(Z, method="MAP", cutoff=i*0.01)
      elnn.map$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      elnn.map$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      elnn.map$tn[m,k,i] <- n*(1-pp[m]) - elnn.map$fp[m,k,i]
      elnn.map$fn[m,k,i] <- n*pp[m] - elnn.map$tp[m,k,i]
      elnn.map$sens[m,k,i] <- elnn.map$tp[m,k,i] / (elnn.map$tp[m,k,i] + elnn.map$fn[m,k,i])
      elnn.map$spec[m,k,i] <- elnn.map$tn[m,k,i] / (elnn.map$tn[m,k,i] + elnn.map$fp[m,k,i])
      elnn.map$fdr[m,k,i] <- elnn.map$fp[m,k,i] / (elnn.map$tp[m,k,i] + elnn.map$fp[m,k,i])
    }

    for (i in 1:100) {
      res <- summary.pp(Z, method="FDR", cutoff=i*0.01)
      elnn.fdr$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      elnn.fdr$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      elnn.fdr$tn[m,k,i] <- n*(1-pp[m]) - elnn.fdr$fp[m,k,i]
      elnn.fdr$fn[m,k,i] <- n*pp[m] - elnn.fdr$tp[m,k,i]
      elnn.fdr$sens[m,k,i] <- elnn.fdr$tp[m,k,i] / (elnn.fdr$tp[m,k,i] + elnn.fdr$fn[m,k,i])
      elnn.fdr$spec[m,k,i] <- elnn.fdr$tn[m,k,i] / (elnn.fdr$tn[m,k,i] + elnn.fdr$fp[m,k,i])
      elnn.fdr$fdr[m,k,i] <- i*0.01
    }





    ##### limma

    design <- model.matrix(~ -1+factor(c(1,1,1,2,2,2)))
    colnames(design) <- c("group1", "group2")  
    limma.fit <- lmFit(log(data), design)
    contrast.matrix <- makeContrasts(group2-group1, levels=design)
    limma.fit2 <- contrasts.fit(limma.fit, contrast.matrix)
    limma.fit2 <- eBayes(limma.fit2)
    p.val <- p.adjust(limma.fit2$p.value, method="fdr")
 
    for (i in 1:100) {
      limma.fdr$tp[m,k,i] <- sum( (p.val < i*0.01)[1:(n*pp[m])] )
      limma.fdr$fp[m,k,i] <- sum( (p.val < i*0.01)[(n*pp[m]+1):n] )
      limma.fdr$tn[m,k,i] <- n*(1-pp[m]) - limma.fdr$fp[m,k,i]
      limma.fdr$fn[m,k,i] <- n*pp[m] - limma.fdr$tp[m,k,i]
      limma.fdr$sens[m,k,i] <- limma.fdr$tp[m,k,i] / (limma.fdr$tp[m,k,i] + limma.fdr$fn[m,k,i])
      limma.fdr$spec[m,k,i] <- limma.fdr$tn[m,k,i] / (limma.fdr$tn[m,k,i] + limma.fdr$fp[m,k,i])
      limma.fdr$fdr[m,k,i] <- limma.fdr$fp[m,k,i] / (limma.fdr$tp[m,k,i] + limma.fdr$fp[m,k,i])
    }





    ##### SAM

    y <- c(1,1,1,2,2,2)
    data.list <- list(x=data, y=y, geneid=as.character(1:nrow(data)), genenames=paste("g",as.character(1:nrow(data)), sep=""), logged2=F)
    res <- samr(data.list, resp.type="Two class unpaired", s0=NULL, s0.perc=NULL, nperms=100, center.arrays=FALSE, testStatistic="standard")
    delta.table <- samr.compute.delta.table(res, nvals=500)

    j <- 1
    for (i in 100:1) {
      while(!is.na(delta.table[j,5]) && j<500 && delta.table[j,5] > i*0.01)    j <- j+1
      delta <- delta.table[j,1]
      siggenes.table <- samr.compute.siggenes.table(res, delta, data.list, delta.table) 
      de <- sort(as.integer( c(siggenes.table$genes.up[,3], siggenes.table$genes.lo[,3]) ))

      samr.fdr$tp[m,k,i] <- sum(de <= n*pp[m])
      samr.fdr$fp[m,k,i] <- length(de) - samr.fdr$tp[m,k,i]
      samr.fdr$tn[m,k,i] <- n*(1-pp[m]) - samr.fdr$fp[m,k,i]
      samr.fdr$fn[m,k,i] <- n*pp[m] - samr.fdr$tp[m,k,i]
      samr.fdr$sens[m,k,i] <- samr.fdr$tp[m,k,i] / (samr.fdr$tp[m,k,i] + samr.fdr$fn[m,k,i])
      samr.fdr$spec[m,k,i] <- samr.fdr$tn[m,k,i] / (samr.fdr$tn[m,k,i] + samr.fdr$fp[m,k,i])
      samr.fdr$fdr[m,k,i] <- i*0.01
    }






    ##### bridge

    for (i in 1:6) data[,i] <- exp(log(data[,i])-mean(log(data[,i])))
    res <- bridge.2samples(data[,1:3],data[,4:6],B=21000,min.iter=1000,batch=10,all.out=FALSE,affy=TRUE, verbose=FALSE,log=FALSE)
    Z <- res$post.p


    for (i in 1:100) {
      res <- summary.pp(Z, method="MAP", cutoff=i*0.01)
      bridge.map$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      bridge.map$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      bridge.map$tn[m,k,i] <- n*(1-pp[m]) - bridge.map$fp[m,k,i]
      bridge.map$fn[m,k,i] <- n*pp[m] - bridge.map$tp[m,k,i]
      bridge.map$sens[m,k,i] <- bridge.map$tp[m,k,i] / (bridge.map$tp[m,k,i] + bridge.map$fn[m,k,i])
      bridge.map$spec[m,k,i] <- bridge.map$tn[m,k,i] / (bridge.map$tn[m,k,i] + bridge.map$fp[m,k,i])
      bridge.map$fdr[m,k,i] <- bridge.map$fp[m,k,i] / (bridge.map$tp[m,k,i] + bridge.map$fp[m,k,i])
    }

    for (i in 1:100) {
      res <- summary.pp(Z, method="FDR", cutoff=i*0.01)
      bridge.fdr$tp[m,k,i] <- sum(res$ind1[1:(n*pp[m])])
      bridge.fdr$fp[m,k,i] <- sum(res$ind1[(n*pp[m]+1):n])
      bridge.fdr$tn[m,k,i] <- n*(1-pp[m]) - bridge.fdr$fp[m,k,i]
      bridge.fdr$fn[m,k,i] <- n*pp[m] - bridge.fdr$tp[m,k,i]
      bridge.fdr$sens[m,k,i] <- bridge.fdr$tp[m,k,i] / (bridge.fdr$tp[m,k,i] + bridge.fdr$fn[m,k,i])
      bridge.fdr$spec[m,k,i] <- bridge.fdr$tn[m,k,i] / (bridge.fdr$tn[m,k,i] + bridge.fdr$fp[m,k,i])
      bridge.fdr$fdr[m,k,i] <- i*0.01
    }




  }

}



save(gg.map, gg.fdr, lnn.map, lnn.fdr, gg.pr.map, gg.pr.fdr, elnn.map, elnn.fdr, limma.fdr, samr.fdr, bridge.map, bridge.fdr, file="dat_gg_variable_variance3.Rdata")