### R code from vignette source 'curatedOvarianData_vignette.rnw'

###################################################
### code chunk number 1: <style-Sweave
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: preliminaries
###################################################
library(sva)
library(logging)
options(width=60)


###################################################
### code chunk number 3: example1tcgastep1
###################################################
library(curatedOvarianData)


###################################################
### code chunk number 4: example1tcgastep2 (eval = FALSE)
###################################################
## data(package="curatedOvarianData")


###################################################
### code chunk number 5: example1tcgastep2
###################################################
data(TCGA_eset)
TCGA_eset


###################################################
### code chunk number 6: example2loadstep1
###################################################
source(system.file("extdata", 
"patientselection.config",package="curatedOvarianData"))
ls()


###################################################
### code chunk number 7: showls
###################################################
#remove.samples and duplicates are too voluminous:
sapply(ls(), function(x) if(!x %in% c("remove.samples", "duplicates")) print(get(x)))


###################################################
### code chunk number 8: example2loadstep2
###################################################
source(system.file("extdata", "createEsetList.R", package =
"curatedOvarianData"))


###################################################
### code chunk number 9: example2loadstep3
###################################################
names(esets)


###################################################
### code chunk number 10: example3prepare
###################################################
esets[[1]]$y
forestplot <- function(esets, y="y", probeset, formula=y~probeset,
mlab="Overall", rma.method="FE", at=NULL,xlab="Hazard Ratio",...) {
    require(metafor)
    esets <- esets[sapply(esets, function(x) probeset %in% featureNames(x))]
    coefs <- sapply(1:length(esets), function(i) {
        tmp   <- as(phenoData(esets[[i]]), "data.frame")
        tmp$y <- esets[[i]][[y]]
        tmp$probeset <- exprs(esets[[i]])[probeset,]

        summary(coxph(formula,data=tmp))$coefficients[1,c(1,3)]
    })  

    res.rma <- metafor::rma(yi = coefs[1,], sei = coefs[2,], 
        method=rma.method)

    if (is.null(at)) at <- log(c(0.25,1,4,20))
    forest.rma(res.rma, xlab=xlab, slab=gsub("_eset$","",names(esets)),
    atransf=exp, at=at, mlab=mlab,...)
    return(res.rma)
}


###################################################
### code chunk number 11: example3plot
###################################################
res <- forestplot(esets=esets,probeset="CXCL12",at=log(c(0.5,1,2,4)))


###################################################
### code chunk number 12: filterdatasets
###################################################
idx.tumorstage <- sapply(esets, function(X) 
    sum(!is.na(X$tumorstage)) > 0 & length(unique(X$tumorstage)) > 1)

idx.debulking <- sapply(esets, function(X) 
    sum(X$debulking=="suboptimal",na.rm=TRUE)) > 0 


###################################################
### code chunk number 13: example4plot
###################################################
res <- forestplot(esets=esets[idx.debulking & idx.tumorstage],
    probeset="CXCL12",formula=y~probeset+debulking+tumorstage, 
    at=log(c(0.5,1,2,4)))


###################################################
### code chunk number 14: example5plot
###################################################
res <- forestplot(esets=esets,probeset="CXCR4",at=log(c(0.5,1,2,4)))


###################################################
### code chunk number 15: combine2
###################################################
combine2 <- function(X1, X2) {
    fids <- intersect(featureNames(X1), featureNames(X2))
    X1 <- X1[fids,]
    X2 <- X2[fids,]
    ExpressionSet(cbind(exprs(X1),exprs(X2)), 
        AnnotatedDataFrame(rbind(as(phenoData(X1),"data.frame"), 
                                 as(phenoData(X2),"data.frame")))
    )
}


###################################################
### code chunk number 16: boxplot1
###################################################
data(E.MTAB.386_eset)
data(GSE30161_eset)
X <- combine2(E.MTAB.386_eset, GSE30161_eset)
boxplot(exprs(X))


###################################################
### code chunk number 17: combat
###################################################
mod <- model.matrix(~as.factor(tumorstage), data=X)
batch <- as.factor(grepl("DFCI",sampleNames(X)))
combat_edata <- ComBat(dat=exprs(X), batch=batch, mod=mod)


