### R code from vignette source 'ropls.Rnw'

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


###################################################
### code chunk number 2: ropls.Rnw:86-87
###################################################
library(ropls)


###################################################
### code chunk number 3: ropls.Rnw:98-100
###################################################
data(sacurine)
names(sacurine)


###################################################
### code chunk number 4: ropls.Rnw:105-109
###################################################
attach(sacurine)
strF(dataMatrix)
strF(sampleMetadata)
strF(variableMetadata)


###################################################
### code chunk number 5: ropls.Rnw:117-118 (eval = FALSE)
###################################################
## sacurine.pca <- opls(dataMatrix)


###################################################
### code chunk number 6: ropls.Rnw:121-122
###################################################
sacurine.pca <- opls(dataMatrix, plotL = FALSE)


###################################################
### code chunk number 7: figPca
###################################################
layout(matrix(1:4, nrow = 2, byrow = TRUE))
for(typeC in c("overview", "outlier", "x-score", "x-loading"))
plot(sacurine.pca, typeVc = typeC, parDevNewL = FALSE)


###################################################
### code chunk number 8: ropls.Rnw:147-150 (eval = FALSE)
###################################################
## genderFc <- sampleMetadata[, "gender"]
## plot(sacurine.pca, typeVc = "x-score",
## parAsColFcVn = genderFc, parEllipsesL = TRUE)


###################################################
### code chunk number 9: figPcaMah
###################################################
genderFc <- sampleMetadata[, "gender"]
plot(sacurine.pca, typeVc = "x-score",
parAsColFcVn = genderFc, parEllipsesL = TRUE, parDevNewL = FALSE)


###################################################
### code chunk number 10: ropls.Rnw:173-174 (eval = FALSE)
###################################################
## sacurine.plsda <- opls(dataMatrix, genderFc)


###################################################
### code chunk number 11: ropls.Rnw:176-177
###################################################
sacurine.plsda <- opls(dataMatrix, genderFc, plotL = FALSE)


###################################################
### code chunk number 12: figPls
###################################################
layout(matrix(1:4, nrow = 2, byrow = TRUE))
for(typeC in c("permutation", "overview", "outlier", "x-score"))
plot(sacurine.plsda, typeVc = typeC, parDevNewL = FALSE)


###################################################
### code chunk number 13: ropls.Rnw:209-211 (eval = FALSE)
###################################################
## sacurine.oplsda <- opls(dataMatrix, genderFc,
## predI = 1, orthoI = NA)


###################################################
### code chunk number 14: ropls.Rnw:213-215
###################################################
sacurine.oplsda <- opls(dataMatrix, genderFc,
predI = 1, orthoI = NA, plotL = FALSE)


###################################################
### code chunk number 15: figOpl
###################################################
layout(matrix(1:4, nrow = 2, byrow = TRUE))
for(typeC in c("permutation", "overview", "outlier", "x-score"))
plot(sacurine.oplsda, typeVc = typeC, parDevNewL = FALSE)


###################################################
### code chunk number 16: ropls.Rnw:232-233 (eval = FALSE)
###################################################
## sacurine.oplsda <- opls(dataMatrix, genderFc, predI = 1, orthoI = NA, subset = "odd")


###################################################
### code chunk number 17: ropls.Rnw:235-236
###################################################
sacurine.oplsda <- opls(dataMatrix, genderFc, predI = 1, orthoI = NA, permI = 0, subset = "odd", plotL = FALSE)


###################################################
### code chunk number 18: ropls.Rnw:241-243
###################################################
trainVi <- getSubsetVi(sacurine.oplsda)
table(genderFc[trainVi], fitted(sacurine.oplsda))


###################################################
### code chunk number 19: ropls.Rnw:248-250
###################################################
table(genderFc[-trainVi],
      predict(sacurine.oplsda, dataMatrix[-trainVi, ]))


###################################################
### code chunk number 20: figOverfit
###################################################
set.seed(123)
obsI <- 20
featVi <- c(2, 20, 200)
featMaxI <- max(featVi)
xRandMN <- matrix(runif(obsI * featMaxI), nrow = obsI)
yRandVn <- sample(c(rep(0, obsI / 2), rep(1, obsI / 2)))

layout(matrix(1:4, nrow = 2, byrow = TRUE))

for(featI in featVi) {
    randPlsi <- opls(xRandMN[, 1:featI], yRandVn,
                  predI = 2,
                  permI = ifelse(featI == featMaxI, 100, 0),
                  plotL = FALSE)
    plot(randPlsi, typeVc = "x-score", parDevNewL = FALSE,
         parCexN = 1.3, parTitleL = FALSE)
    mtext(featI/obsI, font = 2, line = 2)
    if(featI == featMaxI)
         plot(randPlsi, typeVc = "permutation", parDevNewL = FALSE,
           parCexN = 1.3)
    }
mtext(" obs./feat. ratio:", adj = 0, at = 0, font = 2, line = -2, outer = TRUE)


###################################################
### code chunk number 21: ropls.Rnw:293-315 (eval = FALSE)
###################################################
## set.seed(123)
## obsI <- 20
## featVi <- c(2, 20, 200)
## featMaxI <- max(featVi)
## xRandMN <- matrix(runif(obsI * featMaxI), nrow = obsI)
## yRandVn <- sample(c(rep(0, obsI / 2), rep(1, obsI / 2)))
## 
## layout(matrix(1:4, nrow = 2, byrow = TRUE))
## 
## for(featI in featVi) {
##     randPlsi <- opls(xRandMN[, 1:featI], yRandVn,
##                   predI = 2,
##                   permI = ifelse(featI == featMaxI, 100, 0),
##                   plotL = FALSE)
##     plot(randPlsi, typeVc = "x-score", parDevNewL = FALSE,
##          parCexN = 1.3, parTitleL = FALSE)
##     mtext(featI/obsI, font = 2, line = 2)
##     if(featI == featMaxI)
##          plot(randPlsi, typeVc = "permutation", parDevNewL = FALSE,
##            parCexN = 1.3)
##     }
## mtext(" obs./feat. ratio:", adj = 0, at = 0, font = 2, line = -2, outer = TRUE)


