# mixOmics In mixOmics: Omics Data Integration Project

BiocStyle::markdown()

library(knitr)
# global options
knitr::opts_chunkset(dpi = 100, echo=TRUE, warning=FALSE, message=FALSE, eval = TRUE, fig.show=TRUE, fig.width= 7,fig.height= 6,fig.align='center', out.width = '80%'#, #fig.path= 'Figures/' ) # knitr::opts_chunkset(
#   collapse = TRUE,
#   comment = "#>"
# )


# Bookdown version

This vignette is also available as a bookdown format on our gitHub page https://mixomicsteam.github.io/Bookdown/.

# Introduction {#intro}

mixOmics is an R toolkit dedicated to the exploration and integration of biological data sets with a specific focus on variable selection. The package currently includes nineteen multivariate methodologies, mostly developed by the mixOmics team (see some of our references in \@ref(intro:pubs)). Originally, all methods were designed for omics data, however, their application is not limited to biological data only. Other applications where integration is required can be considered, but mostly for the case where the predictor variables are continuous (see also \@ref(intro:datatypes)). In mixOmics, a strong focus is given to graphical representation to better translate and understand the relationships between the different data types and visualize the correlation structure at both sample and variable levels.

## Input data {#intro:datatypes}

Note the data pre-processing requirements before analysing data with mixOmics:

• Types of data. Different types of biological data can be explored and integrated with mixOmics. Our methods can handle molecular features measured on a continuous scale (e.g. microarray, mass spectrometry-based proteomics and metabolomics) or sequenced-based count data (RNA-seq, 16S, shotgun metagenomics) that become continuous' data after pre-processing and normalisation.

• Normalisation. The package does not handle normalisation as it is platform specific and we cover a too wide variety of data! Prior to the analysis, we assume the data sets have been normalised using appropriate normalisation methods and pre-processed when applicable.

• Prefiltering. While mixOmics methods can handle large data sets (several tens of thousands of predictors), we recommend pre-filtering the data to less than 10K predictor variables per data set, for example by using Median Absolute Deviation [@Ten16] for RNA-seq data, by removing consistently low counts in microbiome data sets [@Lec16] or by removing near zero variance predictors. Such step aims to lessen the computational time during the parameter tuning process.

• Data format. Our methods use matrix decomposition techniques. Therefore, the numeric data matrix or data frames have $n$ observations or samples in rows and $p$ predictors or variables (e.g. genes, proteins, OTUs) in columns.

• Covariates. In the current version of mixOmics, covariates that may confound the analysis are not included in the methods. We recommend correcting for those covariates beforehand using appropriate univariate or multivariate methods for batch effect removal. Contact us for more details as we are currently working on this aspect.

## Methods

### Some background knowledge {#intro:background}

We list here the main methodological or theoretical concepts you need to know to be able to efficiently apply mixOmics:

• Individuals, observations or samples: the experimental units on which information are collected, e.g. patients, cell lines, cells, faecal samples ...

• Variables, predictors: read-out measured on each sample, e.g. gene (expression), protein or OTU (abundance), weight ...

• Variance: measures the spread of one variable. In our methods we estimate the variance of components rather that variable read-outs. A high variance indicates that the data points are very spread out from the mean, and from one another (scattered).

• Covariance: measures the strength of the relationship between two variables, i.e whether they co-vary. A high covariance value indicates a strong relationship, e.g weight and height in individuals frequently vary roughly in the same way; roughly, the heaviest are the tallest. A covariance value has no lower or upper bound.

• Correlation: a standardized version of the covariance that is bounded by -1 and 1.

• Linear combination: variables are combined by multiplying each of of them by a coefficient and adding the results. A linear combination of height and weight could be 2 $$weight - 1.5$$ height with the coefficients 2 and -1.5 assigned with weight and height respectively.

• Component: an artificial variable built from a linear combination of the observed variables in a given data set. Variable coefficients are optimally defined based on some statistical criterion. For example in Principal Component Analysis, the coefficients in the (principal) component is defined so as to maximise the variance of the component.

• Sample plot: representation of the samples projected in a small space spanned (defined) by the components. Samples coordinates are determined by their components values, or scores.

• Correlation circle plot: representation of the variables in a space spanned by the components. Each variable coordinate is defined as the correlation between the original variable value and each component. A correlation circle plot enables to visualise the correlation between variables - negative or positive correlation, defined by the cosine angle between the centre of the circle and each variable point) and the contribution of each variable to each component - defined by absolute value of the coordinate on each component. For this interpretation, data need to be centred and scaled (by default in most of our methods except PCA). For more details on this insightful graphic, see Figure 1 in [@Gon12].

• Unsupervised analysis: the method does not take into account any known sample groups and the analysis is exploratory. Examples of unsupervised methods covered in this vignette are Principal Component Analysis (PCA, Chapter \@ref(pca)), Projection to Latent Structures (PLS, Chapter \@ref(pls)), and also Canonical Correlation Analysis (CCA, not covered here).

• Supervised analysis: the method includes a vector indicating the class membership of each sample. The aim is to discriminate sample groups and perform sample class prediction. Examples of supervised methods covered in this vignette are PLS Discriminant Analysis (PLS-DA, Chapter \@ref(plsda)), DIABLO (Chapter \@ref(diablo)) and also MINT (not covered here [@Roh16]).

### Overview {#intro:overview}

Here is an overview of the most widely used methods in mixOmics that will be further detailed in this vignette, with the exception of rCCA and MINT. We depict them along with the type of data set they can handle.

library(mixOmics)
coul <- color.mixo(1:3)

plot(0, type="n", xlim=c(0,100), ylim=c(0,100), axes=FALSE,
xlab="",ylab="", main="mixOmics overview")
box()

