Introduction to multivariate analysis with mixOmics

Let’s get started

Installation

First, download the latest version from Bioconductor (Rohart et al. 2017):

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

The mixOmics package should directly import the following packages: igraph, rgl, ellipse, corpcor, RColorBrewer, plyr, parallel, dplyr, tidyr, reshape2, methods, matrixStats, rARPACK, gridExtra.

For Apple mac users: if you are unable to install the imported package rgl, you will need to install the XQuartz software first.

Load the package

library(mixOmics)

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

Upload data

The examples in this workshop use data that are already part of the package. 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 or via one of these command lines:

# from csv file
data <- read.csv("your_data.csv", row.names = 1, header = TRUE)

# from txt file
data <- read.table("your_data.txt", header = TRUE)

For more details about the arguments used to modify those functions, type ?read.csv or ?read.table in the R console.

Quick start in mixOmics

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!

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

data(nutrimouse)       #load data
X <- nutrimouse$gene    #store data in object called X

Then use the following steps:

MyResult.pca <- pca(X)  # 1 Run the method PCA on X
plotIndiv(MyResult.pca) # 2 Plot the samples from a PCA result

plotVar(MyResult.pca)   # 3 Plot the variables from a PCA result

This is only a first quick-start. The package proposes several methods to perform variable selection and integration.

Following our example here, sparse PCA can be applied to select the top 5 variables contributing to each of the two components in PCA. The user specifies the number of variables to selected on each component, for example, here 5 variables are selected on each of the first two components (keepX=c(5,5)):

MyResult.spca <- spca(X, keepX=c(5,5)) # 1 Run the method
plotIndiv(MyResult.spca)               # 2 Plot the samples

plotVar(MyResult.spca)                 # 3 Plot the variables

You can see know that we have considerably reduced the number of genes in the plotVar correlation circle plot.

1. PCA on the SRBCT case study

The Small Round Blue Cell Tumours (SRBCT) data set from (Khan et al. 2001) includes the expression levels of 2,308 genes measured on 63 samples. The samples are divided into four classes: 8 Burkitt Lymphoma (BL), 23 Ewing Sarcoma (EWS), 12 neuroblastoma (NB), and 20 rhabdomyosarcoma (RMS). The data are directly available in a processed and normalised format from the mixOmics package and contains the following:

• $gene: A data frame with 63 rows and 2,308 columns. These are the expression levels of 2,308 genes in 63 subjects,

• $class: A vector containing the class of tumour for 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.

Load the data

We first load the data from the package and store the gene expression data in a \(\boldsymbol X\) object.

library(mixOmics)
data(srbct)
X <- srbct$gene
dim(X)  # check dimension

Exploration with PCA

PCA is a useful tool to explore the gene expression data and to assess sample similarities between tumour types (information available in srbct$class). Remember that PCA is an unsupervised approach, but we can colour the samples by their tumour subtype to assist in interpreting the PCA. Here we center (default argument) and scale the data:

pca.srbct <- pca(X, ncomp = 10, scale = TRUE)

pca.srbct

plot(pca.srbct)

# simple version
plotIndiv(pca.srbct, 
          pch = 1,    # pch is to show symbols
          title = 'SRBCT, PCA comp 1 - 2')

# advanced version
plotIndiv(pca.srbct, 
          group = srbct$class, # asking to color according to class of tumour
          ind.names = FALSE,   # not showing the sample names
          legend = TRUE, 
          title = 'SRBCT, PCA comp 1 - 2')

Exercise 1: Interpretation

Run the code above, line by line.

1. How many principal components would you retain, and why.

2. How is a sample plot obtained in PCA? What does it represent?

3. Interpret the sample plot outputs (the simple version first, then the advanced version, as shown in the code). What is the major source of variation in the data?

Show solutions

Answers

1. The first two components explain 21.96% of the total variance. Based on the elbow method, we would choose 2 components.

2. A sample plot is obtained by plotting the components as a way to reduce the dimension of the data. It means that samples are projected into the space spanned by the principal components 1 and 2.

