simulatorZ-vignette.R
In simulatorZ: Simulator for Collections of Independent Genomic Data Sets

### R code from vignette source 'simulatorZ-vignette.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: simulatorZ-vignette.Rnw:57-58
###################################################
library(simulatorZ)


###################################################
### code chunk number 3: simulatorZ-vignette.Rnw:93-97 (eval = FALSE)
###################################################
## library(curatedOvarianData)
## source(system.file("extdata",
##                    "patientselection.config",package="curatedOvarianData"))
## source(system.file("extdata", "createEsetList.R", package="curatedOvarianData"))


###################################################
### code chunk number 4: simulatorZ-vignette.Rnw:104-109
###################################################
library(curatedOvarianData)
data(GSE17260_eset)
data(E.MTAB.386_eset)
data(GSE14764_eset)
esets.orig <- list(GSE17260=GSE17260_eset, E.MTAB.386=E.MTAB.386_eset, GSE14764=GSE14764_eset)


###################################################
### code chunk number 5: cleanEsets
###################################################
cleanEsets <- function(obj, keep.common.only=TRUE, meta.required=NULL){
    if(keep.common.only){
        intersect.features <- Reduce(intersect, lapply(obj, featureNames))
        for (i in seq_along(obj))
            obj[[i]] <- obj[[i]][intersect.features, ]
    }
    if(!is.null(meta.required)){
        obj <- lapply(obj, function(obj1){
            for (x in meta.required)
                obj1 <- obj1[, !is.na(obj1[[x]])]
            if(ncol(obj1) > 0){
                return(obj1)
            }else{
                return(NULL)
            }
        })
        return(obj[!sapply(obj, is.null)])
    }
}


###################################################
### code chunk number 6: simulatorZ-vignette.Rnw:143-146
###################################################
esets <- cleanEsets(esets.orig, meta.required=c("days_to_death", "vital_status"))
esets <- lapply(esets, function(eset) eset[1:1000, ])
sapply(esets, dim)


###################################################
### code chunk number 7: simulatorZ-vignette.Rnw:158-161
###################################################
set.seed(8)
sim.esets <- simData(esets, n.samples=500)
names(sim.esets)


###################################################
### code chunk number 8: simulatorZ-vignette.Rnw:183-184
###################################################
sim.esets <- simData(esets, 500, type="one-step")


###################################################
### code chunk number 9: simulatorZ-vignette.Rnw:194-195
###################################################
cleaned.esets <- cleanEsets(esets.orig, meta.required="tumorstage")


###################################################
### code chunk number 10: simulatorZ-vignette.Rnw:199-203
###################################################
sim.sets <- simData(cleaned.esets, 500, 
                    balance.variables="tumorstage")
sim.sets$prob.desired
sim.sets$prob.real[[1]]


###################################################
### code chunk number 11: simulatorZ-vignette.Rnw:209-210
###################################################
cleaned.esets <- cleanEsets(esets.orig, meta.required=c("tumorstage", "grade"))


###################################################
### code chunk number 12: simulatorZ-vignette.Rnw:214-218
###################################################
sim.sets <- simData(cleaned.esets, 500, 
                    balance.variables=c("tumorstage", "grade"))
sim.sets$prob.desired
sim.sets$prob.real[[1]]


###################################################
### code chunk number 13: simulatorZ-vignette.Rnw:228-231
###################################################
y.list <- lapply(esets, function(eset){
    return( Surv(eset$days_to_death, eset$vital_status=="deceased") )
  })


###################################################
### code chunk number 14: simulatorZ-vignette.Rnw:233-234
###################################################
true.mod <- getTrueModel(esets, y.list, 100)


###################################################
### code chunk number 15: simulatorZ-vignette.Rnw:236-237
###################################################
names(true.mod)


###################################################
### code chunk number 16: simulatorZ-vignette.Rnw:240-242
###################################################
simmodel <- simData(esets, 500, y.list)
new.esets <- simTime(simmodel, y.list, true.mod)


###################################################
### code chunk number 17: simulatorZ-vignette.Rnw:249-252
###################################################
true.mod <- getTrueModel(esets, y.list, 100)
sim.sets <- simData(esets, 500, y.list)
sim.sets <- simTime(sim.sets, y.list, true.mod)


