## @knitr env, include=FALSE, echo=FALSE, cache=FALSE
library("knitr")
opts_chunk$set(fig.align = 'center', 
               fig.show = 'hold', 
               par = TRUE,
               prompt = TRUE,
               eval = TRUE,
               comment = NA)
options(replace.assign = TRUE, 
        width = 55)

suppressPackageStartupMessages(library("MSnbase"))
suppressWarnings(suppressPackageStartupMessages(library("pRoloc")))
suppressPackageStartupMessages(library("pRolocdata"))
suppressPackageStartupMessages(library("class"))
set.seed(1)


## @knitr libraries
library("MSnbase")
library("pRoloc")
library("pRolocdata")


## @knitr setcols
setStockcol(paste0(getStockcol(), 70))


## @knitr readCsvData0
## The original data for replicate 1, available
## from the pRolocdata package
f0 <- dir(system.file("extdata", package = "pRolocdata"), 
          full.names = TRUE, pattern = "pr800866n_si_004-rep1.csv")
csv <- read.csv(f0)


## @knitr showOrgCsv
head(csv, n=3)


## @knitr readCsvData1, tidy = FALSE
## The quantitation data, from the original data
f1 <- dir(system.file("extdata", package = "pRoloc"), 
          full.names = TRUE, pattern = "exprsFile.csv")
exprsCsv <- read.csv(f1)
## Feature meta-data, from the original data
f2 <- dir(system.file("extdata", package = "pRoloc"), 
          full.names = TRUE, pattern = "fdataFile.csv")
fdataCsv <- read.csv(f2)
## Sample meta-data, a new file
f3 <- dir(system.file("extdata", package = "pRoloc"), 
          full.names = TRUE, pattern = "pdataFile.csv")
pdataCsv <- read.csv(f3)


## @knitr showExprsFile
head(exprsCsv, n=3)


## @knitr showFdFile
head(fdataCsv, n=3)


## @knitr showPdFile
pdataCsv


## @knitr makeMSnSet
tan2009r1 <- readMSnSet(exprsFile = f1,
                        featureDataFile = f2,
                        phenoDataFile = f3,
                        sep = ",")
tan2009r1


## @knitr showSubset
smallTan <- tan2009r1[1:5, 1:2]
dim(smallTan)
exprs(smallTan)


## @knitr rmtan, echo=FALSE
## remove manullay constructred data
rm(tan2009r1)
data(tan2009r1)


## @knitr loadTan1
data(tan2009r1)


## @knitr lookAtTan
getMarkers(tan2009r1)
getMarkers(tan2009r1, fcol = "PLSDA")


## @knitr realtiveQuants
summary(rowSums(exprs(tan2009r1)))


## @knitr norm, echo=TRUE, message=FALSE, cache=TRUE
## create a small illustrative test data
data(itraqdata)
itraqdata <- quantify(itraqdata, method = "trap", 
                      reporters = iTRAQ4, 
                      verbose = FALSE, parallel = FALSE)
## the quantification data
head(exprs(itraqdata), n = 3)
summary(rowSums(exprs(itraqdata)))
## normalising to the sum of feature intensitites
itraqnorm <- normalise(itraqdata, method = "sum")
processingData(itraqnorm)
head(exprs(itraqnorm), n = 3)
summary(rowSums(exprs(itraqnorm)))


## @knitr featurenames
head(featureNames(tan2009r1))


## @knitr showplotdist, echo=TRUE, eval=FALSE, prompt=FALSE
## ## indices of the mito markers
## j <- which(fData(tan2009r1)$markers == "mitochondrion")
## ## indices of all proteins assigned to the mito
## i <- which(fData(tan2009r1)$PLSDA == "mitochondrion")
## plotDist(tan2009r1[i, ],
##          markers = featureNames(tan2009r1)[j])


## @knitr plotdist1, dev='pdf', fig.width=9, fig.height=7, echo=FALSE
par(mfrow = c(2,2), mar = c(4, 4, 2, 0.1))
cls <- getStockcol()[1:4]
plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "mitochondrion"), ],
         markers = featureNames(tan2009r1)
         [which(fData(tan2009r1)$markers == "mitochondrion")],
         mcol = cls[1], 
         main = "mitochondrion")
plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "ER/Golgi"), ],
         markers = featureNames(tan2009r1)
         [which(fData(tan2009r1)$markers == "ER")],
         mcol = cls[2],
         main = "ER")
plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "ER/Golgi"), ],
         markers = featureNames(tan2009r1)
         [which(fData(tan2009r1)$markers == "Golgi")],
         mcol = cls[3],
         main = "Golgi")
plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "PM"), ],
         markers = featureNames(tan2009r1)
         [which(fData(tan2009r1)$markers == "PM")],
         mcol = cls[4],
         main = "PM")