3. The tumour types are not clustered, meaning that the major source of variation cannot be explained by tumour types. Instead, samples seem to cluster according to an unknown source of variation. We observe almost no separation between the different tumour types in the PCA sample plot, with perhaps the exception of the NB samples that tend to cluster with other samples. This preliminary exploration teaches us two important findings:

• The major source of variation is not attributable to tumour type, but an unknown source (we tend to observe clusters of samples but those are not explained by tumour type).

• We need a more ‘directed’ (supervised) analysis to separate the tumour types, and we should expect that the amount of variance explained by the dimensions in PLS-DA analysis will be small since they will be focusing on discriminating the sample groups.

2. PLS-DA on the SRBCT case study

Set the outcome variable

We have already loaded the SRBCT gene expression data and stored in \(\boldsymbol X\). We need to store the factor indicating the sample class membership in \(\boldsymbol Y\)

Y <- srbct$class 
length(Y)

First pass: PLS-DA

We run a PLS-DA model that includes three components:

plsda.srbct <- plsda(X,Y, ncomp = 3)

plotIndiv(plsda.srbct, ind.names = FALSE, legend=TRUE,
          comp=c(1,2), 
          title = 'PLS-DA on SRBCT comp 1-2')

plotLoadings(plsda.srbct, 
             contrib = 'max',  # max = which class has maximum mean (median) expression (vs. min)
             method = 'mean',  # choose either mean or median
             comp = 1,
             ndisplay = 50)

plotVar(plsda.srbct, cutoff = 0.7)  # variable plot (correlation circle plot)

biplot(plsda.srbct, cutoff = 0.7)   # we set a correlation cutoff to only show the most contributing variables   

boxplot(X[, 'g1932']~ Y)

Demo and interpretation

Let’s run all the command lines above and discuss about the following points:

1. Interpretation of the sample plot. What do we notice regarding the amount of explained variance per component, compared to PCA. Is this expected? Let’s also produce the sample plot for component 1 vs. 3.

2. Loadings plot: what does each bar represent? Let’s change the parameters contrib = 'max' to contrib = 'min'. Let’s now change the parameter comp = 1 to comp = 2. What do we observe? Now, let’s interpret the loading plot in combination with the sample plot. Does the loading plot make sense in light of this additional information?

3. Interpretation of the variable plot (correlation circle plot, refer to slides).

4. Interpretation of the biplot.

5. As a follow-up to the biplot interpretation, let’s focus on one gene and plot a boxplot with respect to the tumour subtype group (as shown in code). What is the link between how the gene is located on the biplot, and its expression levels per sample group?

Show solutions

Interpretation

1.1 Amount of explained variance: while PCA maximises the variance of the components, PLSDA maximises the covariance (akin to correlation) between the data and the outcome, hence the amount of explained variance can be lower than PCA (especially in this case, as we have seen already that the major source of variation in the data is not attributable to tumour subtypes).

1.2 Sample plot: By default, the samples are coloured according to their class membership. Component 1 discriminates RMS + EWS vs. NB + BL, component 2 discriminates RMS + NB vs. EWS + BL, while component 3 discriminates further the NB and BL groups. It is the combination of all three components that enables us to discriminate all classes.

2.1 Loading plots: each bar represent the contribution (= coefficient) of each gene in each component. Here we only show the top 50 genes with the highest contribution (positive or negative sign, going from the bottom of the graph to the top). The argument contrib = 'max' colors the bars according to the tumour subtype for which the mean expression value is maximal. The argument contrib = 'min' colors the bars according to the tumour subtype for which the mean expression value is minimal. When we change to the argument comp = 2, we observe a new set of genes and the color patterns are different.

2.2 Loading plots and sample plot: from the sample plot we have observed that component 1 discriminates RMS + EWS vs. NB + BL. We observe the same trend at the gene level from the loading plot, showing that these genes are the one driving the separation between this group of samples. (similar interpretation for component 2). Hence, the loading plot confirms what we have observed at the sample level on the sample plot.

