#### Test 1

  test.1.sam <- read.csv("test-1-sam.csv",header=T)
      data.1 <- test.1.sam[,3:43]; n.a <- 20; n.b<-21; genes <- test.1.sam[,2]
      data.es <- data.1
      data.es <- data.1[, c(1:10,21:30)]; n.a <- 10; n.b<-10; genes <- test.1.sam[,2]
      data.es <- data.1[, c(1:10,32:41)]; n.a <- 10; n.b<-10; genes <- test.1.sam[,2]

    #### make 2 tests for 10 and 20 samples
     class.label <- c(rep(0,10), rep(1,10))    
     class.label <- c(rep(0,20), rep(1,21))   

    ##### SAM
       sam.output<-sam(data.es,class.label,med.fdr=FALSE,rand=123)
         sam.res <- sam.plot(sam.output,delta=1.8,file.out="sam-res.txt",gene.names=as.character(genes))
      intersect(sam.res$sam.output[,6],genes[1:12])
    
  #### BADGE:
           es.ln    <- badge.lnorm.dv(data.es,genes,n.a,n.b)
           es.gamma <- badge.gamma(data.es,genes,n.a,n.b,1)
           pos.wei <- posterior.weights(data.es)
           b.fac <- exp(pos.wei[1,]-pos.wei[2,])
            wei.gamma <-  b.fac/(1+b.fac)
            wei.lnorm <-  1-wei.gamma

  theta.es  <- wei.gamma* es.gamma[, 2]+wei.lnorm* es.ln[,2]
  prob.es   <- (wei.gamma* es.gamma[,3]+wei.lnorm* es.ln[, 3]) # prob (theta > 1)
    alpha <- c(.025,.01,.005,.001)
     for(i in 1:4){
       fpr  <- union(which(prob.es>(1-alpha[i])), which(prob.es< alpha[i]))
          print(length(fpr))
      print(length(intersect(genes[fpr], genes[1:12])))}

     ### use false discovery rate
       fdr <- c(); tdr <- c()
          ind.sort <- sort.list(prob.es); prob.sort <- prob.es[ind.sort]
         prob.sort.b <- 1-prob.sort
            alpha <- c(.2,.1,.05,.025);
           n.test<- length(theta.es)
           for(ind.a in 1:length(alpha)){
            selection <- union( which( (prob.sort.b) < c(1:n.test)*alpha[ind.a]/(n.test)) , 
                                 which(    prob.sort < c(1:n.test)*alpha[ind.a]/(n.test)) )
                fdr <- c(fdr, length(selection))
      selection.top <-  which( (prob.sort.b) < c(1:n.test)*alpha[ind.a]/(n.test)) 
      selection.bot <-  which(    prob.sort < c(1:n.test)*alpha[ind.a]/(n.test)) 
                  print(length(intersect(genes[ind.sort][union(selection.top,selection.bot)], genes[1:12])))}


  ##### EBAM
   find.out<-find.a0(data.es,class.label,delta=.9,alpha=c(0,0.05,0.1),rand=123)
                    ebam.out<-ebam(find.out,gene.names=genes,R.fold=T)

  #### T-test on log-transformed data
       t.ln <- ttest.log(data.es,genes,n.a,n.b)
       prob.tln <-  (pt((t.ln), (n.a+n.b)))
               ind.sort <- sort.list(prob.tln); prob.sort <- prob.tln[ind.sort]
    for(i in 1:4){
       fpr  <- union(which(prob.sort>(1-alpha[i])), which(prob.sort< alpha[i]))
          print(length(fpr))
      print(length(intersect(genes[fpr], genes[1:12])))}
    

     ### use false discovery rate
       fdr <- c(); tdr <- c()
          ind.sort <- sort.list(prob.tln); prob.sort <- prob.tln[ind.sort]
         prob.sort.b <- 1-prob.sort
            alpha <- c(.2,.1,.05,.025);
           n.test<- length(theta.es)
           for(ind.a in 1:length(alpha)){
            selection <- union( which( (prob.sort.b) < c(1:n.test)*alpha[ind.a]/(n.test)) , 
                                 which(    prob.sort < c(1:n.test)*alpha[ind.a]/(n.test)) )
                fdr <- c(fdr, length(selection))
      selection.top <-  which( (prob.sort.b) < c(1:n.test)*alpha[ind.a]/(n.test)) 
      selection.bot <-  which(    prob.sort < c(1:n.test)*alpha[ind.a]/(n.test)) 
                  print(length(intersect(genes[ind.sort][union(selection.top,selection.bot)], genes[1:12])))}
      fdr

  #### T-test on original data

      tval <- ttest(data.es,genes,n.a,n.b)
                     prob.t <- pt((tval), (n.a+n.b))
    for(i in 1:4){
       fpr  <- union(which(prob.t>(1-alpha[i])), which(prob.t< alpha[i]))
          print(length(fpr))
      print(length(intersect(genes[fpr], genes[1:12])))}

    ### use false discovery rate
       fdr <- c(); tdr <- c()
          ind.sort <- sort.list(prob.tln); prob.sort <- prob.t[ind.sort]
         prob.sort.b <- 1-prob.sort
            alpha <- c(.2,.1,.05,.025);
           n.test<- length(theta.es)
           for(ind.a in 1:length(alpha)){
            selection <- union( which( (prob.sort.b) < c(1:n.test)*alpha[ind.a]/(n.test)) , 
                                 which(    prob.sort < c(1:n.test)*alpha[ind.a]/(n.test)) )
                fdr <- c(fdr, length(selection))
      selection.top <-  which( (prob.sort.b) < c(1:n.test)*alpha[ind.a]/(n.test)) 
      selection.bot <-  which(    prob.sort < c(1:n.test)*alpha[ind.a]/(n.test)) 
                  print(length(intersect(genes[ind.sort][union(selection.top,selection.bot)], genes[1:12])))}
      fdr