## @knitr plot2d, dev='pdf', fig.width=5, fig.height=5, echo=TRUE
plot2D(tan2009r1, fcol = "PLSDA", fpch = "markers")
addLegend(tan2009r1, fcol = "PLSDA", 
          where = "bottomright", bty = "n", cex = .7)


## @knitr plot2dnull, dev='pdf', fig.width=5, fig.height=5, echo=TRUE
plot2D(tan2009r1, fcol = NULL)


## @knitr plotKmeans, dev='pdf', fig.width=8, fig.height=6, echo=TRUE
kcl <- MLearn( ~ ., tan2009r1,  kmeansI, centers=5)
plot(kcl, exprs(tan2009r1))


## @knitr plotHclust, dev='pdf', fig.width=8, fig.height=6, echo=TRUE, tidy = FALSE
hcl <- MLearn( ~ ., tan2009r1,  
              hclustI(distFun = dist, cutParm = list(k = 5)))
plot(hcl, labels = FALSE)


## @knitr plotPam, dev='pdf', fig.width=8, fig.height=6, echo=TRUE
pcl <- MLearn( ~ ., tan2009r1,  pamI(dist), k = 5)
plot(pcl, data = exprs(tan2009r1))


## @knitr svmParamOptim, cache = TRUE, message = FALSE
params <- svmOptimisation(tan2009r1, times = 10, xval = 5, verbose = FALSE)
params


## @knitr params, dev='pdf', fig.width=4, fig.height=4, echo=TRUE, out.width='.49\\linewidth'
plot(params)
levelPlot(params)


## @knitr f1count
f1Count(params)
getParams(params)


## @knitr svmRes0, warning=FALSE, eval = FALSE, tidy = FALSE
## ## manual setting of parameters
## svmres <- svmClassification(tan2009r1, sigma = 1, cost = 1)


## @knitr svmRes, warning=FALSE, tidy = FALSE
## using default best parameters
svmres <- svmClassification(tan2009r1, params)
processingData(svmres)
tail(fvarLabels(svmres), 4)


## @knitr getPredictions
p1 <- getPredictions(svmres, fcol = "svm")
minprob <- median(fData(svmres)$svm.scores)
p2 <- getPredictions(svmres, fcol = "svm", t = minprob)
table(p1, p2)


## @knitr predscoresPlot, eval=FALSE
## boxplot(svm.scores ~ svm,data = fData(svmres))
## abline(h = minprob)


## @knitr svmresfig, dev='pdf', fig.width=5, fig.height=5, echo=TRUE
ptsze <- exp(fData(svmres)$svm.scores) - 1
plot2D(svmres, fcol = "svm", fpch = "markers", cex = ptsze)
addLegend(svmres, fcol = "svm", 
          where = "bottomright", bty = "n", cex = .5)


## @knitr runPhenoDisco, eval=FALSE, warning=FALSE
## pdres <- phenoDisco(tan2009r1, GS = 10, times = 100, fcol = "PLSDA")


## @knitr loadpdres, echo=FALSE
fn <- dir(system.file("extdata", package = "pRoloc"), 
          full.names = TRUE, pattern = "pdres.rda")
load(fn)
rm(fn)


## @knitr phenoDiscoFvar
processingData(pdres)
tail(fvarLabels(pdres), 3)


## @knitr pdresfig, dev='pdf', fig.width=5, fig.height=5, echo=TRUE
plot2D(pdres, fcol = "pd")
addLegend(pdres, fcol = "pd", ncol = 2, 
          where = "bottomright", bty = "n", cex = .5)


## @knitr pdmarkers
getMarkers(tan2009r1, fcol = "pd.markers")


## @knitr pdsvm, cache = TRUE, message = FALSE, warning = FALSE, tidy = FALSE
w <- table(fData(tan2009r1)[, "pd.markers"])
(w <- 1/w[names(w) != "unknown"])
params <- svmOptimisation(tan2009r1, fcol = "pd.markers",                           
                          times = 10, xval = 5, 
                          class.weights = w,
                          verbose = FALSE)
tan2009r1 <- svmClassification(tan2009r1, params, 
                               class.weights = w,
                               fcol = "pd.markers")


## @knitr pdres2fig, dev='pdf', fig.width=6, fig.height=6, echo=TRUE
ptsze <- exp(fData(tan2009r1)$svm.scores) - 1
plot2D(tan2009r1, fcol = "svm", cex = ptsze)
addLegend(tan2009r1, fcol = "svm", where = "bottomright", ncol = 2, bty = "n", cex = .5)


## @knitr sessioninfo, results='asis', echo=FALSE
toLatex(sessionInfo())