###################################################
### code chunk number 22: figVip
###################################################
ageVn <- sampleMetadata[, "age"]

pvaVn <- apply(dataMatrix, 2,
               function(feaVn) cor.test(ageVn, feaVn)[["p.value"]])

vipVn <- getVipVn(opls(dataMatrix, ageVn, predI = 1, orthoI = NA, plot = FALSE))

quantVn <- qnorm(1 - pvaVn / 2)
rmsQuantN <- sqrt(mean(quantVn^2))

par(font = 2, font.axis = 2, font.lab = 2, las = 1,
    mar = c(5.1, 4.6, 4.1, 2.1),
    lwd = 2, pch = 16)

plot(pvaVn, vipVn,
     col = "red",
     pch = 16,
     xlab = "p-value", ylab = "VIP", xaxs = "i", yaxs = "i")

box(lwd = 2)

curve(qnorm(1 - x / 2) / rmsQuantN, 0, 1, add = TRUE, col = "red", lwd = 3)

abline(h = 1, col = "blue")
abline(v = 0.05, col = "blue")


###################################################
### code chunk number 23: ropls.Rnw:361-386 (eval = FALSE)
###################################################
## ageVn <- sampleMetadata[, "age"]
## 
## pvaVn <- apply(dataMatrix, 2,
##                function(feaVn) cor.test(ageVn, feaVn)[["p.value"]])
## 
## vipVn <- getVipVn(opls(dataMatrix, ageVn, predI = 1, orthoI = NA, plot = FALSE))
## 
## quantVn <- qnorm(1 - pvaVn / 2)
## rmsQuantN <- sqrt(mean(quantVn^2))
## 
## par(font = 2, font.axis = 2, font.lab = 2, las = 1,
##     mar = c(5.1, 4.6, 4.1, 2.1),
##     lwd = 2, pch = 16)
## 
## plot(pvaVn, vipVn,
##      col = "red",
##      pch = 16,
##      xlab = "p-value", ylab = "VIP", xaxs = "i", yaxs = "i")
## 
## box(lwd = 2)
## 
## curve(qnorm(1 - x / 2) / rmsQuantN, 0, 1, add = TRUE, col = "red", lwd = 3)
## 
## abline(h = 1, col = "blue")
## abline(v = 0.05, col = "blue")


###################################################
### code chunk number 24: ropls.Rnw:425-426
###################################################
detach(sacurine)


###################################################
### code chunk number 25: ropls.Rnw:432-433 (eval = FALSE)
###################################################
## sessionInfo()


###################################################
### code chunk number 26: ropls.Rnw:435-436
###################################################
toLatex(sessionInfo())


###################################################
### code chunk number 27: ropls.Rnw:446-450
###################################################
library(faahKO)
cdfpath <- system.file("cdf", package = "faahKO")
cdffiles <- list.files(cdfpath, recursive = TRUE, full.names = TRUE)
basename(cdffiles)


###################################################
### code chunk number 28: ropls.Rnw:455-456
###################################################
library(xcms)


###################################################
### code chunk number 29: ropls.Rnw:458-459
###################################################
xset <- xcmsSet(cdffiles)


###################################################
### code chunk number 30: ropls.Rnw:461-463
###################################################
xset
xset <- group(xset)


###################################################
### code chunk number 31: ropls.Rnw:465-466 (eval = FALSE)
###################################################
## xset2 <- retcor(xset, family = "symmetric", plottype = "mdevden")


###################################################
### code chunk number 32: ropls.Rnw:468-469
###################################################
xset2 <- retcor(xset, family = "symmetric", plottype = "none")


###################################################
### code chunk number 33: ropls.Rnw:471-472
###################################################
xset2 <- group(xset2, bw = 10)


###################################################
### code chunk number 34: ropls.Rnw:474-475
###################################################
xset3 <- fillPeaks(xset2)


###################################################
### code chunk number 35: ropls.Rnw:480-483
###################################################
library(CAMERA)
diffreport <- annotateDiffreport(xset3, quick=TRUE)
diffreport[1:4, ]


###################################################
### code chunk number 36: ropls.Rnw:488-495
###################################################
sampleVc <- grep("^ko|^wt", colnames(diffreport), value = TRUE)
dataMatrix <- t(as.matrix(diffreport[, sampleVc]))
dimnames(dataMatrix) <- list(sampleVc, diffreport[, "name"])
sampleMetadata <- data.frame(row.names = sampleVc,
genotypeFc = substr(sampleVc, 1, 2))
variableMetadata <- diffreport[, !(colnames(diffreport) %in% c("name", sampleVc))]
rownames(variableMetadata) <- diffreport[, "name"]


###################################################
### code chunk number 37: ropls.Rnw:500-503
###################################################
library(ropls)
opls(dataMatrix)
opls(dataMatrix, sampleMetadata[, "genotypeFc"], orthoI = NA)


###################################################
### code chunk number 38: ropls.Rnw:506-507
###################################################
rm(list = ls())