3. Variable plot: here we project the variables in the space spanned by the first two components by calculating the correlation between the variables, and the components (refer to slides). We only show the top genes with a correlation cutoff above 0.7 across both components. We observe a groups of genes that are positively correlated together (left and right hand side on the x-axis), and contribute to component 1 either negatively (left) or positively (right). These two groups are negatively correlated to each other. The trend is not as clear for component 2 (less genes and contribute less to component 2 on the y-axis).

4. Biplot: overlays both samples (as dots) and variables (as arrows). Here we only show the top genes with a correlation cutoff above 0.7 across both components. We can now visualise the same set of genes as in the correlation circle plot, and how they relate (are over / under expressed) with the samples.

5. Genes that ‘points’ towards a specific group of samples tend to be overexpressed in these samples.

Sparse PLS-DA for variable selection

We run similar command lines, but this time we use sparse PLS-DA that will enable us to select the top discriminant variables.

ⓘ Sparse means zeroes in the loading vectors

Sparse stands for loading vectors that are sparse (include zeroes) because of we use lasso regularisation in the method (Lê Cao, Boitard, and Besse 2011) so that only a subset of variables are used to calculate the components in the linear combination.

The variable selection process is done within PLS-DA. For this, we need to specify the number of variables to select per component.

splsda.srbct <- splsda(X,Y, ncomp = 3,
                       keepX = c(20, 10, 10) # select 20 genes on the first component, 10 on the second, and 10 on the third component
                       )

plotIndiv(splsda.srbct, ind.names = FALSE, legend=TRUE,
          comp=c(1,2), 
          title = 'PLS-DA on SRBCT comp 1-2')

plotLoadings(splsda.srbct, 
             contrib = 'max',  # max = which class has maximum mean (median) expression (vs. min)
             method = 'mean',  # choose either mean or median
             comp = 1)

plotVar(splsda.srbct)  # variable plot (correlation circle plot)

selectVar(splsda.srbct) #extracts the list of variables selected, and their loading weights

Exercise 2 on sPLS-DA: Interpretation

Run the code above, line by line.

1. What are the differences that you observe across all outputs, compared to the full PLS-DA run previously (i.e with no variable selection).

2. Plot a sample plot for component 1 vs 3, what do you observe? Also plot the loadings for component 3 for more insights.

3. Change the values of the keepX argument, do you notice large changes?

4. Conclude on the advantages / disadvantages of using sparse PLS-DA for variable selection.

Sorry! we currently have a bug with the biplot that we need to fix!

Show solutions

Answers

1. On the sample plot, the NB and RMS are clustered together, whereas the other groups are well separated from each other. Component 1 discriminates BL vs. the other classes. Component 2 discriminates mostly BL vs EWS vs the others. The loading plot is now much clearer (only one color shown on each side of the positive / negative sign). The loading plot only shows the variables that were selected on each component. Similarly for the correlation circle plot, genes are now better aligned on either the x- or y- axis with stronger contribution.

2. Sample plot of comp 1 vs comp 3 shows a clearer separation of NB vs. RMS compared to Comp 1 vs Comp 2. This is also reflected in the loading plot on component 3.

3. If we stick to values that are not too large (i.e up to 50 or so), then the plots and trends remain the same. When the keepX values get very large, then we tend to obtain the same results as PLS-DA.

4. Conclusion: sPLS-DA can help identifying key discriminative genes for each sample group. It does so by adopting a ‘one group’ vs ‘all’ (or 1 group vs another group) per component. This enables us to identify which genes discriminate which class. The output results are usually neater than a classic PLS-DA. The main challenge is to decide how many genes we should select per component. We do this using repeated cross-validation and estimating the classification error. This is outside the scope of this workshop, but you can have a look at our examples in http://mixomics.org/case-studies/splsda-srbct-case-study/. There are also functions in mixOmics to estimate the classification performance of the final sPLS-DA.

3. DIABLO to integrate multi-omics datasets

