BiocStyle::markdown()
library(knitr) # global options knitr::opts_chunk$set(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_chunk$set( # collapse = TRUE, # comment = "#>" # )
This vignette is also available as a bookdown format on our gitHub page https://mixomicsteam.github.io/Bookdown/.
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.
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.
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.
Loadings: variable coefficients used to define a 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]).
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")
The methods implemented in mixOmics
are described in detail in the following publications. A more extensive list can be found at this link.
DIABLO (New): Singh A, Gautier B, Shannon C, Vacher M, Rohart F, Tebbutt S, K-A. Le Cao (2019). DIABLO: an integrative approach for identifying key molecular drivers from multi-omics assays. Bioinformatics bty1054.
Overview and recent integrative methods: Rohart F., Gautier, B, Singh, A, Le Cao, K. A. mixOmics: an R package for omics feature selection and multiple data integration. PLoS Comput Biol 13(11): e1005752.
Graphical outputs for integrative methods: [@Gon12] Gonzalez I., Le Cao K.-A., Davis, M.D. and Dejean S. (2012) Insightful graphical outputs to explore relationships between two omics data sets. BioData Mining 5:19.
sparse PLS: Le Cao K.-A., Martin P.G.P, Robert-Granie C. and Besse, P. (2009) Sparse Canonical Methods for Biological Data Integration: application to a cross-platform study. BMC Bioinformatics, 10:34.
sparse PLS-DA:Le Cao K.-A., Boitard S. and Besse P. (2011) Sparse PLS Discriminant Analysis: biologically relevant feature selection and graphical displays for multiclass problems. BMC Bioinformatics, 22:253.
sPLS-DA for microbiome data: Le Cao K-A$^$, Costello ME $^$, Lakis VA , Bartolo F, Chua XY, Brazeilles R and Rondeau P. (2016) MixMC: Multivariate insights into Microbial Communities.PLoS ONE 11(8): e0160169
Multilevel approach for repeated measurements: Liquet B, Le Cao K-A, Hocini H, Thiebaut R (2012). A novel approach for biomarker selection and the integration of repeated measures experiments from two assays. BMC Bioinformatics, 13:325
For those interested in how these methods have been applied to omics biological problems:
Lee AH, Shanon CS, [...], some members of mixOmics team and Kollman T (2019). Dynamic molecular changes during the first week of human life follow a robust developmental trajectory Nature Communications 10: 1092. We used DIABLO and other multi omics integrative methods to integratve 3 - 5 omics datasets.
Gavin PG, [...], and Hamilton-Williams EE (2018). Intestinal metaproteomics reveals host-microbiota interactions in subjects at risk for type 1 diabetes Diabetes care 41: 10. We used DIABLO to integrate microbiome, proteomics and meta-proteomics.
Each of the methods chapter has the following outline:
Other methods not covered in this document are described on our website and the following references:
regularised Canonical Correlation Analysis, see the Methods and Case study tabs, and [@Gon08] that describes CCA for large data sets.
Microbiome (16S, shotgun metagenomics) data analysis, see also [@Lec16] and kernel integration for microbiome data. The latter is in collaboration with Drs J Mariette and Nathalie Villa-Vialaneix (INRA Toulouse, France), an example is provided for the Tara ocean metagenomics and environmental data, see also [@Mar17].
MINT or P-integration to integrate independently generated transcriptomics data sets. An example in stem cells studies, see also [@Roh16].
First, download the latest \texttt{mixOmics} version from Bioconductor.:
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") # then load 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 data <- read.csv("your_data.csv", row.names = 1, header = TRUE) # from txt file data <- read.table("your_data.txt", header = TRUE) # 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!
mixOmics
{#start:PCA}Each analysis should follow this workflow:
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) X <- nutrimouse$gene
Then use the following steps:
MyResult.pca <- pca(X) # 1 Run the method plotIndiv(MyResult.pca) # 2 Plot the samples plotVar(MyResult.pca) # 3 Plot the variables
This is only a first quick-start, there will be many avenues you can take to deepen your exploratory and integrative analyses. The package proposes several methods to perform variable, or feature selection to identify the relevant information from rather large omics data sets. The sparse methods are listed in the Table in \@ref(intro:overview).
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.
Do not stop here! We are not done yet. You can enhance your analyses with the following:
Have a look at our manual and each of the functions and their examples, e.g. ?pca
, ?plotIndiv
, ?sPCA
, ...
Run the examples from the help file using the example
function: example(pca)
, example(plotIndiv)
, ...
Have a look at out website that features many tutorials and case studies,
Keep reading this vignette, this is just the beginning!
library(mixOmics) coul <- color.mixo(1:3) plot(0, type="n", xlim=c(0,100), ylim=c(83,97), axes=FALSE, xlab="",ylab="", main="PCA overview") box() # PCA rect(xleft = 20, ybottom = 85, xright = 60, ytop = 95, col=coul[1]) text(5, 90, "PCA") # legend rect(xleft = 85, ybottom = 90, xright = 87, ytop = 92, col=coul[1]) text(90, 94, "Quantitative", cex=0.75) #rect(xleft = 75, ybottom = 90, xright = 77, ytop = 92, col=coul[2]) #text(90, 91, "Qualitative", cex=0.75)
I would like to identify the major sources of variation in my data and itdentify whether such sources of variation correspond to biological conditions, or experimental bias. I would like to visualise trends or patterns between samples,whether they 'naturally' cluster according to known biological conditions.
liver.toxicity
studyThe liver.toxicity
is a list in the package that contains:
gene
: a data frame with 64 rows and 3116 columns, corresponding to the expression levels of 3,116 genes measured on 64 rats.
clinic
: a data frame with 64 rows and 10 columns, corresponding to the measurements of 10 clinical variables on the same 64 rats.
treatment
: data frame with 64 rows and 4 columns, indicating the treatment information of the 64 rats, such as doses of acetaminophen and times of necropsy.
gene.ID
: a data frame with 3116 rows and 2 columns, indicating geneBank IDs of the annotated genes.
More details are available at ?liver.toxicity
.
To illustrate PCA, we focus on the expression levels of the genes in the data frame liver.toxicity$gene
. Some of the terms mentioned below are listed in \@ref(intro:background).
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.
We first load the data from the package. See \@ref(start:upload) to upload your own data.
library(mixOmics) data(liver.toxicity) X <- liver.toxicity$gene
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.
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.
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.
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')
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 plotLoadings(MyResult.spca, comp = 2)
For this set of methods, two parameters need to be chosen:
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.
Additional reading in [@She08].
Should I scale my data before running PCA? (scale = TRUE
in pca
)
Can I perform PCA with missing values?
When should I apply a multilevel approach in PCA? (multilevel
argument in PCA
)
When should I apply a CLR transformation in PCA? (logratio = 'CLR'
argument in PCA
)
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)
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.
srbct
studyThe 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.
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.
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.
We first load the data from the package. See \@ref(start:upload) to upload your own 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)
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
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')
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.
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)
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)
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")
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 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]
Can I discriminate more than two groups of samples (multiclass classification)?
Can I have a hierarchy between two factors (e.g. diet nested into genotype)?
Can I have missing values in my data?
tune, perf, predict
)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)
I would like to integrate two data sets measured on the same samples by extracting correlated information, or by highlighing commonalities between data sets.
nutrimouse
studyThe 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
).
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.
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].
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.
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)
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).
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, pch = nutrimouse$genotype, rep.space = "XY-variate", legend = TRUE, legend.title = 'Diet', legend.title.pch = 'Genotype', ind.names = FALSE, title = 'Nutrimouse: sPLS')
See (example(plotVar)
) for more examples. Here we change the size of the labels. By default the colours are assigned to each type of variable. The coordinates of the variables can also be saved as follows:
plotVar(MyResult.spls, cex=c(3,2), legend = TRUE) coordinates <- plotVar(MyResult.spls, plot = FALSE)
We extended other types of plots, based on clustered image maps and relevance networks to ease the interpretation of the relationships between two types of variables. A similarity matrix is calculated from the outputs of PLS and represented with those graphics, see [@Gon12] for more details, and our website
A clustered image map can be produced using the cim
function. You may experience figures margin issues in RStudio. Best is to either use X11()
or save the plot as an external file. For example to show the correlation structure between the X and Y variables selected on component 1:
X11() cim(MyResult.spls, comp = 1) cim(MyResult.spls, comp = 1, save = 'jpeg', name.save = 'PLScim')
Using the same similarity matrix input in CIM, we can also represent relevance bipartite networks. Those networks only represent edges between on type of variable from X
and the other type of variable, from Y
. Whilst we use sPLS to narrow down to a few key correlated variables, our keepX
and keepY
values might still be very high for this kind of output. A cut-off can be set based on the correlation coefficient between the different types of variables.
Other arguments such as interactive = TRUE
enables a scrollbar to change the cut-off value interactively, see other options in ?network
. Additionally, the graph object can be saved to be input into Cytoscape for an improved visualisation.
X11() network(MyResult.spls, comp = 1) network(MyResult.spls, comp = 1, cutoff = 0.6, save = 'jpeg', name.save = 'PLSnetwork') # save as graph object for cytoscape myNetwork <- network(MyResult.spls, comp = 1)$gR
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')
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))
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:
tune.spls.MAE$choice.keepX
Tuning keepX
and keepY
conjointly is still work in progress. What we advise in the meantime is either to adopt an arbitrary approach by setting those parameters arbitrarily, depending on the biological question, or tuning one parameter then the other.
You may have noticed the mode
argument in PLS. We can calculate the residual matrices at each PLS iteration differently. Note: this is for advanced users.
Regression mode
The PLS regression mode models uni-directional (or 'causal') relationship between two data sets. The Y
matrix is deflated with respect to the information extracted/modelled from the local regression on X. Here the goal is to predict Y
from X
(Y
and X
play an assymmetric role). Consequently the latent variables computed to predict Y from X are different from those computed to predict X from Y. More details about the model can be found in[ Appendix [@Lec08].
PLS regression mode, also called PLS2, is commonly applied for the analysis of biological data [@Bou05; @Byl07] due to the biological assumptions or the biological dogma. In general, the number of variables in Y
to predict are fewer in number than the predictors in X
.
Canonical mode
Similar to a Canonical Correlation Analysis (CCA) framework, this mode is used to model a bi-directional (or symmetrical) relationship between the two data sets. The Y
matrix is deflated with respect to the information extracted or modelled from the local regression on Y
. Here X
and Y
play a symmetric role and the goal is similar to CCA. More details about the model can be found in [@Lec09a].
PLS canonical mode is not well known (yet), but is applicable when there is no a priori relationship between the two data sets, or in place of CCA but when variable selection is required in large data sets. In [@Lec09a], we compared the measures of the same biological samples on different types of microarrays, cDNA and Affymetrix arrays, to highlight complementary information at the transcripts levels. Note however that for this mode we do not provide any tuning function.
Other modes
The 'invariant' mode performs a redundancy analysis, where the Y
matrix is not deflated. The 'classic' mode is similar to a regression mode. It gives identical results for the variates and loadings associated to the X
data set, but differences for the loadings vectors associated to the Y data set (different normalisations are used). Classic mode is the PLS2 model as defined by [@Ten98], Chap 9.
Difference between PLS modes For the first PLS dimension, all PLS modes will output the same results in terms of latent variables and loading vectors. After the first dimension, these vectors will differ, as the matrices are deflated differently.
Additional examples are provided in example(spls)
and in our case studies on our website in the Methods and Case studies sections, see also [@Lec08; @Lec09a].
Can PLS handle missing values?
perf
or tune
is not possible with missing values.Can PLS deal with more than 2 data sets?
DIABLO
(Chapter \@ref(diablo)) for multi-block analysesWhat are the differences between sPLS and Canonical Correlation Analysis (CCA, see ?rcca
in mixOmics)?
Can I perform PLS with more variables than observations?
Can I perform PLS with 2 data sets that are highly unbalanced (thousands of variables in one data set and less than 10 in the other ?
keepX
.library(mixOmics) coul <- color.mixo(1:3) plot(0, type="n", xlim=c(0,100), ylim=c(83,97), axes=FALSE, xlab="",ylab="", main="DIABLO overview") box() rect(xleft = 12, ybottom = 85, xright = 30, ytop = 95, col=coul[1]) rect(xleft = 32, ybottom = 85, xright = 45, ytop = 95, col=coul[1]) points(x=47, y=90, pch=16, col='black', cex=0.5) points(x=48.5, y=90, pch=16, col='black', cex=0.5) points(x=50, y=90, pch=16, col='black', cex=0.5) rect(xleft = 52, ybottom = 85, xright = 61, ytop = 95, col=coul[1]) rect(xleft = 63, ybottom = 85, xright = 65, ytop = 95, col=coul[2]) text(5, 90, "DIABLO") # 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)
DIABLO is our new framework in mixOmics
that extends PLS for multiple data sets integration and PLS-discriminant analysis. The acronyms stands for Data Integration Analysis for Biomarker discovery using a Latent cOmponents [@Sin16].
I would like to identify a highly correlated multi-omics signature discriminating known groups of samples.
breast.TCGA
studyHuman 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 levels of mRNA expression [@Sor01. Here we consider a subset of data generated by The Cancer Genome Atlas Network [@TCGA12]. For the package, data were normalised and drastically prefiltered for illustrative purposes but DIABLO can handle larger data sets, see [@mixomics] Table 2. 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 (proteomics 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
is a list containing training and test sets of omics data data.train
and data.test
which both 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 only for the protein abundance.
subtype
: a factor indicating the breast cancer subtypes in the training (length of 150) and test (length of 70) sets.
More details can be found in ?breast.TCGA
.
To illustrate DIABLO, 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 test samples in the test set.
The core DIABLO method extends Generalised Canonical Correlation Analysis [@Ten11], which contrary to what its name suggests, generalises PLS for multiple matching datasets, and the sparse sGCCA method [@Ten14]. Starting from the R package RGCCA, we extended these methods for different types of analyses, including unsupervised N-integration (block.pls
, block.spls
) and supervised analyses (block.plsda
, block.splsda
).
The aim of N-integration with our sparse methods is to identify correlated (or co-expressed) variables measured on heterogeneous data sets which also explain the categorical outcome of interest (supervised analysis). The multiple data integration task is not trivial, as the analysis can be strongly affected by the variation between manufacturers or omics technological platforms despite being measured on the same biological samples. Before you embark on data integration, we strongly suggest individual or paired analyses with sPLS-DA and PLS to first understand the major sources of variation in each data set and to guide the integration process.
More methodological details can be found in [@Sin16].
We consider as input a list X
of data frames with $n$ rows (the number of samples) and different number of variations in each data frame. Y
is a factor vector of length $n$ that indicates the class of each sample. Internally and similar to PLS-DA in Chapter \@ref(plsda) it will be coded as a dummy matrix.
DIABLO 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 DIABLO.
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 DIABLO. 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 each data set and associated to each component if sparse DIABLO is applied.
We first set up the input data as a list of data frames X
expression matrix and Y
as a factor indicating the class membership of each sample. Each data frame in X
should be named consistently to match with the keepX
parameter.
We check that the dimensions are correct and match. We then set up arbitrarily the number of variables keepX
that we wish to select in each data set and each component.
library(mixOmics) data(breast.TCGA) # extract training data and name each data frame X <- list(mRNA = breast.TCGA$data.train$mrna, miRNA = breast.TCGA$data.train$mirna, protein = breast.TCGA$data.train$protein) Y <- breast.TCGA$data.train$subtype summary(Y) list.keepX <- list(mRNA = c(16, 17), miRNA = c(18,5), protein = c(5, 5))
MyResult.diablo <- block.splsda(X, Y, keepX=list.keepX) plotIndiv(MyResult.diablo) plotVar(MyResult.diablo)
Similar to PLS (Chapter \@ref(pls)), DIABLO generates a pair of components, each associated to each data set. This is why we can visualise here 3 sample plots. As DIABLO is a supervised method, samples are represented with different colours depending on their known class.
The variable plot suggest some correlation structure between proteins, mRNA and miRNA. We will further customize these plots in Sections \@ref(diablo:plotIndiv) and \@ref(diablo:plotVar).
If you were to run block.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 for data integration);We focused here on the sparse version as would like to identify a minimal multi-omics signature, however, the non sparse version could also be run with block.plsda
:
MyResult.diablo2 <- block.plsda(X, Y)
Here is an example of an improved plot, see also Section \@ref(splsda:plotIndiv) for additional sources of inspiration.
plotIndiv(MyResult.diablo, ind.names = FALSE, legend=TRUE, cex=c(1,2,3), title = 'BRCA with DIABLO')
Labels can be omitted in some data sets to improve readability. For example her we only show the name of the proteins:
plotVar(MyResult.diablo, var.names = c(FALSE, FALSE, TRUE), legend=TRUE, pch=c(16,16,1))
Several plots were added for the DIABLO framework.
plotDiablo
A global overview of the correlation structure at the component level can be represented with the plotDiablo
function. It plots the components across the different data sets for a given dimension. Colours indicate the class of each sample.
plotDiablo(MyResult.diablo, ncomp = 1)
Here, we can see that a strong correlation is extracted by DIABLO between the mRNA and protein data sets. Other dimensions can be plotted with the argument comp
.
circosPlot
The circos plot represents the correlations between variables of different types, represented on the side quadrants. Several display options are possible, to show within and between connexions between blocks, expression levels of each variable according to each class (argument line = TRUE
). The circos plot is built based on a similarity matrix, which was extended to the case of multiple data sets from [@Gon12]. A cutoff
argument can be included to visualise correlation coefficients above this threshold in the multi-omics signature.
circosPlot(MyResult.diablo, cutoff=0.7)
cimDiablo
The cimDiablo
function is a clustered image map specifically implemented to represent the multi-omics molecular signature expression for each sample. It is very similar to a classic hierarchical clustering:
# minimal example with margins improved: # cimDiablo(MyResult.diablo, margin=c(8,20)) # extended example: cimDiablo(MyResult.diablo, color.blocks = c('darkorchid', 'brown1', 'lightgreen'), comp = 1, margin=c(8,20), legend.position = "right")
plotLoadings
The plotLoadings
function visualises the loading weights of each selected variables on each component (default is comp = 1
) and each data set. The color indicates the class in which the variable has the maximum level of expression (contrib = "max"
) or minimum (contrib ="min"
), on average (method="mean"
) or using the median (method ="median"
). We only show the last plot here:
#plotLoadings(MyResult.diablo, contrib = "max") plotLoadings(MyResult.diablo, comp = 2, contrib = "max")
Another visualisation of the correlation between the different types of variables is the relevance network, which is also built on the similarity matrix [@Gon12]. Each color represents a type of variable. A threshold can also be set using the argument cutoff
.
See also Section \@ref(pls:network) to save the graph and the different options, or ?network
.
network(MyResult.diablo, blocks = c(1,2,3), color.node = c('darkorchid', 'brown1', 'lightgreen'), cutoff = 0.6, save = 'jpeg', name.save = 'DIABLOnetwork')
Similar to what is described in Section \@ref(tuning:sPLSDA) we use repeated cross-validation with perf
to assess the prediction of the model. For this complex classification problems, often a centroid distance is suitable, see details in [@mixomics] Suppl. Material S1.
set.seed(123) # for reproducibility in this vignette MyPerf.diablo <- perf(MyResult.diablo, validation = 'Mfold', folds = 5, nrepeat = 10, dist = 'centroids.dist') #MyPerf.diablo # lists the different outputs # Performance with Majority vote #MyPerf.diablo$MajorityVote.error.rate
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)
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)
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 examples are provided in example(block.splsda)
and in our DIABLO tab in http://www.mixomics.org. Also have a look at [@Sin16]
When performing a multi-block analysis, how do I choose my design?
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
?
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
).
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...
My Y is continuous, what can I do?
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.sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.