# PCA
rect(xleft = 20, ybottom = 75, xright = 40, ytop = 95, col=coul[1])
text(5, 85, "PCA")

# PLS-DA
rect(xleft = 20, ybottom = 50, xright = 40, ytop = 70, col=coul[1])
rect(xleft = 43, ybottom = 50, xright = 45, ytop = 70, col=coul[2])
text(5, 60, "PLS-DA")

# PLS
rect(xleft = 20, ybottom = 25, xright = 40, ytop = 45, col=coul[1])
rect(xleft = 43, ybottom = 25, xright = 60, ytop = 45, col=coul[1])
text(5, 35, "PLS")

# DIABLO
rect(xleft = 20, ybottom = 0, xright = 40, ytop = 20, col=coul[1])
rect(xleft = 43, ybottom = 0, xright = 60, ytop = 20, col=coul[1])
points(x=61, y=10, pch=16, col=coul[3], cex=0.5)
points(x=62.5, y=10, pch=16, col=coul[3], cex=0.5)
points(x=64, y=10, pch=16, col=coul[3], cex=0.5)
rect(xleft = 65, ybottom = 0, xright = 80, ytop = 20, col=coul[1])
rect(xleft = 85, ybottom = 0, xright = 88, ytop = 20, col=coul[2])
text(5, 10, "DIABLO")

# legend
rect(xleft = 75, ybottom = 95, xright = 77, ytop = 97, col=coul[1])
text(90, 96, "Quantitative", cex=0.75)
rect(xleft = 75, ybottom = 90, xright = 77, ytop = 92, col=coul[2])
text(90, 91, "Qualitative", cex=0.75)

knitr::include_graphics("XtraFigs/Methods.png")

knitr::include_graphics("XtraFigs/cheatsheet.png")


### Key publications {#intro:pubs}

#### Methods papers

The methods implemented in mixOmics are described in detail in the following publications. A more extensive list can be found at this link.

#### Biological application papers

For those interested in how these methods have been applied to omics biological problems:

## Outline of this Vignette

• Chapter \@ref(start) details some practical aspects to get started
• Chapter \@ref(pca): Principal Components Analysis (PCA)
• Chapter \@ref(plsda): Projection to Latent Structure - Discriminant Analysis (PLS-DA)
• Chapter \@ref(pls): Projection to Latent Structures (PLS)
• Chapter \@ref(diablo): Integrative analysis for multiple data sets (DIABLO)

Each of the methods chapter has the following outline:

1. Type of biological question to be answered
2. Brief description of an illustrative data set
3. Principle of the method
4. Quick start of the method with the main functions and arguments
5. To go further: customized plots, additional graphical outputs and tuning parameters
6. FAQ

## Other methods not covered in this vignette

Other methods not covered in this document are described on our website and the following references:

# Let's get started {#start}

## Installation

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("mixOmics")


Alternatively, you can install the latest development version of the package from Github.

install.packages("devtools")
library(devtools)
install_github("mixOmicsTeam/mixOmics")


The mixOmics package will directly import the following packages: r Rpackage("igraph"), r Rpackage("rgl"), r Rpackage("ellipse"), r Rpackage("corpcor"), r Rpackage("RColorBrewer"), r Rpackage("parallel"), r Rpackage("dplyr"), r Rpackage("tidyr"), r Rpackage("reshape2"), r Rpackage("methods"), r Rpackage("matrixStats"), r Rpackage("rARPACK"), r Rpackage("gridExtra").

For apple mac users, if you are unable to install the imported library rgl, you will need to install the XQuartz software first.

library(mixOmics)


Check that there is no error when loading the package, especially for the rgl library (see above).

The examples we give in this vignette use data that are already part of the package, which means that the objects are stored as list. However, you do not need to store your data as list, just data frames will do. To upload your own data, check first that your working directory is set, then read your data from a .txt or .csv format, either by using File > Import Dataset in RStudio (although be mindful this option may pose problem with row and col name headers) or via one of these command lines:

# from csv file

# from txt file

# then in the argument in the mixOmics functions, just fill with:
# X = data


For more details about the arguments used to modify those functions, type ?read.csv or ?read.table in the R console. In this vignette we assume you are sufficiently proficient in R to input your data to carry on with the mixOmics analyses!

## Quick start in mixOmics {#start:PCA}

Each analysis should follow this workflow:

1. Run the method
2. Graphical representation of the samples
3. Graphical representation of the variables

Then use your critical thinking and additional functions and visual tools to make sense of your data! (some of which are listed in \@ref(intro:overview)) and will be described in the next Chapters.

For instance, for Principal Components Analysis, we first load the data:

data(nutrimouse)

## Principle of PCA

The aim of PCA [@Jol05] is to reduce the dimensionality of the data whilst retaining as much information as possible. 'Information' is referred here as variance. The idea is to create uncorrelated artificial variables called principal components (PCs) that combine in a linear manner the original (possibly correlated) variables (e.g. genes, metabolites, etc.).

Dimension reduction is achieved by projecting the data into the space spanned by the principal components (PC). In practice, it means that each sample is assigned a score on each new PC dimension - this score is calculated as a linear combination of the original variables to which a weight is applied. The weights of each of the original variables are stored in the so-called loading vectors associated to each PC. The dimension of the data is reduced by projecting the data into the smaller subspace spanned by the PCs, while capturing the largest sources of variation between samples.

The principal components are obtained so that their variance is maximised. To that end, we calculate the eigenvectors/eigenvalues of the variance-covariance matrix, often via singular value decomposition when the number of variables is very large. The data are usually centred (center = TRUE), and sometimes scaled (scale = TRUE) in the method. The latter is especially advised in the case where the variance is not homogeneous across variables.