DIABLO is a a method to integrate multiple data sets while explaining their relationship with a categorical outcome variable. DIABLO stands for Data Integration Analysis for Biomarker discovery using Latent variable approaches for Omics studies (Singh et al. 2019). It can also be referred to as Multiblock (s)PLS-DA.

ⓘ We need the same samples across all omics

One pre-requisite of DIABLO is that the same samples / individuals should be measured across all omics. This is because we are looking for correlation / covariance between the data sets by calculating the covariance between components (all of length \(N\), the number of samples).

TCGA case study

Human breast cancer is a heterogeneous disease in terms of molecular alterations, cellular composition, and clinical outcome. Breast tumours can be classified into several subtypes, according to their levels of mRNA expression (Sørlie et al. 2001). Here we consider a subset of data generated by The Cancer Genome Atlas Network (Cancer Genome Atlas Network et al. 2012). For the package, data were normalised, and then drastically prefiltered for illustrative purposes.

The data were divided into a training set with a subset of 150 samples from the mRNA, miRNA and proteomics data, and a test set including 70 samples, but only with mRNA and miRNA data (the proteomics data are missing). The aim of this integrative analysis is to identify a highly correlated multi-omics signature discriminating the breast cancer subtypes Basal, Her2 and LumA.

The breast.TCGA (more details can be found in ?breast.TCGA) is a list containing training and test sets of omics data data.train and data.test which include:

• $miRNA: A data frame with 150 (70) rows and 184 columns in the training (test) data set for the miRNA expression levels,
• $mRNA: A data frame with 150 (70) rows and 520 columns in the training (test) data set for the mRNA expression levels,
• $protein: A data frame with 150 rows and 142 columns in the training data set for the protein abundance (there are no proteomics in the test set),
• $subtype: A factor indicating the breast cancer subtypes in the training (for 150 samples) and test sets (for 70 samples).

This case study covers an interesting scenario where one omic data set is missing in the test set, but because the method generates a set of components per training data set, we can still assess the prediction or performance evaluation using majority or weighted prediction vote.

Load the data

We will integrate the expression levels of miRNA, mRNA and the abundance of proteins while discriminating the subtypes of breast cancer, then predict the subtypes of the samples in the test set.

Each omics data matrix is stored into a list of matrices \(\boldsymbol X\). Each data frame in \(\boldsymbol X\) should be named as we will match these names with the keepX parameter for variable selection. A factor indicating the class membership of each sample is stored in \(\boldsymbol Y\).

data(breast.TCGA)

# Extract training data and name each data frame
# Store as list
X <- list(mRNA = breast.TCGA$data.train$mrna, 
          miRNA = breast.TCGA$data.train$mirna, 
          protein = breast.TCGA$data.train$protein)

# Outcome
Y <- breast.TCGA$data.train$subtype
summary(Y)

## Basal  Her2  LumA 
##    45    30    75

Set the design matrix

In DIABLO, we need to specify a design matrix that indicates what is the level of correlation we wish to extract between datasets. In the example below, we ask for a small amount of correlation between these data sets.

design <- matrix(0.1, ncol = length(X), nrow = length(X), 
                dimnames = list(names(X), names(X)))
diag(design) <- 0
design 

##         mRNA miRNA protein
## mRNA     0.0   0.1     0.1
## miRNA    0.1   0.0     0.1
## protein  0.1   0.1     0.0

More details on how to choose the design, either based on apriori knowledge / expectation, or on data-driven approach:

• https://mixomics.org/methods/n-integration-methods/
• https://mixomics.org/mixdiablo/diablo-tcga-case-study/ (initial analysis part)

DIABLO

Here we decide to run a model with 2 components. The number of variables to keep is listed below in keepX (optimal number according to repeated cross-validation. More details available on our website).

# number of variables to select per dataset and per component (here 2 comps)
list.keepX <- list( mRNA = c(8, 25), miRNA = c(14,5), protein = c(10, 5))

