Golubsim.R

Code snipits and programs written in R, S or S-PLUS{{#security:*|Sven}} nCr <- function(n,r) {

 factorial(n) / (factorial(n-r) * factorial(r))

} golubsim <- function(outfile = "tmp", iiter=nCr(11,6), alpha = 0.05, seed=1,

                    datadir = "Data", lambdaseq = seq(0, 0.95, by=0.05),
                    ncs.value="paper", use.siggenes=TRUE, use.convest = TRUE
                    ) {

1. ----------------------------------------------------------------------------- #
2. itter = numeric number of iterations
3. alpha = numeric number [0,1] type I probability test cutoff
4. use.siggenes = logical to use pi0.est function in siggenes package
5. lambdaseq = numeric vector lambda range for pi0.est smoothing splines
6. ncs.value = character string "paper " or "max"
7. use.convest = logical to use convest function in limma package
8. ----------------------------------------------------------------------------- #

1. 1) ---------------------------------- Setup --------------------------------- #
2. Check and create directories for save(...) of objects

 if(length(dir(datadir))==0) {
   create.dir(datadir)
 }

1. Check and install multtest

 if(!require(multtest, quietly=TRUE)) {
   cat("Package:Multtest not available attempting to install\n")
   install.packages("multtest")
 }

 data(golub)
 m <- nrow(golub)

1. AML slides=11

 n <- 11
 r <-  6

 AML.samples <- list()
 i <- 1
 while(i <= iiter) {
   value <- sort(sample(n, r, replace=FALSE))
   key   <- paste(value, collapse=" ")
   if(is.null(AML.sampleskey)) {
     AML.sampleskey <- value
     i <-i+1

1. if(!(i %% 10)){ # print every 10th
2. print(i)
3. }

   }
 }

1. ALL slides=27

 n <- 27
 r <-  6
 
 ALL.samples <- list()
 i <- 1
 while(i <= iiter) {
   value <- sort(sample(n, r, replace=FALSE))
   key   <- paste(value, collapse=" ")
   if(is.null(ALL.sampleskey)) {
     ALL.sampleskey <- value
     i <-i+1

1. if(!(i %% 10)){ # print every 10th
2. print(i)
3. }

   }
 }

cat("Column permutations sampled\n")

1. Need to randomise the way the keys are drawn

 ALLdraw <- sample(iiter,iiter, replace=FALSE)
 AMLdraw <- sample(iiter,iiter, replace=FALSE)

 golub.ALL <- golub[,golub.cl==0]
 golub.AML <- golub[,golub.cl==1]
 
 golub.clSample <- rep(0:1, each=6)
 
 design.Sample <- cbind("ALL"=1, "AMLvsALL"=golub.clSample)
 

1. Check and install limma

 if(!require(limma, quietly=TRUE)) {
   cat("Package:limma not available attempting to install\n")
   install.packages("limma")
 }

1. Check and install siggenes

 if(!require(siggenes, quietly=TRUE)) {
   cat("Package:siggenes not available attempting to install\n")
   install.packages("siggenes")
 }

 set.seed(seed)
 starttime <- proc.time()[3]

1. Output lists

 FDlist <- list(m=m, alpha=alpha, p0=rep(NA,iiter), p0NW=rep(NA,iiter), convest=rep(NA,iiter))

 alphastar <- alpha # Do we need this?

1. 3) --------------------------- iiter loop ----------------------------- #

 for(i in seq(iiter)){

1. 4) ------------------------ Permutation analysis ---------------------------- #

   x <- cbind(golub.ALL[,ALL.samples[[ALLdraw[i]]]],
              golub.AML[,AML.samples[[AMLdraw[i]]]]
              )

   fit <- lmFit(x, design.Sample) 
   fit <- eBayes(fit)
   pvalues <- fit$p.value[,"AMLvsALL"]

1. p.adjusted <- p.adjust(fit$p.value[,"AMLvsALL"], method="BH")
2. pGenes <- sum(p.adjusted < 0.05, na.rm=TRUE)

1. 5) ----------------------------- Saving statistics -------------------------- #
2. pi0 estimation
3. Storey uses seq(0,0.95, by=0.05) for lambda, estimates at pi(0.95)
4. Maximum value of lambda needs to be less than 1 estimated by (m-1)/m

   if(use.convest) {# Convest estimation
     FDlist"convest"[i] <- convest(pvalues, niter=50)
   }
 # Siggenes or first principles
   if( use.siggenes ) {
     require(siggenes, quietly=TRUE)
         # A) weighting by 1-lambda
     pi0hat <- pi0.est(pvalues, lambdaseq, ncs.value=ncs.value,
                       ncs.weights=1-lambdaseq)$p0
     FDlist"p0"[i] <- round(pi0hat,8) 
         # B) no 1-lambda weighting
     pi0hat <- pi0.est(pvalues, lambdaseq, ncs.value=ncs.value,
                       ncs.weights=NULL)$p0
     FDlist"p0NW"[i] <- round(pi0hat,8)
   } else {
     pi0lambda <- tapply(lambdaseq, lambdaseq,
                         function(x){sum(pvalues>x)/(m*(1-x))})
         # A) weighting by 1-lambda
     pi0hat <- min(smooth.spline(lambdaseq, pi0lambda, w=1-lambdaseq, 
                                 df=3)$y[length(lambdaseq)],1)
     FDlist"p0"[i] <- round(pi0hat,8) 
         # B) no 1-lambda weighting
     pi0hat <- min(smooth.spline(lambdaseq, pi0lambda,
                                 df=3)$y[length(lambdaseq)],1)
     FDlist"p0NW"[i] <- round(pi0hat,8) 
   }
   save(FDlist, file=file.path(datadir, outfile), compress=T)
 } # End of main loop

1. Return simulation time

 endtime <- proc.time()[3] - starttime
 return(endtime)

}

1. load(file.path("Data", "tmp"))

1. source("plotfun.R")
2. plotfun(FDlist)