Note: this case study is quite detailed for illustrative purposes. You can skip some parts for your analyses. For more details about the methods and performance measures we use, have a look at our list of references, and any of the other tabs focusing on specific methods.
Here we apply mixMC to the 16S study of Koren et al. (2011)  who examined the link between oral, gut and plaque microbial communities in patients with atherosclerosis and controls. Only healthy individuals were retained in the analysis. This study contains measures from different bodysites, saliva and stool and arterial plaque. The design of the experiment is partially repeated (15 unique patients for saliva and stool and 13 for arterial plaque). For simplicity here we consider the experiment as non-repeated. See the HMP example in the Tab for a repeated measures analysis.
The data were downloaded from the Qiita database. After pre-filtering, the data included 980 OTU for 43 samples. Data were TSS normalized and log ratio transformations are applied in PCA and PLS-DA.
Using the Koren data in this example, we assume the data are offset and pre-filtered, as described in mixMC pre-processing steps.
library(mixOmics) data("Koren.16S") #ls(Koren.16S) # offset data, filtered data = Koren.16S$data.raw
We set up the outcome Y as a factor indicating the body sites of each sample.
Y = Koren.16S$bodysite summary(Y)
Unsupervised Analysis with PCA
PCA allows for dimension reduction of the data and visualisation of diversity patterns in microbiome studies. Because we are dealing with compositional data, spurious results can arise if the data are not transformed into an Euclidean space. Therefore, we use CLR transformation , using the argument
logratio = 'CLR' in our functions.
We first run a PCA with a sufficiently large number of components
ncomp to choose the final reduced dimension of the model.
koren.pca = pca(data, ncomp = 10, logratio = 'CLR') plot(koren.pca)
The barplot above depicts the percentage of explained variance per component for a PCA with 10 components. It is a good indicator of how many components to retain in the final PCA model (here two components will be sufficient).
The sample plot below displays the samples coloured according to Y, but note that PCA does not take Y into account in the analysis!
plotIndiv(koren.pca, comp = c(1,2), # the components to plot ind.names = F, # not showing sample names group = Y, # color according to Y legend = TRUE, title = 'Koren, PCA comp 1 - 2')
Supervised Analysis and Selection of Discriminative OTUs with sPLS-DA
The mixMC framework uses the sPLS-DA multivariate analysis from mixOmics . The sPLS-DA selects the most discriminative OTUs that best characterise each body site, see sPLS-DA.
Note that with PLS-DA and sPLS-DA we can only choose a CLR transformation (see details here.
We first run the
perf function with a PLS-DA model with no variable selection. Often
ncomp = K-1 where K is the number of categories in the outcome Y is sufficient, but it depends on the data. We choose 5-fold cross-validation repeated 10 times. To obtain a more reliable estimation of the error rate, the number of repeats should be increased (between 50 to 100).
Here we assess the performance of the PLS-DA on 3 components.
koren.plsda = plsda(X = data, Y, ncomp = nlevels(Y), logratio = 'CLR') koren.perf.plsda = perf(koren.plsda, validation = 'Mfold', folds = 5, progressBar = FALSE, nrepeat = 10) plot(koren.perf.plsda, overlay = 'measure', sd=TRUE)
The plot indicates a decrease in the classification error rate (i.e. an increase in classification performance) from one component to 3 components in the model. The BER stands for Balanced Error Rate and should be considered when there is an unbalanced number of samples per group. Here the number of samples per group is very similar, which explains why both overall and BER overlap. The performance reaches its best for 3 components, which suggests
ncomp = 3 for a final PLS-DA model. Note that for the sparse PLSDA we may obtain a different
For more details on the prediction distances, refer to Supplementary information in Rohart et al. 2017 . Here all distances are similar, we choose 'max.dist' for sPLS-DA.
#First two components plotIndiv(koren.plsda , comp = c(1,2), group = Y, ind.names = FALSE, ellipse = TRUE, legend = TRUE, title = 'Koren, PLSDA comp 1 - 2')
It is possible that the third component may not visually add much to the model and could be ignored.
The parameters to choose in sPLS-DA are the number of variables to select (
keepX) and the number of components (
ncomp). To do this we use the function
tune.splsda. In this example the sPLS-DA tuning step is performed on 3 components, using 5-fold validation repeated 10 times.
tune.splsda needs to be performed prior to the sPLS-DA analysis to choose the parameters on a grid of
keepX values (make sure you choose the appropriate M fold cross-validation and provide sufficient
nrepeat in the evaluation model, except for loo where it can only be run on 1 repeat). Also check the stability of the features selected during the cross-validation process.
# first, set a grid of values to test: grid.keepX = c(seq(5,150, 5)) # if you dont understand what this means, type # grid.keepX # adjust this grid as necessary for your own data # this chunk takes ~2 min to run set.seed(33) # for reproducible results for this code, remove for your own code koren.tune.splsda = tune.splsda(data, Y = Y, ncomp = 3, logratio = 'CLR', test.keepX = grid.keepX, validation = c('Mfold'), folds = 5, dist = 'max.dist', # prediction distance can be chosen according to tune.plsda results nrepeat = 10, progressBar = FALSE) # may show some convergence issues for some of the cases, it is ok for tuning
The graphic above shows that the error rate decreases when 2 components are included in the sPLS-DA, whereas the third component seems to add noise. The diamonds indicate the optimal
keepX variables to select on each component based on the balanced error rate.
# optimal number of variables to select on 2 comps: select.keepX = koren.tune.splsda$choice.keepX[1:2] select.keepX
## comp1 comp2 ## 130 40
We now run a classic sPLS-DA.
Note: if you have not tuned your sPLS-DA model (or you are unhappy with the size of the selection) you can set the optimal
keepX manually, e.g.:
# select.keepX = c(50, 65) # manually choose size of selection
# the sPLS-DA koren.splsda = splsda(X = data, Y = Y, ncomp = 2, keepX = select.keepX, logratio= "CLR")
The sample plot below shows the first two components for the sPLS-DA. The ellipse are 0.95 confidence interval ellipse for each body site.
plotIndiv(koren.splsda, comp = c(1,2), ind.names = F, ellipse = TRUE, legend = TRUE, title = 'Koren, sPLS-DA comp 1-2')
The classification performance of the sPLS-DA multilevel model can be assessed using the function
perf. The mean error rates per component and the type of distance are output. It can be beneficial to increase the number of repeats for more accurate estimations.
set.seed(34) # for reproducible results for this code, remove for your own code koren.perf.splsda = perf(koren.splsda, validation = 'Mfold', folds = 5, progressBar = FALSE, nrepeat = 10, dist = 'max.dist')
The mean overall and Balanced Error Rate are output for each component.
## $overall ## max.dist ## comp1 0.3046512 ## comp2 0.0000000 ## ## $BER ## max.dist ## comp1 0.3297436 ## comp2 0.0000000
OTU selection and plots
The sPLS-DA selects the most discriminative OTUs that best characterize each body site. The contribution plots below display the abundance of each OTU (large abundance = large absolute value) and in which body site they are the most abundant for each sPLS-DA component. The contribution plots need to be interpreted in combination with the sample plot below to understand the similarities between body sites, in addition to answer the question “Which bacteria characterise those body sites?”
The code below outputs the first selected OTUs and their coefficient (from the loading vector) on the first component. Consider the absolute value as an indication of the importance of the OTU in the microbial signature, while the sign generally indicates positive or negative correlations between the OTUs, relative to the proportions of the others (see also with a
For example on comp 2:
head(selectVar(koren.splsda, comp = 2)$value)
## value.var ## 406145 -0.3725344 ## 370772 -0.3191646 ## 1089121 0.2676551 ## 2356875 -0.2651746 ## 388506 -0.2557496 ## 925707 0.2444380
To this list of selected OTUs displayed from the most important to the least important, we can combine their stability measure from the
perf output (i.e. how often were they selected across the different CV runs).
selected.OTU.comp2 = selectVar(koren.splsda, comp = 2)$name # stability of OTUs selected on comp 2 koren.perf.splsda$features$stable[][selected.OTU.comp2]
plotLoadings plot shows that all OTUs selected on the first component are mostly highly abundant in faecal samples (based on their mean per body site). Let's display only the top 50 selected OTUs:
plotLoadings(koren.splsda, comp = 1, method = 'mean', contrib = 'max', size.title = 1, ndisplay = 50, size.name = 0.5, size.legend = 0.3)
plotLoadings plot shows that all OTUs selected on the second component are highly abundant in the two other body sites. The sign indicates the opposition between body sites.
plotLoadings(koren.splsda, comp = 2, method = 'mean', contrib = 'max', size.title = 1, size.name = 0.5, size.legend = 0.3)
We could improve the plot by showing only the top OTUs, and with their names at the Family level, followed by their OTU ID:
name.var = paste(Koren.16S$taxonomy[, 'Family'], colnames(data), sep = '|') plotLoadings(koren.splsda, comp = 2, method = 'mean', contrib = 'max', name.var = name.var, size.title = 1, size.name = 0.5, size.legend = 0.3)
Clustered Image Map
A heatmap will also help to understand the microbial signature. We represent Clustered Image Maps (with Euclidian distance, Ward linkage set by default) for the OTUs selected on each sPLS-DA component. The abundance values that are displayed are the normalised, log ratio transformed values. All OTUs selected by the sPLS-DA model are displayed, other options can include a specific component, or a specific cutoff of 'association', see
Note, if the margins are too large, use either
X11() prior, or use the
name.save in the
#CIM plot cim(koren.splsda, row.sideColors = color.mixo(Y))
Koren, O., Knights, D., Gonzalez, A., Waldron, L., Segata, N., Knight, R., Huttenhower, C. and Ley, R.E., 2013. A guide to enterotypes across the human body: meta-analysis of microbial community structures in human microbiome datasets. PLoS Comput Biol, 9(1), p.e1002863