diablo.tcga <- block.splsda(X, Y, ncomp = 2, 
                            keepX = list.keepX, design = design)

## Design matrix has changed to include Y; each block will be
##             linked to Y.

# the message tells us that each data set will be linked to the outcome Y for maximal discrimination

# sample plots:
plot(diablo.tcga) # pairs of components across datasets and their correlation

plotIndiv(diablo.tcga, ind.names = FALSE, legend = TRUE)

# variable plots:
plotVar(diablo.tcga, legend = TRUE)

circosPlot(diablo.tcga, , cutoff = 0.7)

plotLoadings(diablo.tcga, contrib = 'max',  # max = which class has maximum mean (median) expression (vs. min)
             method = 'mean',  # choose either mean or median
             comp = 1)

# variable selection
selectVar(diablo.tcga)

Exercise 3: Interpretation

Run the code above, line by line.

1. Sample plots: What do the sample plots plot and plotIndiv tell you about the ability of DIABLO to extract correlated information between data sets and to discriminate sample groups? Are there datasets that are more noisy than others? (remember that we decompose each data matrix with a set of components, and loading vectors, while maximising the correlation between the datasets, see slides).

2. Variable plots: What do the variable plots plotVar and circosPlot tell you about the correlation between the variables that were selected by DIABLO?

3. Loading plots: inspect the loading plots for components 1 and 2. What are the discriminative properties of the variables selected on component 1? and component 2? (the function selectVar outputs the list of variables selected, and their loading weights)

Show solutions

Answers

1.1 The plot output tells us that the mRNA and protein data sets are the most correlated. The least correlated pair is miRNA and protein. On component 1, DIABLO is able to extract correlated information between the data sets, and we observe some separation between the sample groups.

1.2 The plotIndiv output plots the set of components per data set. The mRNA and the protein datasets have some good discriminative properties. The miRNA dataset is a bit more noisy.

2.1 The variable plot plotVar shows a high positive correlation between miRNAs and a few proteins and one mRNA (x-axis, right hand side) or between mRNA and proteins (x-axis, left hand side). These groups are negatively correlated together, as they are either overexpressed in Basal or LumA (confirmed in the loading plots).

2.2 The variable plot circosPlot shows mostly a negative correlation between some of the variables, with a cutoff of 0.7.

3 According to the loading plot, the variables selected by DIABLO on component 1 discriminate Basal (positive loading value) vs. LumA (negative loading value), whereas the variables selected by DIABLO on component 2 discriminate Her2 (positive loading value) vs. the two other groups (negative loading value).

Prediction on the test set (optional)

The predict function associated predicts the class of samples from an external test set. In our specific case, one data set (proteomics) is missing in the test set but the method can still be applied for sample prediction. We need to ensure the names of the blocks correspond exactly to those from the training set:

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

predict.diablo.tcga <- predict(diablo.tcga, newdata = data.test.tcga)
# The warning message will inform us that one block is missing

predict.diablo.tcga # Lists the different outputs

The following output is a confusion matrix that compares the real subtypes with the predicted subtypes from our trained DIABLO model on the second component for the prediction distance centroids.dist and the prediction scheme WeightedVote (each dataset casts a vote):

confusion.mat.tcga <- get.confusion_matrix(truth = breast.TCGA$data.test$subtype, 
                     predicted = predict.diablo.tcga$WeightedVote$centroids.dist[,2]) # 2nd component
confusion.mat.tcga

##       predicted.as.Basal predicted.as.Her2 predicted.as.LumA
## Basal                 20                 1                 0
## Her2                   0                13                 1
## LumA                   0                 3                32

From this table, we see that one Basal and one Her2 sample are wrongly predicted as Her2 and Lum A respectively, and 3 LumA samples are wrongly predicted as Her2. The balanced prediction error rate can be obtained as:

get.BER(confusion.mat.tcga)

## [1] 0.06825397

It would be worthwhile at this stage to revisit the chosen design of DIABLO to assess the influence of the design on the prediction performance on this test set - even though this back and forth analysis is a biased criterion to choose the design!