###################################################
### code chunk number 18: boxplot2
###################################################
boxplot(combat_edata)


###################################################
### code chunk number 19: expand
###################################################
expandProbesets <- function (eset, sep = "///") 
{
    x <- lapply(featureNames(eset), function(x) strsplit(x, sep)[[1]])
    eset <- eset[order(sapply(x, length)), ]
    x <- lapply(featureNames(eset), function(x) strsplit(x, sep)[[1]])
    idx <- unlist(sapply(1:length(x), function(i) rep(i, length(x[[i]]))))
    xx <- !duplicated(unlist(x))
    idx <- idx[xx]
    x <- unlist(x)[xx]
    eset <- eset[idx, ]
    featureNames(eset) <- x
    eset
}

X <- TCGA_eset[head(grep("///", featureNames(TCGA_eset))),]
exprs(X)[,1:3]
exprs(expandProbesets(X))[,1:3]


###################################################
### code chunk number 20: featureData
###################################################
head(pData(featureData(GSE18520_eset)))


###################################################
### code chunk number 21: annotationeg
###################################################
annotation(GSE18520_eset)


###################################################
### code chunk number 22: loadallsamples
###################################################
rm(list=ls())
source(system.file("extdata", 
    "patientselection_all.config",package="curatedOvarianData"))
source(system.file("extdata", "createEsetList.R", package =
    "curatedOvarianData"))


###################################################
### code chunk number 23: heatmap
###################################################
.esetsStats <- function(esets) {
    res <- lapply(varLabels(esets[[1]]), function(covar) unlist(sapply(esets, 
        function(X) sum(!is.na(X[[covar]]))>0)))
    names(res) <- varLabels(esets[[1]])    
    do.call(rbind, res)
}

df.r <- .esetsStats(esets)
M <- as.matrix(apply(df.r,c(1,2),ifelse,0,1))
colnames(M) <- gsub("_eset$", "", colnames(M))
# no need to show the sample ids
M <- M[-(1:2),]
heatmap(M[nrow(M):1,],scale="none",margins=c(8,10),Rowv=NA)


###################################################
### code chunk number 24: esetToTableFuns
###################################################
source(system.file("extdata", "summarizeEsets.R", package =
    "curatedOvarianData"))


###################################################
### code chunk number 25: esettable
###################################################
summary.table <- t(sapply(esets, getEsetData))
rownames(summary.table) <- sub("_eset", "", rownames(summary.table))


###################################################
### code chunk number 26: writeesettable
###################################################
(myfile <- tempfile())
write.table(summary.table, file=myfile, row.names=FALSE, quote=TRUE, sep=",")


###################################################
### code chunk number 27: xtable
###################################################
library(xtable)
print(xtable(summary.table[, c(2, 3, 4, 5, 7)], 
             caption="Datasets provided by curatedOvarianData.  This is an abbreviated version of Table 1 of the manuscript; the full version is written by the write.table command above.  Stage column is early/late/unknown, histology column is ser/clearcell/endo/mucinous/other/unknown.", 
             table.placement="p", caption.placement="bottom"),
      floating.environment='sidewaystable')


###################################################
### code chunk number 28: simplygetdata (eval = FALSE)
###################################################
## library(curatedOvarianData)
## data(GSE30161_eset)
## write.csv(exprs(GSE30161_eset), file="GSE30161_eset_exprs.csv")
## write.csv(pData(GSE30161_eset), file="GSE30161_eset_clindata.csv")


###################################################
### code chunk number 29: simplyseveraldatasets (eval = FALSE)
###################################################
## data.to.fetch <- c("GSE30161_eset", "E.MTAB.386_eset")
## for (onedata in data.to.fetch){
##     print(paste("Fetching", onedata))
##     data(list=onedata)
##     write.csv(exprs(get(onedata)), file=paste(onedata, "_exprs.csv", sep=""))
##     write.csv(pData(get(onedata)), file=paste(onedata, "_clindata.csv", sep=""))
## }


###################################################
### code chunk number 30: sessioninfo
###################################################
toLatex(sessionInfo())