The first PC is defined as the linear combination of the original variables that explains the greatest amount of variation. The second PC is then defined as the linear combination of the original variables that accounts for the greatest amount of the remaining variation subject of being orthogonal (uncorrelated) to the first component. Subsequent components are defined likewise for the other PCA dimensions. The user must therefore report how much information is explained by the first PCs as these are used to graphically represent the PCA outputs.

library(mixOmics)
data(liver.toxicity)
X <- liver.toxicity$gene  ## Quick start MyResult.pca <- pca(X) # 1 Run the method plotIndiv(MyResult.pca) # 2 Plot the samples plotVar(MyResult.pca) # 3 Plot the variables  If you were to run pca with this minimal code, you would be using the following default values: • ncomp =2: the first two principal components are calculated and are used for graphical outputs; • center = TRUE: data are centred (mean = 0) • scale = FALSE: data are not scaled. If scale = TRUE standardizes each variable (variance = 1). Other arguments can also be chosen, see ?pca. This example was shown in Chapter \@ref(start:PCA). The two plots are not extremely meaningful as specific sample patterns should be further investigated and the variable correlation circle plot contains too many variables to be easily interpreted. Let's improve those graphics as shown below to improve interpretation. ## To go further ### Customize plots Plots can be customized using numerous options in plotIndiv and plotVar. For instance, even if PCA does not take into account any information regarding the known group membership of each sample, we can include such information on the sample plot to visualize any natural' cluster that may corresponds to biological conditions. Here is an example where we include the sample groups information with the argument group: plotIndiv(MyResult.pca, group = liver.toxicity$treatment$Dose.Group, legend = TRUE)  Additionally, two factors can be displayed using both colours (argument group) and symbols (argument pch). For example here we display both Dose and Time of exposure and improve the title and legend: plotIndiv(MyResult.pca, ind.names = FALSE, group = liver.toxicity$treatment$Dose.Group, pch = as.factor(liver.toxicity$treatment$Time.Group), legend = TRUE, title = 'Liver toxicity: genes, PCA comp 1 - 2', legend.title = 'Dose', legend.title.pch = 'Exposure')  By including information related to the dose of acetaminophen and time of exposure enables us to see a cluster of low dose samples (blue and orange, top left at 50 and 100mg respectively), whereas samples with high doses (1500 and 2000mg in grey and green respectively) are more scattered, but highlight an exposure effect. To display the results on other components, we can change the comp argument provided we have requested enough components to be calculated. Here is our second PCA with 3 components: MyResult.pca2 <- pca(X, ncomp = 3) plotIndiv(MyResult.pca2, comp = c(1,3), legend = TRUE, group = liver.toxicity$treatment$Time.Group, title = 'Multidrug transporter, PCA comp 1 - 3')  Here, the 3rd component on the y-axis clearly highlights a time of exposure effect. ### Amount of variance explained and choice of number of components The amount of variance explained can be extracted with the following: a screeplot or the actual numerical proportions of explained variance, and cumulative proportion. plot(MyResult.pca2) MyResult.pca2  There are no clear guidelines on how many components should be included in PCA: it is data dependent and level of noise dependent. We often look at the elbow' on the screeplot above as an indicator that the addition of PCs does not drastically contribute to explain the remainder variance. ### Other useful plots We can also have a look at the variable coefficients in each component with the loading vectors. The loading weights are represented in decreasing order from bottom to top in plotLoadings. Their absolute value indicates the importance of each variable to define each PC, as represented by the length of each bar. See ?plotLoadings to change the arguments. # a minimal example plotLoadings(MyResult.pca)  # a customized example to only show the top 100 genes # and their gene name plotLoadings(MyResult.pca, ndisplay = 100, name.var = liver.toxicity$gene.ID[, "geneBank"],
size.name = rel(0.3))


Such representation will be more informative once we select a few variables in the next section \@ref(sPCA).

Plots can also be displayed in 3 dimensions using the option style="3d", and interactively (we use the rgl package for this).

plotIndiv(MyResult.pca2,
group = liver.toxicity$treatment$Dose.Group, style="3d",
legend = TRUE, title = 'Liver toxicity: genes, PCA comp 1 - 2 - 3')


## Variable selection with sparse PCA {#sPCA}

### Biological question

I would like to apply PCA but also be able to identify the key variables that contribute to the explanation of most variance in the data set.

Variable selection can be performed using the sparse version of PCA implemented in spca [@She08]. The user needs to provide the number of variables to select on each PC. Here for example we ask to select the top 15 genes contributing to the definition of PC1, the top 10 genes contributing to PC2 and the top 5 genes for PC3 (keepX=c(15,10,5)).

MyResult.spca <- spca(X, ncomp = 3, keepX = c(15,10,5))                 # 1 Run the method
plotIndiv(MyResult.spca, group = liver.toxicity$treatment$Dose.Group,   # 2 Plot the samples
pch = as.factor(liver.toxicity$treatment$Time.Group),
legend = TRUE, title = 'Liver toxicity: genes, sPCA comp 1 - 2',
legend.title = 'Dose', legend.title.pch = 'Exposure')
plotVar(MyResult.spca, cex = 1)                                        # 3 Plot the variables
# cex is used to reduce the size of the labels on the plot


Selected variables can be identified on each component with the selectVar function. Here the coefficient values are extracted, but there are other outputs, see ?selectVar:

selectVar(MyResult.spca, comp = 1)$value  Those values correspond to the loading weights that are used to define each component. A large absolute value indicates the importance of the variable in this PC. Selected variables are ranked from the most important (top) to the least important. We can complement this output with plotLoadings. We can see here that all coefficients are negative. plotLoadings(MyResult.spca)  If we look at component two, we can see a mix of positive and negative weights (also see in the plotVar), those correspond to variables that oppose the low and high doses (see from the plotIndiv): selectVar(MyResult.spca, comp=2)$value


## Tuning parameters

For this set of methods, two parameters need to be chosen:

• The number of components to retain,
• The number of variables to select on each component for sparse PCA.

The function tune.pca calculates the percentage of variance explained for each component, up to the minimum between the number of rows, or column in the data set. The optimal' number of components can be identified if an elbow appears on the screeplot. In the example below the cut-off is not very clear, we could choose 2 components.

tune.pca(X)


Regarding the number of variables to select in sparse PCA, there is not clear criterion at this stage. As PCA is an exploration method, we prefer to set arbitrary thresholds that will pinpoint the key variables to focus on during the interpretation stage.

Additional examples are provided in example(pca) and in our case studies on our website in the Methods and Case studies sections.

## FAQ

• Should I scale my data before running PCA? (scale = TRUE in pca)

• Without scaling: a variable with high variance will solely drive the first principal component
• With scaling: one noisy variable with low variability will be assigned the same variance as other meaningful variables
• Can I perform PCA with missing values?

• NIPALS (Non-linear Iterative PArtial Least Squares, implemented in mixOmics) can impute missing values but must be built on many components. The proportion of NAs should not exceed 20% of total data.
• When should I apply a multilevel approach in PCA? (multilevel argument in PCA)

• When the unique individuals are measured more than once (repeated measures)
• When the individual variation is less than treatment or time variation. This means that samples from each unique individual will tend to cluster rather than the treatments.
• When a multilevel vs no multilevel seems to visually make a difference on a PCA plot
• More details in this case study
• When should I apply a CLR transformation in PCA? (logratio = 'CLR' argument in PCA)

• When data are compositional, i.e. expressed as relative proportions. This is usually the case with microbiome studies as a result of pre-processing and normalisation, see more details here and in our case studies in the same tab.

# PLS - Discriminant Analysis (PLS-DA) {#plsda}

library(mixOmics)
coul <- color.mixo(1:3)

plot(0, type="n", xlim=c(0,100), ylim=c(83,97), axes=FALSE,
xlab="",ylab="", main="PLSDA overview")
box()

# PLS-DA
rect(xleft = 20, ybottom = 85, xright = 60, ytop = 95, col=coul[1])
rect(xleft = 63, ybottom = 85, xright = 65, ytop = 95, col=coul[2])
text(5, 90, "PLS-DA")

# legend
rect(xleft = 85, ybottom = 92, xright = 87, ytop = 94, col=coul[1])
text(90, 95, "Quantitative", cex=0.75)
rect(xleft = 85, ybottom = 87, xright = 87, ytop = 89, col=coul[2])
text(90, 90, "Qualitative", cex=0.75)


## Biological question

I am analysing a single data set (e.g. transcriptomics data) and I would like to classify my samples into known groups and predict the class of new samples. In addition, I am interested in identifying the key variables that drive such discrimination.

## The srbct study

The data are directly available in a processed and normalised format from the package. The Small Round Blue Cell Tumours (SRBCT) dataset from [@Kha01] includes the expression levels of 2,308 genes measured on 63 samples. The samples are classified into four classes as follows: 8 Burkitt Lymphoma (BL), 23 Ewing Sarcoma (EWS), 12 neuroblastoma (NB), and 20 rhabdomyosarcoma (RMS).

The srbct dataset contains the following:

$gene: a data frame with 63 rows and 2308 columns. The expression levels of 2,308 genes in 63 subjects.$class: a class vector containing the class tumour of each individual (4 classes in total).

$gene.name: a data frame with 2,308 rows and 2 columns containing further information on the genes. More details can be found in ?srbct. To illustrate PLS-DA, we will analyse the gene expression levels of srbct$gene to discriminate the 4 groups of tumours.

## Principle of sparse PLS-DA

Although Partial Least Squares was not originally designed for classification and discrimination problems, it has often been used for that purpose [@Ngu02a;@Tan04]. The response matrix Y is qualitative and is internally recoded as a dummy block matrix that records the membership of each observation, i.e. each of the response categories are coded via an indicator variable (see [@mixomics] Suppl. Information S1 for an illustration). The PLS regression (now PLS-DA) is then run as if Y was a continuous matrix. This PLS classification trick works well in practice, as demonstrated in many references [@Bar03;@Ngu02a;@Bou07;@Chung10].

Sparse PLS-DA [@Lec11] performs variable selection and classification in a one step procedure. sPLS-DA is a special case of sparse PLS described later in \@ref(pls), where $\ell_1$ penalization is applied on the loading vectors associated to the X data set.

## Inputs and outputs

We use the following data input matrices: X is a $n \times p$ data matrix, Y is a factor vector of length $n$ that indicates the class of each sample, and $Y^*$ is the associated dummy matrix ($n \times K$) with $n$ the number of samples (individuals), $p$ the number of variables and $K$ the number of classes. PLS-DA main outputs are:

• A set of components, also called latent variables. There are as many components as the chosen dimension of the PLS-DA model.

• A set of loading vectors, which are coefficients assigned to each variable to define each component. Those coefficients indicate the importance of each variable in PLS-DA. Importantly, each loading vector is associated to a particular component. Loading vectors are obtained so that the covariance between a linear combination of the variables from X (the X-component) and the factor of interest Y (the $Y^*$-component) is maximised.

• A list of selected variables from X and associated to each component if sPLS-DA is applied.

## Set up the data