To go further

We have not discussed how to choose the number of components, or the number of variables keepX to select on each dataset. This requires repeated cross-validation. More details on the tuning, and the classification performance of the method can be found here: http://mixomics.org/mixdiablo/diablo-tcga-case-study/

Further details, analyses and vignette are avaiable on www.mixOmics.org

Session Info

sessionInfo()

## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.5
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Australia/Melbourne
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] mixOmics_6.32.0 ggplot2_4.0.2   lattice_0.22-7  MASS_7.3-65    
## [5] knitr_1.50     
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10         generics_0.1.4      tidyr_1.3.2        
##  [4] stringi_1.8.7       digest_0.6.37       magrittr_2.0.4     
##  [7] evaluate_1.0.4      grid_4.5.2          RColorBrewer_1.1-3 
## [10] fastmap_1.2.0       plyr_1.8.9          jsonlite_2.0.0     
## [13] Matrix_1.7-4        ggrepel_0.9.6       RSpectra_0.16-2    
## [16] gridExtra_2.3       purrr_1.2.1         scales_1.4.0       
## [19] codetools_0.2-20    jquerylib_0.1.4     cli_3.6.6          
## [22] rlang_1.1.7         withr_3.0.2         cachem_1.1.0       
## [25] yaml_2.3.10         ellipse_0.5.0       tools_4.5.2        
## [28] parallel_4.5.2      reshape2_1.4.5      BiocParallel_1.42.2
## [31] dplyr_1.2.1         corpcor_1.6.10      vctrs_0.7.1        
## [34] R6_2.6.1            matrixStats_1.5.0   lifecycle_1.0.5    
## [37] stringr_1.6.0       pkgconfig_2.0.3     pillar_1.11.0      
## [40] bslib_0.10.0        gtable_0.3.6        glue_1.8.0         
## [43] rARPACK_0.11-0      Rcpp_1.1.0          xfun_0.53          
## [46] tibble_3.3.0        tidyselect_1.2.1    rstudioapi_0.18.0  
## [49] farver_2.1.2        htmltools_0.5.9     igraph_2.1.4       
## [52] labeling_0.4.3      rmarkdown_2.29      compiler_4.5.2     
## [55] S7_0.2.1

Cancer Genome Atlas Network et al. 2012. “Comprehensive Molecular Portraits of Human Breast Tumours.” Nature 490 (7418): 61–70.

Khan, Javed, Jun S Wei, Markus Ringner, Lao H Saal, Marc Ladanyi, Frank Westermann, Frank Berthold, et al. 2001. “Classification and Diagnostic Prediction of Cancers Using Gene Expression Profiling and Artificial Neural Networks.” Nature Medicine 7 (6): 673–79.

Lê Cao, Kim-Anh, Simon Boitard, and Philippe Besse. 2011. “Sparse PLS Discriminant Analysis: Biologically Relevant Feature Selection and Graphical Displays for Multiclass Problems.” BMC Bioinformatics 12 (1): 253.

Rohart, Florian, Benoit Gautier, Amrit Singh, and Kim-Anh Lê Cao. 2017. “mixOmics: An r Package for ‘Omics Feature Selection and Multiple Data Integration.” PLoS Computational Biology 13 (11): e1005752.

Singh, Amrit, Casey P Shannon, Benoı̂t Gautier, Florian Rohart, Michaël Vacher, Scott J Tebbutt, and Kim-Anh Lê Cao. 2019. “DIABLO: An Integrative Approach for Identifying Key Molecular Drivers from Multi-Omics Assays.” Bioinformatics 35 (17): 3055–62.

Sørlie, Therese, Charles M Perou, Robert Tibshirani, Turid Aas, Stephanie Geisler, Hilde Johnsen, Trevor Hastie, et al. 2001. “Gene Expression Patterns of Breast Carcinomas Distinguish Tumor Subclasses with Clinical Implications.” Proceedings of the National Academy of Sciences 98 (19): 10869–74.