###################################################
### code chunk number 18: simulatorZ-vignette.Rnw:256-257
###################################################
sim.sets <- simBootstrap(esets, y.list, 500, 100)


###################################################
### code chunk number 19: simulatorZ-vignette.Rnw:267-281
###################################################
y.vars=y.list
balance.variables=c("tumorstage", "grade")
X.list <- lapply(esets, function(eset){
  newx <- t(exprs(eset))
  return(newx)
}) 
all.X <- do.call(rbind, X.list)
cov.list <- lapply(esets, function(eset){
  return(pData(eset)[, balance.variables])
})  
all.cov <- as(data.frame(do.call(rbind, cov.list)), "AnnotatedDataFrame")
rownames(all.cov) <- colnames(t(all.X))
all.yvars <- do.call(rbind, y.vars)
whole_eset <- ExpressionSet(t(all.X), all.cov)


###################################################
### code chunk number 20: simulatorZ-vignette.Rnw:286-304
###################################################
lp.index <- c()
for(i in seq_along(esets)){
  lp.index <- c(lp.index, rep(i, length(y.vars[[i]][, 1])))
}
  
truemod <- getTrueModel(list(whole_eset), list(all.yvars), parstep=100)

simmodels <- simData(esets, 150, y.vars)
beta <- grid <- survH <- censH <- lp <- list()
for(listid in seq_along(esets)){
  beta[[listid]] <- truemod$beta[[1]]
  grid[[listid]] <- truemod$grid[[1]]
  survH[[listid]] <- truemod$survH[[1]]
  censH[[listid]] <- truemod$censH[[1]]
  lp[[listid]] <- truemod$lp[[1]][which(lp.index==listid)]
}
res <- list(beta=beta, grid=grid, survH=survH, censH=censH, lp=lp)
simmodels <- simTime(simmodels, y.vars, res)


###################################################
### code chunk number 21: simulatorZ-vignette.Rnw:323-325
###################################################
library(parathyroidSE)
data("parathyroidGenesSE")


###################################################
### code chunk number 22: simulatorZ-vignette.Rnw:331-334
###################################################
sim.sets <- simData(list(parathyroidGenesSE), 100)
names(sim.sets)
sim.sets$obj[[1]]   #The simulated RangedSummarizedExperiment


###################################################
### code chunk number 23: simulatorZ-vignette.Rnw:339-340
###################################################
table(colData(parathyroidGenesSE)$run)


###################################################
### code chunk number 24: simulatorZ-vignette.Rnw:345-346
###################################################
table(colData(sim.sets$obj[[1]])$run)


###################################################
### code chunk number 25: simulatorZ-vignette.Rnw:354-356
###################################################
z <- zmatrix(obj=esets, y.vars=y.list,
             fold=3,trainingFun=plusMinus)


###################################################
### code chunk number 26: simulatorZ-vignette.Rnw:358-359
###################################################
print(z)


###################################################
### code chunk number 27: simulatorZ-vignette.Rnw:362-378
###################################################
Z.list <- list()
CV <- CSV <- c()
for(b in 1:10){
  print(paste("iteration: ", b, sep=""))
  sim2.esets <- simBootstrap(obj=esets, n.samples=150, y.vars=y.list,
                             parstep=100, type="two-steps")
  Z.list[[b]] <- zmatrix(obj=sim2.esets$obj.list, 
                         y.vars=sim2.esets$y.vars.list, fold=4,
                         trainingFun=plusMinus)
  sum.cv <- 0
  for(i in seq_along(esets)){
    sum.cv <- sum.cv + Z.list[[b]][i, i]
  }
  CV[b] <- sum.cv / length(esets)
  CSV[b] <- (sum(Z.list[[b]]) - sum.cv) / (length(esets)*(length(esets)-1))
}


###################################################
### code chunk number 28: box
###################################################
average.Z <- Z.list[[1]]
for(i in 2:length(Z.list)){
  average.Z <- average.Z + Z.list[[i]]
}
average.Z <- average.Z / 10
print(average.Z)

resultlist <- list(CSV=CSV, CV=CV)
boxplot(resultlist, col=c("white", "grey"), ylab="C-Index", 
        boxwex = 0.25, xlim=c(0.5, 2.5))


###################################################
### code chunk number 29: simulatorZ-vignette.Rnw:396-397
###################################################
sessionInfo()