inst/doc/ropls-vignette.R

In ropls: PCA, PLS(-DA) and OPLS(-DA) for multivariate analysis and feature selection of omics data

## ----global_options, include=FALSE--------------------------------------------
knitr::opts_chunk$set(fig.width = 6,
                      fig.height = 6,
                      fig.path = 'figures/')

## ----load, message = FALSE----------------------------------------------------
library(ropls)

## ----sacurine-----------------------------------------------------------------
data(sacurine)
names(sacurine)

## ----attach_code, message = FALSE---------------------------------------------
attach(sacurine)

## ----strF---------------------------------------------------------------------
view(dataMatrix)
view(sampleMetadata)
view(variableMetadata)

## ----pca_code, eval = FALSE---------------------------------------------------
#  sacurine.pca <- opls(dataMatrix)

## ----pca_result, echo = FALSE-------------------------------------------------
sacurine.pca <- opls(dataMatrix, fig.pdfC = "none")

## ----pca_figure, echo = FALSE, fig.show = 'hold'------------------------------
plot(sacurine.pca)

## ----pca-col------------------------------------------------------------------
genderFc <- sampleMetadata[, "gender"]
plot(sacurine.pca,
     typeVc = "x-score",
     parAsColFcVn = genderFc)

## ----pca-col-personalized-----------------------------------------------------
plot(sacurine.pca,
     typeVc = "x-score",
     parAsColFcVn = genderFc,
     parLabVc = as.character(sampleMetadata[, "age"]),
     parPaletteVc = c("green4", "magenta"))

## ----plsda--------------------------------------------------------------------
sacurine.plsda <- opls(dataMatrix, genderFc)

## ----oplsda-------------------------------------------------------------------
sacurine.oplsda <- opls(dataMatrix, genderFc,
                        predI = 1, orthoI = NA)

## ----oplsda_subset, warning=FALSE---------------------------------------------
sacurine.oplsda <- opls(dataMatrix, genderFc,
                        predI = 1, orthoI = NA,
                        subset = "odd")

## ----train--------------------------------------------------------------------
trainVi <- getSubsetVi(sacurine.oplsda)
table(genderFc[trainVi], fitted(sacurine.oplsda))

## ----test---------------------------------------------------------------------
table(genderFc[-trainVi],
      predict(sacurine.oplsda, dataMatrix[-trainVi, ]))

## ----overfit, echo = 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),
                   fig.pdfC = "none",
                   info.txtC = "none")
  plot(randPlsi, typeVc = "x-score",
       parCexN = 1.3, parTitleL = FALSE,
       parCexMetricN = 0.5)
  mtext(featI/obsI, font = 2, line = 2)
  if (featI == featMaxI)
    plot(randPlsi,
         typeVc = "permutation",
         parCexN = 1.3)
}
mtext(" obs./feat. ratio:",
      adj = 0, at = 0, font = 2,
      line = -2, outer = TRUE)

## ----vip----------------------------------------------------------------------
ageVn <- sampleMetadata[, "age"]

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

vipVn <- getVipVn(opls(dataMatrix, ageVn,
                       predI = 1, orthoI = NA,
                       fig.pdfC = "none"))

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

opar <- 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")

par(opar)

## ----expressionset_code, message = FALSE, warning=FALSE-----------------------
library(Biobase)
sacSet <- ExpressionSet(assayData = t(dataMatrix),
                        phenoData = new("AnnotatedDataFrame",
                                        data = sampleMetadata))
view(sacSet)
opls(sacSet, "gender", orthoI = NA)

## ----fromW4M------------------------------------------------------------------
sacSet <- fromW4M(file.path(path.package("ropls"), "extdata"))
sacSet

## ----toW4M, eval = FALSE------------------------------------------------------
#  toW4M(sacSet, paste0(getwd(), "/out_"))

## ----detach-------------------------------------------------------------------
detach(sacurine)

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()