We will mainly focus on sparse PLS-DA that is more suited for large biological data sets where the aim is to identify molecular signatures, as well as classifying samples. We first set up the data as X expression matrix and Y as a factor indicating the class membership of each sample. We also check that the dimensions are correct and match:

library(mixOmics)
data(srbct)
X <- srbct$gene Y <- srbct$class
summary(Y)
dim(X); length(Y)


## Quick start

For a quick start we arbitrarily set the number of variables to select to 50 on each of the 3 components of PLS-DA (see section \@ref(tuning:sPLSDA) for tuning these values).

MyResult.splsda <- splsda(X, Y, keepX = c(50,50)) # 1 Run the method
plotIndiv(MyResult.splsda)                          # 2 Plot the samples
plotVar(MyResult.splsda)                            # 3 Plot the variables
selectVar(MyResult.splsda, comp=1)$name # Selected variables on component 1  As PLS-DA is a supervised method, the sample plot automatically displays the group membership of each sample. We can observe a clear discrimination between the BL samples and the others on the first component (x-axis), and EWS vs the others on the second component (y-axis). Remember that this discrimination spanned by the first two PLS-DA components is obtained based on a subset of 100 variables (50 selected on each component). From the plotIndiv the axis labels indicate the amount of variation explained per component. Note that the interpretation of this amount is not the same as in PCA. In PLS-DA, the aim is to maximise the covariance between X and Y, not only the variance of X as it is the case in PCA! If you were to run splsda with this minimal code, you would be using the following default values: • ncomp = 2: the first two PLS components are calculated and are used for graphical outputs; • scale = TRUE: data are scaled (variance = 1, strongly advised here); It is important to note that all Discriminat Analyses (e.g. plsda, mint.splsda etc.) use mode='regression' for deflation. PLS-DA without variable selection can be performed as: MyResult.plsda <- plsda(X,Y) # 1 Run the method plotIndiv(MyResult.plsda) # 2 Plot the samples plotVar(MyResult.plsda) # 3 Plot the variables  ## To go further {#plsda-tgf} ### Customize sample plots {#splsda:plotIndiv} The sample plots can be improved in various ways. First, if the names of the samples are not meaningful at this stage, they can be replaced by symbols (ind.names=TRUE). Confidence ellipses can be plotted for each sample (ellipse = TRUE, confidence level set to 95\% by default, see the argument ellipse.level), Additionally, a star plot displays arrows from each group centroid towards each individual sample (star = TRUE). A 3D plot is also available, see plotIndiv for more details. plotIndiv(MyResult.splsda, ind.names = FALSE, legend=TRUE, ellipse = TRUE, star = TRUE, title = 'sPLS-DA on SRBCT', X.label = 'PLS-DA 1', Y.label = 'PLS-DA 2')  ### Customize variable plots The name of the variables can be set to FALSE (var.names=FALSE): plotVar(MyResult.splsda, var.names=FALSE)  In addition, if we had used the non sparse version of PLS-DA, a cut-off can be set to display only the variables that mostly contribute to the definition of each component. Those variables should be located towards the circle of radius 1, far from the centre. plotVar(MyResult.plsda, cutoff=0.7)  In this particular case, no variable selection was performed. Only the display was altered to show a subset of variables. ### Other useful plots #### Background prediction A 'prediction' background can be added to the sample plot by calculating a background surface first, before overlaying the sample plot. See ?background.predict for more details. More details about prediction, prediction distances can be found in [@mixomics] in the Suppl. Information. background <- background.predict(MyResult.splsda, comp.predicted=2, dist = "max.dist") plotIndiv(MyResult.splsda, comp = 1:2, group = srbct$class,
ind.names = FALSE, title = "Maximum distance",
legend = TRUE,  background = background)


#### ROC

As PLS-DA acts as a classifier, we can plot a ROC Curve to complement the sPLS-DA classification performance results detailed in \@ref(tuning:sPLSDA). The AUC is calculated from training cross-validation sets and averaged. Note however that ROC and AUC criteria may not be particularly insightful, or may not be in full agreement with the PLSDA performance, as the prediction threshold in PLS-DA is based on specified distance as described in [@mixomics].

auc.plsda <- auroc(MyResult.splsda)


### Variable selection outputs

First, note that the number of variables to select on each component does not need to be identical on each component, for example:

MyResult.splsda2 <- splsda(X,Y, ncomp=3, keepX=c(15,10,5))


Selected variables are listed in the selectVar function:

selectVar(MyResult.splsda2, comp=1)$value  and can be visualised in plotLoadings with the arguments contrib = 'max' that is going to assign to each variable bar the sample group colour for which the mean (method = 'mean') is maximum. See example(plotLoadings) for other options (e.g. min, median) plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean')  Interestingly from this plot, we can see that all selected variables on component 1 are highly expressed in the BL (orange) class. Setting contrib = 'min' would highlight that those variables are lowly expressed in the NB grey class, which makes sense when we look at the sample plot. Since 4 classes are being discriminated here, samples plots in 3d may help interpretation: plotIndiv(MyResult.splsda2, style="3d")  ### Tuning parameters and numerical outputs {#tuning:sPLSDA} For this set of methods, three parameters need to be chosen: 1 - The number of components to retain ncomp. The rule of thumb is usually$K - 1$where$K$is the number of classes, but it is worth testing a few extra components. 2 - The number of variables keepX to select on each component for sparse PLS-DA, 3 - The prediction distance to evaluate the classification and prediction performance of PLS-DA. For item 1, the perf evaluates the performance of PLS-DA for a large number of components, using repeated k-fold cross-validation. For example here we use 3-fold CV repeated 10 times (note that we advise to use at least 50 repeats, and choose the number of folds that are appropriate for the sample size of the data set): MyResult.plsda2 <- plsda(X,Y, ncomp=10) set.seed(30) # for reproducibility in this vignette, otherwise increase nrepeat MyPerf.plsda <- perf(MyResult.plsda2, validation = "Mfold", folds = 3, progressBar = FALSE, nrepeat = 10) # we suggest nrepeat = 50 # type attributes(MyPerf.plsda) to see the different outputs # slight bug in our output function currently see the quick fix below #plot(MyPerf.plsda, col = color.mixo(5:7), sd = TRUE, legend.position = "horizontal") # quick fix matplot(MyPerf.plsda$error.rate$BER, type = 'l', lty = 1, col = color.mixo(1:3), main = 'Balanced Error rate') legend('topright', c('max.dist', 'centroids.dist', 'mahalanobis.dist'), lty = 1, col = color.mixo(5:7))  The plot outputs the classification error rate, or Balanced classification error rate when the number of samples per group is unbalanced, the standard deviation according to three prediction distances. Here we can see that for the BER and the maximum distance, the best performance (i.e. low error rate) seems to be achieved for ncomp = 3. In addition (item 3 for PLS-DA), the numerical outputs listed here can be reported as performance measures: MyPerf.plsda  Regarding item 2, we now use tune.splsda to assess the optimal number of variables to select on each component. We first set up a grid of keepX values that will be assessed on each component, one component at a time. Similar to above we run 3-fold CV repeated 10 times with a maximum distance prediction defined as above. list.keepX <- c(5:10, seq(20, 100, 10)) list.keepX # to output the grid of values tested set.seed(30) # for reproducbility in this vignette, otherwise increase nrepeat tune.splsda.srbct <- tune.splsda(X, Y, ncomp = 3, # we suggest to push ncomp a bit more, e.g. 4 validation = 'Mfold', folds = 3, dist = 'max.dist', progressBar = FALSE, measure = "BER", test.keepX = list.keepX, nrepeat = 10) # we suggest nrepeat = 50  We can then extract the classification error rate averaged across all folds and repeats for each tested keepX value, the optimal number of components (see ?tune.splsda for more details), the optimal number of variables to select per component which is summarised in a plot where the diamond indicated the optimal keepX value: error <- tune.splsda.srbct$error.rate
ncomp <- tune.splsda.srbct$choice.ncomp$ncomp # optimal number of components based on t-tests on the error rate
ncomp
select.keepX <- tune.splsda.srbct$choice.keepX[1:ncomp] # optimal number of variables to select select.keepX plot(tune.splsda.srbct, col = color.jet(ncomp))  Based on those tuning results, we can run our final and tuned sPLS-DA model: MyResult.splsda.final <- splsda(X, Y, ncomp = ncomp, keepX = select.keepX) plotIndiv(MyResult.splsda.final, ind.names = FALSE, legend=TRUE, ellipse = TRUE, title="SPLS-DA, Final result")  Additionally we can run perf for the final performance of the sPLS-DA model. Also note that perf will output features that lists the frequency of selection of the variables across the different folds and different repeats. This is a useful output to assess the confidence of your final variable selection, see a more detailed example here. ## Additional resources Additional examples are provided in example(splsda) and in our case studies on our website in the Methods and Case studies sections, and in particular here. Also have a look at [@Lec11] ## FAQ • Can I discriminate more than two groups of samples (multiclass classification)? • Yes, this is one of the advantage of PLS-DA, see this example above • Can I have a hierarchy between two factors (e.g. diet nested into genotype)? • Unfortunately no, sparse PLS-DA only allows to discriminate all groups at once (i.e. 4 x 2 groups when there are 4 diets and 2 genotypes) • Can I have missing values in my data? • Yes in the X data set, but you won't be able to do any prediction (i.e. tune, perf, predict) • No in the Y factor # Projection to Latent Structure (PLS) {#pls} library(mixOmics) coul <- color.mixo(1:3) plot(0, type="n", xlim=c(0,100), ylim=c(83,97), axes=FALSE, xlab="",ylab="", main="PLS overview") box() # PLS rect(xleft = 20, ybottom = 85, xright = 50, ytop = 95, col=coul[1]) rect(xleft = 52, ybottom = 85, xright = 73, ytop = 95, col=coul[1]) text(5, 90, "PLS") # legend rect(xleft = 85, ybottom = 92, xright = 87, ytop = 94, col=coul[1]) text(90, 95, "Quantitative", cex=0.75)  ## Biological question I would like to integrate two data sets measured on the same samples by extracting correlated information, or by highlighing commonalities between data sets. ## The nutrimouse study The nutrimouse study contains the expression levels of genes potentially involved in nutritional problems and the concentrations of hepatic fatty acids for forty mice. The data sets come from a nutrigenomic study in the mouse from our collaborator [@Mar07], in which the effects of five regimens with contrasted fatty acid compositions on liver lipids and hepatic gene expression in mice were considered. Two sets of variables were measured on 40 mice: • gene: the expression levels of 120 genes measured in liver cells, selected among (among about 30,000) as potentially relevant in the context of the nutrition study. These expressions come from a nylon microarray with radioactive labelling. • lipid: concentration (in percentage) of 21 hepatic fatty acids measured by gas chromatography. • diet: a 5-level factor. Oils used for experimental diets preparation were corn and colza oils (50/50) for a reference diet (REF), hydrogenated coconut oil for a saturated fatty acid diet (COC), sunflower oil for an Omega6 fatty acid-rich diet (SUN), linseed oil for an Omega3-rich diet (LIN) and corn/colza/enriched fish oils for the FISH diet (43/43/14). • genotype 2-levels factor indicating either wild-type (WT) and PPAR$\alpha$-/- (PPAR). More details can be found in ?nutrimouse. To illustrate sparse PLS, we will integrate the gene expression levels (gene) with the concentrations of hepatic fatty acids (lipid). ## Principle of PLS Partial Least Squares (PLS) regression [@Wol66; @Wol01] is a multivariate methodology which relates (\textit{integrates}) two data matrices X (e.g. transcriptomics) and Y (e.g. lipids). PLS goes beyond traditional multiple regression by modelling the structure of both matrices. Unlike traditional multiple regression models, it is not limited to uncorrelated variables. One of the many advantages of PLS is that it can handle many noisy, collinear (correlated) and missing variables and can also simultaneously model several response variables in Y. PLS is a multivariate projection-based method that can address different types of integration problems. Its flexibility is the reason why it is the backbone of most methods in mixOmics. PLS is computationally very efficient when the number of variables$p + q >> n$the number of samples. It performs successive local regressions that avoid computational issues due to the inversion of large singular covariance matrices. Unlike PCA which maximizes the variance of components from a single data set, PLS maximizes the covariance between components from two data sets. The mathematical concepts of covariance and correlation are similar, but the covariance is an unbounded measure and covariance has a unit measure (see \@ref(intro:background)). In PLS, linear combination of variables are called latent variables or latent components. The weight vectors used to calculate the linear combinations are called the loading vectors. Latent variables and loading vectors are thus associated, and come in pairs from each of the two data sets being integrated. ## Principle of sparse PLS Even though PLS is highly efficient in a high dimensional context, the interpretability of PLS needed to be improved. sPLS has been recently developed by our team to perform simultaneous variable selection in both data sets X and Y data sets, by including LASSO$\ell_1$penalizations in PLS on each pair of loading vectors [@Lec08]. ## Inputs and outputs We consider the data input matrices: X is a$n \times p$data matrix and Y a$n \times q$data matrix, where$n$the number of samples (individuals),$p$and$q$are the number of variables in each data set. PLS main outputs are: • A set of components, also called latent variables associated to each data set. There are as many components as the chosen dimension of the PLS. • A set of loading vectors, which are coefficients assigned to each variable to define each component. Those coefficients indicate the importance of each variable in PLS. Importantly, each loading vector is associated to a particular component. Loading vectors are obtained so that the covariance between a linear combination of the variables from X (the X-component) and from Y (the$Y$-component) is maximised. • A list of selected variables from both X and Y and associated to each component if sPLS is applied. ## Set up the data We first set up the data as X expression matrix and Y as the lipid abundance matrix. We also check that the dimensions are correct and match: library(mixOmics) data(nutrimouse) X <- nutrimouse$gene
Y <- nutrimouse$lipid dim(X); dim(Y)  ## Quick start We will mainly focus on sparse PLS for large biological data sets where variable selection can help the interpretation of the results. See ?pls for a model with no variable selection. Here we arbitrarily set the number of variables to select to 50 on each of the 2 components of PLS (see section \@ref(tuning:PLS) for tuning these values). MyResult.spls <- spls(X,Y, keepX = c(25, 25), keepY = c(5,5)) plotIndiv(MyResult.spls) plotVar(MyResult.spls)  If you were to run spls with this minimal code, you would be using the following default values: • ncomp = 2: the first two PLS components are calculated and are used for graphical outputs; • scale = TRUE: data are scaled (variance = 1, strongly advised here); • mode = "regression": by default a PLS regression mode should be used (see \@ref(PLS:details) for more details) . Because PLS generates a pair of components, each associated to each data set, the function plotIndiv produces 2 plots that represent the same samples projected in either the space spanned by the X-components, or the Y-components. A single plot can also be displayed, see section \@ref(pls:plotIndiv). ## To go further {#pls-tgf} ### Customize sample plots {#pls:plotIndiv} Some of the sample plot additional arguments were described in \@ref(splsda:plotIndiv). In addition, you can choose the representation space to be either the components from the X-data set, the Y- data set, or an average between both components rep.space = 'XY-variate'. See more examples in examples(plotIndiv) and on our website. Here are two examples with colours indicating genotype or diet: plotIndiv(MyResult.spls, group = nutrimouse$genotype,
rep.space = "XY-variate", legend = TRUE,
legend.title = 'Genotype',
ind.names = nutrimouse$diet, title = 'Nutrimouse: sPLS')  plotIndiv(MyResult.spls, group=nutrimouse$diet,


#### Arrow plots

Instead of projecting the samples into the combined XY representation space, as shown in \@ref(pls:plotIndiv), we can overlap the X- and Y- representation plots. One arrow joins the same sample from the X- space to the Y- space. Short arrows indicate a good agreement found by the PLS between both data sets.

plotArrow(MyResult.spls,group=nutrimouse$diet, legend = TRUE, X.label = 'PLS comp 1', Y.label = 'PLS comp 2')  ### Variable selection outputs The selected variables can be extracted using the selectVar function for further analysis. MySelectedVariables <- selectVar(MyResult.spls, comp = 1) MySelectedVariables$X$name # Selected genes on component 1 MySelectedVariables$Y$name # Selected lipids on component 1  The loading plots help visualise the coefficients assigned to each selected variable on each component: plotLoadings(MyResult.spls, comp = 1, size.name = rel(0.5))  ### Tuning parameters and numerical outputs {#tuning:PLS} For PLS and sPLS, two types of parameters need to be chosen: 1 - The number of components to retain ncomp, 2 - The number of variables to select on each component and on each data set keepX and keepY for sparse PLS. For item 1 we use the perf function and repeated k-fold cross-validation to calculate the Q$^2$criterion used in the SIMCA-P software [@Ume96]. The rule of thumbs is that a PLS component should be included in the model if its value is$\leq 0.0975$. Here we use 3-fold CV repeated 10 times (note that we advise to use at least 50 repeats, and choose the number of folds that are appropriate for the sample size of the data set). We run a PLS model with a sufficient number of components first, then run perf on the object. MyResult.pls <- pls(X,Y, ncomp = 4) set.seed(30) # for reproducbility in this vignette, otherwise increase nrepeat perf.pls <- perf(MyResult.pls, validation = "Mfold", folds = 5, progressBar = FALSE, nrepeat = 10) plot(perf.pls$Q2.total)
abline(h = 0.0975)


This example seems to indicate that up to 3 components could be enough. In a small $p+q$ setting we generally observe a Q$^2$ that decreases, but that is not the case here as $n << p+q$.

Item 2 can be quite difficult to tune. Here is a minimal example where we only tune keepX based on the Mean Absolute Value. Other measures proposed are Mean Square Error, Bias and R2 (see ?tune.spls):

list.keepX <- c(2:10, 15, 20)
# tuning based on MAE
set.seed(30) # for reproducbility in this vignette, otherwise increase nrepeat
tune.spls.MAE <- tune.spls(X, Y, ncomp = 3,
test.keepX = list.keepX,
validation = "Mfold", folds = 5,
nrepeat = 10, progressBar = FALSE,
measure = 'MAE')
plot(tune.spls.MAE, legend.position = 'topright')


Based on the lowest MAE obtained on each component, the optimal number of variables to select in the X data set, including all variables in the Y data set would be:



### AUC

An AUC plot per block is plotted using the function auroc see [@mixomics] for the interpretation of such output as the ROC and AUC criteria are not particularly insightful in relation to the performance evaluation of our methods, but can complement the statistical analysis.

Here we evaluate the AUC for the model that includes 2 components in the miRNA data set.

Myauc.diablo <- auroc(MyResult.diablo, roc.block = "miRNA", roc.comp = 2)


#### Prediction on an external test set

The predict function predicts the class of samples from a test set. In our specific case, one data set is missing in the test set but the method can still be applied. Make sure the name of the blocks correspond exactly.

# prepare test set data: here one block (proteins) is missing
X.test <- list(mRNA = breast.TCGA$data.test$mrna,
miRNA = breast.TCGA$data.test$mirna)

Mypredict.diablo <- predict(MyResult.diablo, newdata = X.test)
# the warning message will inform us that one block is missing
#Mypredict.diablo # list the different outputs


The confusion table compares the real subtypes with the predicted subtypes for a 2 component model, for the distance of interest:

confusion.mat <- get.confusion_matrix(
truth = breast.TCGA$data.test$subtype,
predicted = Mypredict.diablo$MajorityVote$centroids.dist[,2])
kable(confusion.mat)
get.BER(confusion.mat)


### Tuning parameters

For DIABLO, the parameters to tune are:

1 - The design matrix design indicates which data sets, or blocks should be connected to maximise the covariance between components, and to which extend. A compromise needs to be achieved between maximising the correlation between data sets (design value between 0.5 and 1) and maximising the discrimination with the outcome Y (design value between 0 and 0.5), see [@Sin16] for more details.

2 - The number of components to retain ncomp. The rule of thumb is usually $K-1$ where $K$ is the number of classes, but it is worth testing a few extra components.

3 - The number of variables to select on each component and on each data set in the list keepX.

For item 1, by default all data sets are linked as follows:

MyResult.diablo$design  The design can be changed as follows. By default each data set will be linked to the Y outcome. MyDesign <- matrix(c(0, 0.1, 0.3, 0.1, 0, 0.9, 0.3, 0.9, 0), byrow=TRUE, ncol = length(X), nrow = length(X), dimnames = list(names(X), names(X))) MyDesign MyResult.diablo.design <- block.splsda(X, Y, keepX=list.keepX, design=MyDesign)  Items 2 and 3 can be tuned using repeated cross-validation, as we described in Chapter \@ref(plsda). A detailed tutorial is provided on our website in the different DIABLO tabs. ## Additional resources Additional examples are provided in example(block.splsda) and in our DIABLO tab in http://www.mixomics.org. Also have a look at [@Sin16] ## FAQ • When performing a multi-block analysis, how do I choose my design? • We recommend first relying on some prior biological knowledge you may have on the relationship you expect to see between data sets. Conduct a few trials on a non sparse version block.plsda, look at the classification performance with perf and plotDiablo before you can decide on your final design. • I have a small number of samples (n < 10), should I still tune keepX? • It is probably not worth it. Try with a few keepX values and look at the graphical outputs so see if they make sense. With a small$n$you can adopt an exploratory approach that does not require a performance assessment. • During tune or perf the code broke down (system computationally singular). • Check that the$M$value for your M-fold is not too high compared to$n$(you want$n/M > 6 - 8\$ as rule of thumb). Try leave-one-out instead with validation = 'loo' and make sure ncomp is not too large as you are running on empty matrices!
• My tuning step indicated the selection of only 1 miRNA...

• Choose a grid of keepX values that starts at a higher value (e.g. 5). The algorithm found an optimum with only one variable, either because it is highly discriminatory or because the data are noisy, but it does not stop you from trying for more.
• My Y is continuous, what can I do?

• You can perform a multi-omics regression with block.spls. We have not found a way yet to tune the results so you will need to adopt an exploratory approach or back yourself up with down stream analyses once you have identified a list of highly correlated features.

# Session Information

sessionInfo()
`

# References

## Try the mixOmics package in your browser

Any scripts or data that you put into this service are public.

mixOmics documentation built on Nov. 8, 2020, 11:12 p.m.