Multiomics workshop part 3: Multi-omics integration with DIABLO
Quentin D. Read, Matt Burnham, & Matthias Anagho-Mattanovich
July 07, 2026
Background
Up to now, we have been working with a single omics dataset at once. Although we can do multivariate analyses with a single omics dataset, we have greater ambitions. We want to know whether there are patterns in the data that cut across multiple omics and also are associated with our variable of interest (time, in this case study).
So far, in Lessons 1 and 2, we’ve worked with the transcriptome and proteome data, pre-processing them and exploring them individually. There is also a metabolome dataset, but we are not going to go through that dataset in a separate lesson. We are going to go ahead and do multi-omics analysis using DIABLO on all three datasets together. Here’s what we’ll be doing in this lesson (we’ll explain the different steps in more detail as we go):
- Combine the three pre-transformed omics datasets together and do any final pre-processing
- Do partial least squares analysis on each pair of omics to get correlation parameters
- Fit an initial DIABLO model
- Cross-validate the initial model to determine the optimal number of components
- Cross-validate the model again to determine the optimal number of features per component
- Fit the final DIABLO model
- Make a bunch of plots to make sure the model is valid and explore patterns in the data that the model reveals
- One last cross-validation to assess the final model performance
N-Integration
N-Integration is the framework of having multiple datasets which measure different aspects of the same samples. For example, you may have transcriptomic, proteomic, methylation, and microbiome data for the same set of cells. N-integrative methods are built to use the information in all of these dataframes simultaneously.
DIABLO is a novel mixOmics framework for the integration
of 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. It can also be
referred to as Multiblock (s)PLS-DA.
What is sPLS-DA?
DIABLO is based on a technique called sPLS-DA, which stands for sparse partial least squares discriminant analysis. But what is that? Let’s start by explaining what partial least squares is and (vaguely) how it works.
Partial least squares
Partial least squares is a statistical technique. Like many of the statistical techniques we are familiar with, starting with fitting a line through points in a scatter plot and going all the way up to the neural networks underlying AI large language models, it is a form of regression. The objective of PLS is to find the relationships between two matrices X and Y, by identifying the multidimensional direction in X that explains the maximum variance in the Y multidimensional space. Partial least squares discriminan analysis (PLS-DA) is just a variant of PLS where Y is a single categorical variable. We can use PLS-DA in this case because our response variable is just time point, with three discrete categories. The goal of PLS-DA is to summarize the original data in fewer dimensions, maximizing variation between groups and minimizing variation within groups.
In PLS, we typically choose a number of components, meaning the number of dimensions we want to reduce our many-dimensional dataset down to, while preserving as much information about it as possible. You can see in the diagram below that our original dataset X is of dimensions \(N \times P\), where N is the number of samples and P the number of features. Usually in omics analyses, P is much bigger than N: in other words we have a lot more variables (genes, metabolites, microbiome taxa, etc., which could number up to the thousands) than we have independent samples (maybe in the 10-100 range). The goal of PLS is to preserve as much information as possible from the original \(N \times P\) dataset using k components (k = 3 in the diagram). We estimate Q components, which are each N-length vectors, and k loading vectors, which are each of length P. The product of a \(N \times k\) matrix times a \(k \times P\) matrix is a matrix of dimensions \(N \times P\). We estimate the values of the component and loading vectors to get as close as possible to the original data. If P is a really big number, and k is small, and our PLS model does a good job of recapturing the initial data, we have reduced the amount of data points needed to explain the pattern from the original \(NP\) down to \(Nk + kP = Q(N+P)\) which is tiny in comparison.
Sparse partial least squares discriminant analysis
Sparse partial least squares (sPLS), which includes a special variant for categorical responses called sPLS-DA, is a form of PLS where only a subset of the variables in X are used to explain Y. This is useful both for explaining variation in the data and for predicting the response. It helps with explanation because it narrows down the pool of candidate variables, and it helps with prediction because it protects against overfitting.
sPLS-DA uses an algorithm called LASSO to reduce the number of variables from the dataset that are used to construct the loading vectors. In other words, the loading vectors that we estimated in PLS are now forced to have a lot of zeros in them. The amount of information we are using to summarize the information in the original dataset is even less.
Multiblock sPLS-DA
Multiblock sPLS-DA, the technique used by DIABLO, extends sPLS-DA to more than one X dataset. Now we’re trying to find a set of components that does two things: maximizes the correlation between our datasets \(X_1, X_2, ... X_Q\) and, as before, explains the most variation in our categorical variable of interest Y.
A NOTE OF CAUTION: It’s important to note that the analyses we are doing here are biologically naive. We know that there are complex biological processes going on within these cells that involve RNA transcripts, proteins, and various metabolites. We’re not explicitly incorporating any knowledge about those processes into the models we are creating here. We are just mining the data for statistical associations between features across multiple multivariate datasets. I use the word “features” intentionally here; the analytical techniques we are using here are completely ignorant of whether those vectors of numbers represent RNA, protein, or anything else. To the model, they are just numbers. Thus, it is very important to keep in mind that this type of analysis is exploratory in nature. It can help draw attention to patterns existing in the data, that may indicate that certain biological processes are going on. We can use the patterns to speculate about the processes. Of course, these data do come from a manipulative experiment where cells were randomly assigned to treatments; so that does give us some basis to talk about causation when it comes to the treatment groups. But when it comes to the statistical associations between features in the different omics datasets, ultimately it is an exploratory/correlational analysis. We can use it to draw a lot of pretty pictures that are interesting to talk about, but it should be viewed as a jumping-off point for more detailed future studies into the specific biological processes behind them. Rant over!!!
Setup
Load packages
Load packages; as before, you may need to install packages from CRAN or Bioconductor if you are running this code on your own machine.
library(mixOmics)
library(Mus.musculus)
custom_colors <- c("salmon", "darkolivegreen3", "turquoise3")Import pre-transformed data
These data have already been transformed to the \(log_2\) scale, and features with low or no variance have been filtered out. This was done in the first two notebooks. In addition, we are loading some model tuning outputs into our workspace so that we don’t have to wait for the code to run later.
transcriptomics_log2_z <- read.csv('data/transcriptomics_log2.csv', row.names = 1)
proteomics_log2_z <- read.csv('data/proteomics_log2.csv', row.names = 1)
metabolomics_log2_z <- read.csv('data/metabolomics_filtered_log2.csv', row.names = 1)
load('data/diablo_tune_multi.rda')If you are running this code on your own laptop and not on the remote server, you will need to load the data from GitHub Pages:
transcriptomics_log2_z <- read.csv('https://usda-ree-ars.github.io/SEAStatsData/omics2026/transcriptomics_log2.csv', row.names = 1)
proteomics_log2_z <- read.csv('https://usda-ree-ars.github.io/SEAStatsData/omics2026/proteomics_log2.csv', row.names = 1)
metabolomics_log2_z <- read.csv('https://usda-ree-ars.github.io/SEAStatsData/omics2026/metabolomics_filtered_log2.csv', row.names = 1)
load(file = url('https://usda-ree-ars.github.io/SEAStatsData/omics2026/diablo_tune_multi.rda'))Prepare data for integration
Put the three data frames into a list. The data were imported in a
format where each row is a feature and each column is a sample, but
mixOmics requires the opposite: samples in rows and features in columns.
So, we use the t() function to transpose each of the data
frames. Also, mixOmics doesn’t play well with missing data, so we use
the na.omit() function to get rid of any rows with even a
single missing value. That’s a very simplistic way of dealing with
missing data; in practice you may want to play with multiple imputation
or other methods for imputing missing data. In this case, we end up
removing only a few hundred of the ~15,000 transcriptome features, and
none of the proteome features, but we lose about half of the 171
metabolome features.
data <- list(transcripts = t(na.omit(transcriptomics_log2_z)),
metabolites = t(na.omit(metabolomics_log2_z)),
proteins = t(na.omit(proteomics_log2_z)))Check the dimensions of each of the three omics datasets. They all have 18 rows, as expected, one for each sample (6 each from 3 different time points), with differing numbers of columns.
lapply(data, dim)## $transcripts
## [1] 18 14222
##
## $metabolites
## [1] 18 94
##
## $proteins
## [1] 18 3321Check to make sure each of the three omics datasets have their rows in the same order.
lapply(data, row.names)## $transcripts
## [1] "X0h_01" "X0h_02" "X0h_03" "X0h_04" "X0h_05" "X0h_06" "X4h_01"
## [8] "X4h_02" "X4h_03" "X4h_04" "X4h_05" "X4h_06" "X24h_01" "X24h_02"
## [15] "X24h_03" "X24h_04" "X24h_05" "X24h_06"
##
## $metabolites
## [1] "X0h_01" "X0h_02" "X0h_03" "X0h_04" "X0h_05" "X0h_06" "X4h_01"
## [8] "X4h_02" "X4h_03" "X4h_04" "X4h_05" "X4h_06" "X24h_01" "X24h_02"
## [15] "X24h_03" "X24h_04" "X24h_05" "X24h_06"
##
## $proteins
## [1] "X0h_01" "X0h_02" "X0h_03" "X0h_04" "X0h_05" "X0h_06" "X4h_01"
## [8] "X4h_02" "X4h_03" "X4h_04" "X4h_05" "X4h_06" "X24h_01" "X24h_02"
## [15] "X24h_03" "X24h_04" "X24h_05" "X24h_06"The samples are all ordered in the same way in each dataset, in
increasing order by time point. Create a vector Y to assign
the samples to groups (time points in this case).
Y <- factor(rep(c('0h', '4h', '24h'), each = 6), levels = c('0h', '4h', '24h'))One final setup step. As we did in part 1, let’s use the Mus.musculus package, which contains mouse gene annotation information, to replace the numerical gene IDs with annotations where possible. Some tweaking of the code was needed to not have duplicate column names, which can cause problems later on.
genes <- AnnotationDbi::select(
Mus.musculus,
keys = colnames(data$transcripts),
columns = c("SYMBOL", "TXCHROM", "ENTREZID"),
keytype = "ENSEMBL"
)
genes <- genes[!duplicated(genes$ENSEMBL), ]
colnames(data$transcripts)[!duplicated(genes$SYMBOL) & !is.na(genes$SYMBOL)] <- genes$SYMBOL[!duplicated(genes$SYMBOL) & !is.na(genes$SYMBOL)]Pairwise correlations between omics datasets
Before we do the full multiomics analysis, we need to look at the correlations between each pair of omics datasets. This will help us determine the design weights to use in the DIABLO analysis. What are those? Well, we have to assign a weight between 0 and 1 to each pair of omics datasets. The weight we choose depends on the goal of the analysis. Higher weights close to 1 prioritize a descriptive analysis that explains the maximum amount of correlation in the dataset. Lower weights give higher prediction accuracy but do not have as much explanatory power. If prediction is your goal, then you should use a small weight, for example 0.1.
Instead of trying to guess what weight we should use, let’s use a data-driven approach. We will do a quick partial least squares (PLS) analysis, with only one component, on each pair of omics datasets. Then calculate the correlation between the components associated with each of the datasets from the PLS.
pls1 <- pls(data[["transcripts"]], data[["metabolites"]], ncomp = 1)
pls2 <- pls(data[["transcripts"]], data[["proteins"]], ncomp = 1)
pls3 <- pls(data[["metabolites"]], data[["proteins"]], ncomp = 1)
# calculate correlation of transcripts and metabolites
print(cor(pls1$variates$X, pls1$variates$Y))## comp1
## comp1 0.9829229# calculate correlation of transcripts and proteins
print(cor(pls2$variates$X, pls2$variates$Y))## comp1
## comp1 0.9495981# calculate correlation of metabolites and proteins
print(cor(pls3$variates$X, pls3$variates$Y))## comp1
## comp1 0.9241838The pairwise correlations are all very high, at least 0.9. We can use 0.9 as the design weight in our DIABLO analysis.
Tuning the DIABLO model
Create the design matrix with a diagonal of 0s and all off-diagonal entries being 0.9.
design <- matrix(0.9, ncol = length(data), nrow = length(data),
dimnames = list(names(data), names(data)))
diag(design) <- 0
design## transcripts metabolites proteins
## transcripts 0.0 0.9 0.9
## metabolites 0.9 0.0 0.9
## proteins 0.9 0.9 0.0Fit an initial DIABLO model using the block.splsda()
function. We input our list of three omics datasets called
data to the X argument of the function, and
our categorical time point variable as Y. We specify
ncomp = 4 to fit up to four components, and we also pass
the design matrix we created earlier.
basic.diablo.model <- block.splsda(X = data, Y = Y, ncomp = 4, design = design)The next step is to tune the model using cross-validation to
determine how many components is optimal to explain the most variation
between groups in Y without overfitting. Now, because we
have only three groups in Y, it is unlikely that we will
need more than two components — if you have two groups in your data, you
should only need one component to separate them, and you’d need two
components for three groups, three for four groups, etc. So, we could
probably use common sense to say that two components are adequate
without doing the cross-validation. But let’s do it anyway to confirm
that intuition.
The perf() function is used to tune the model. The
arguments validation = 'Mfold', folds = 3 tell it to do
3-fold cross-validation. What this does is divide the data up into three
equal groups and fit the model three times, each time excluding the data
from one group. Then, use the model that was fit to the other two groups
to predict the classes of the observations in the excluded group. Repeat
this process for 1-4 components, then repeat the entire thing 10 times
(nrepeat = 10) to make sure we didn’t get a fluke result
from the way we randomly split up the data.
perf.diablo <- perf(basic.diablo.model, validation = 'Mfold',
folds = 3, nrepeat = 10, seed = 1234)The plot() method for the cross-validation output shows
two types of error rate: overall (ER) and balanced (BER). Because we
have exactly 6 samples at each time point, the data are perfectly
balanced, so the ER and BER lines are identical here. If we had
observational data that was imbalanced, they might differ. It also shows
three types of distance metric that can be used for the prediction.
What does the plot show us? First of all, it shows that the error
rate declines as the number of components increases. This makes sense
because we are allowing the model to have more parameters, the more
components are included. It also shows that the error rate is basically
zero when we have \(k = 2\) components.
As we guessed beforehand, more than 2 is not necessary and will really
only lead to overfitting. We also see that the metric
max.dist is not very good when there’s only one component,
but the difference between distance metrics is negligible at 2
components and up.
plot(perf.diablo) 
The perf.diablo object also tells us the optimal number
of components it identified, using various methods. Here is one of them.
Two of the three distance metrics identify 2 as the optimal number of
components. Mahalanobis distance identifies 4 as optimal, but it
performed less well than the other metrics at 2 components, so we’ll
stick with 2. Simpler is better.
perf.diablo$choice.ncomp$WeightedVote## max.dist centroids.dist mahalanobis.dist
## Overall.ER 2 2 4
## Overall.BER 2 2 4Given the results of the first tuning, let’s choose \(k = 2\) components and use the
centroids.dist metric because it always performed as well
or better than the other metrics.
We will now do a second tuning step, this time to select a number of
features to keep on each component. Remember that DIABLO is based on
sparse PLS, meaning that only a subset of the original data is
used. We are also using cross-validation for this part of the process,
but this time we have to fit the model many more times because of all
the different combinations of numbers of features we are trying out. We
do the cross-validation for all possible combinations out of the list in
test.keepX. Then repeat the whole process 10 times.
This takes quite a bit of time to run. Because it is a relatively
small omics dataset, I did not bother to do anything to speed up the
process. If you are working with a bigger dataset, consider looking into
the BiocParallel package which will allow you to
parallelize the computations so that multiple processors can work
simultaneously, each fitting models with different combinations of
numbers of features.
For teaching purposes, there is no need for us to sit here and wait for the cross-validation to run. Familiarize yourself with the code, but don’t bother to execute it. We already loaded the pre-tuned model output into our workspace earlier, when we imported the data.
test.keepX <- list(
transcripts = c(1, 2, 3, 4, 5, 10, 15, 20, 30, 40, 50),
metabolites = c(1, 2, 3, 4, 5, 10, 15, 20, 30, 40, 50),
proteins = c(1, 2, 3, 4, 5, 10, 15, 20, 30, 40, 50)
)
tune.MULTI <- tune.block.splsda(X = data, Y = Y, ncomp = 2,
test.keepX = test.keepX, design = design,
validation = 'Mfold', folds = 3, nrepeat = 10,
dist = "centroids.dist", seed = 123)How many features on each component gave the best predictive performance in cross-validation? We see that other than gene transcripts, very few variables are needed to capture most of the information in the original dataset. In particular, component 2 has only one variable from each omics dataset with a nonzero loading.
Note that in practice you might want to use more finely spaced numbers of features in the tuning grid. I used less finely spaced numbers to save time for the demo.
list.keepX <- tune.MULTI$choice.keepX
list.keepX## $transcripts
## [1] 40 1
##
## $metabolites
## [1] 50 1
##
## $proteins
## [1] 1 2Final DIABLO model
To review, we have done the following to decide on the parameters of our final DIABLO model:
- PLS on each pair of omics datasets to decide what cross-correlation
values to put in the design matrix (
design) - Cross-validation on an initial DIABLO model to decide how many
components to include in the final model (
ncomp) - Another cross-validation after setting the optimal number of
components, to decide how many features from each omics dataset will
have nonzero loadings on each component in the final model
(
list.keepX)
Input all those values into the block.splsda() function
to fit our final DIABLO model!
final.diablo.model <- block.splsda(X = data, Y = Y, ncomp = 2,
keepX = list.keepX, design = design)Look at selected variables
Which variables from each of the omics datasets had non-zero loadings
on the two components? The function selectVar() pulls these
out. The selected variables have patterns of differential abundance that
maximizes the covariances between all pairs of datasets and the response
variable (time point).
Note that the stability of selected variables, or how consistently
the same variables were selected across multiple cross-validation runs,
can be extracted from the output of the perf()
function.
selectVar(final.diablo.model, block = 'transcripts', comp = 1)$transcripts$name ## [1] "Dgat2" "Mmp11" "Fasn"
## [4] "Pum3" "Nectin2" "Slc5a3"
## [7] "Insig1" "Fads1" "Parm1"
## [10] "Mki67" "Trib1" "Pfkfb1"
## [13] "Slc16a10" "Pals2" "ENSMUSG00000089901"
## [16] "Megf9" "Nudt4" "Cyb5b"
## [19] "Net1" "D5Ertd579e" "Pdzd2"
## [22] "Srf" "Txndc12" "Apln"
## [25] "Cavin2" "Lrrc39" "Itpr1"
## [28] "Plin1" "Usp2" "Aacs"
## [31] "Rasa3" "Bltp3b" "Cry2"
## [34] "Pmm1" "Dbp" "Pex14"
## [37] "Gypc" "Vps13c" "Mid1ip1"
## [40] "Snta1"selectVar(final.diablo.model, block = 'metabolites', comp = 1)$metabolites$name ## [1] "3-hydroxyglutaric acid" "acetoacetic acid"
## [3] "l-malic acid" "uric acid"
## [5] "pantothenic acid" "galactose-1-phosphate"
## [7] "2-ketoglutaric acid" "l-leucine"
## [9] "o-butanoylcarnitine" "citric acid"
## [11] "l-aspartic acid" "fumaric acid"
## [13] "udp-n-acetylglucosamine" "glucose-1-phosphate"
## [15] "(2s)-2-azaniumylpropanoate" "l-phenylalanine"
## [17] "5-methoxy-5-oxopentanoic acid" "oleoylcarnitine"
## [19] "l-methionine" "nicotinamide adenine dinucleotide"
## [21] "4-oxo-l-proline" "succinic acid"
## [23] "hypotaurine" "d-fructose 6-phosphoric acid"
## [25] "l-tyrosine" "riboflavin"
## [27] "hypoxanthine" "l-proline"
## [29] "d-ribofuranose 5-phosphate" "n-acetyl-l-serine"
## [31] "oleic acid" "creatine"
## [33] "l-glutamine" "l-lactic acid"
## [35] "nadide" "3-hydroxyphenylacetic acid"
## [37] "barbituric acid" "aica ribonucleotide"
## [39] "thymidine" "n-acetyl-l-glutamic acid"
## [41] "arachidonic acid" "5-hydroxyindole-3-acetic acid"
## [43] "3-indoleacrylic acid" "sedoheptulose 7-phosphate"
## [45] "n-methylnicotinamide" "l-palmitoylcarnitine"
## [47] "l-carnitine" "d-erythrose 4-phosphate"
## [49] "isovalerylcarnitine" "4-methyl-2-oxopentanoic acid"selectVar(final.diablo.model, block = 'proteins', comp = 1)$proteins$name## [1] "Q3THX5"selectVar(final.diablo.model, block = 'transcripts', comp = 2)$transcripts$name ## [1] "Adrb3"selectVar(final.diablo.model, block = 'metabolites', comp = 2)$metabolites$name ## [1] "inosine"selectVar(final.diablo.model, block = 'proteins', comp = 2)$proteins$name ## [1] "P21956-2" "Q9Z2I8"Plots
Diagnostic correlation plot
The first thing to do is a diagnostic correlation plot that confirms that the correlations between each pair of omics were maximized as specified by the design matrix. The correlations are high and positive which looks good.
plotDiablo(final.diablo.model, ncomp = 1, col.per.group = custom_colors)
Individual scatter plots and arrow plot
Next, we can plot scatter plots of how the individual samples are represented in the two-dimensional space formed by the two components. A separate plot is generated for each of the three omics (or blocks, as they are referred to in the mixOmics package).
The plot here shows that across all three omics, component 1 separates the 0h time point from the other two. Despite only having a few variables with nonzero loadings, component 2 further separates the 4h and 24h time points. The separation is cleanest for transcriptome but the other two omics also separate the time points well.
plotIndiv(final.diablo.model, ind.names = FALSE, legend = TRUE,
title = 'DIABLO Sample Plots', col = custom_colors)
A so-called “arrow plot” shows how well the three omics agree with
one another. This plot has the same information as the plot we just made
by plotIndiv() but with a slightly different twist: each of
the 18 samples is shown as a central point with three arrows radiating
from it, each one pointing to the sample’s coordinate on the two PLS
components with respect to one of the three omics. If the arrows are
short and the representations of each sample in the three different
omics are close to one another, that means the omics datasets agree with
each other closely. As we would expect given the high (> 0.9)
correlations between each pair of omics, we see that the transcriptome,
metabolome, and proteome agree well with each other here. Regardless of
which omics dataset you look at, the clustering pattern of the time
points is virtually identical.
plotArrow(final.diablo.model, ind.names = FALSE, legend = TRUE,
title = 'DIABLO', col = custom_colors)
Correlation circle plot
This is a plot of the correlations between each of the selected variables and the components. Recall that 40, 50, and 1 variables were selected from transcriptome, metabolome, and proteome, respectively, for component 1, and only 1 or 2 from each of the omics for component 2. The transcriptome and proteome features are nearly all perfectly positively or negatively correlated with the component they were selected on. This is expected. The metabolome features are not as cleanly correlated.
plotVar(final.diablo.model, var.names = FALSE,
style = 'lattice', legend = TRUE,
pch = c(16, 17, 15), cex = c(0.8, 0.8, 0.8),
col = custom_colors)
Circos plot
A Circos plot is a fairly common type of visualization that you see in omics integration studies. I am not sure how effective they are at conveying information about the data, because they are somewhat unintuitive. But they are very information-rich, that is for sure.
Let’s break down what we are looking at here.
- You can see a slice of the pie for each of the three omics. Most of the circle is taken up by the transcriptome and metabolome features because 40 to 50 variables were selected from each of them to have nonzero loadings in the final model, compared to only 3 total for proteome.
- Around the outside are lines indicating the relative abundance of each feature at the three time points (there’s no labeled y-axis because it’s all just relative, but the further from the center, the higher the abundance). Of the selected genes from the transcriptome dataset, most were abundant at time 0h and less abundant at times 4h and 24h, and the rest showed the opposite pattern. This matches what we saw in the individual scatterplots: component 1 separates 0h from the other two time points, so it makes sense that the variables selected for that component are the ones where the 0h time point differs the most from the other two. There is more variability in the abundance pattern of the selected metabolites.
- The spiderweb-like pattern in the middle consists of lines
connecting pairs of features that are (a) in two different omics and (b)
correlated with one another. I set the argument
cutoff = 0.97to display a line at the very strict cutoff of \(|r| \ge 0.97\). Even with this strict cutoff, several of the metabolites and all of the proteins are shown as being correlated with many of the genes in the transcriptome. The pairs that are most abundant at the same time point(s) are connected with purple lines (indicating positive correlation) and black lines (indicating negative correlation) connect pairs of features that are abundant at different time points.
circosPlot(final.diablo.model, cutoff = 0.97, line = TRUE,
color.blocks = custom_colors,
color.cor = c("mediumpurple1","grey20"), size.labels = 1,
size.variables = 0.4, linkWidth = 0.3, showIntraLinks = FALSE, var.adj = -0.33)
It is possible to make Circos plots for only individual components,
for example by adding comp = 1 as an argument. There are
many other arguments that you can read about in the function
documentation.
Loading weights plots
This plot shows the selected variables for one or more of the omics, as a horizontal bar plot. The bar for each variable shows the group that has the highest absolute value loading on component 1.
First let’s make a transcriptome loadings plot. We see that the
selected genes have either a positive or a negative loading on component
1. Because positive loadings on component 1 happen to be associated with
high abundance at the 4h time point, the plotLoadings()
function automatically colors the genes with the 4h color (green) if
they have a positive loading on component 1. The genes with a negative
loading are all colored pink because they were more abundant at the 0h
time point.
plotLoadings(final.diablo.model,
comp = 1, block = 'transcripts',
contrib = 'max', method = 'median',
size.name = 0.6, legend.color = custom_colors)
Here is a similar plot for metabolites. Most of the selected variables had highest abundance at the 0h time point.
plotLoadings(final.diablo.model,
comp = 1, block = 'metabolites',
contrib = 'max', method = 'median',
size.name = 0.6, legend.color = custom_colors)
Only one variable was selected from the proteome dataset for component 1, so we aren’t making a loadings plot for proteome.
Heat map
This heat map has a similar interpretation to the one we made in a
previous session. In the heat map, the colors represent standardized
abundance values for the features. Each column is a feature that was
selected to have nonzero loading in the final DIABLO model. Each row is
one of the 18 samples. If we don’t supply a comp argument
to the heat map plotting function cimDiablo(), it will
default to showing us component 1. So as with some of the previous
visualizations we’ve seen, this heat map once again emphasizes that
component 1 separates out the 0h time point from the other two.
Essentially all features shown here are either high in abundance at 0h
and low at 4h and 24h (the columns on the left), or the opposite (the
columns on the right).
color_palette <- colorRampPalette(c("turquoise4", "turquoise3", "paleturquoise1", "lemonchiffon", "peachpuff", "salmon", "darkred"))(100)
cimDiablo(final.diablo.model, margins = c(2, 18), size.legend = 1, color.blocks = custom_colors, color.Y = c("gray60", "gray90", "gray70"), color = color_palette)##
## trimming values to [-3, 3] range for cim visualisation. See 'trim' arg in ?cimDiablo
Performance of the model
We can now look at how well the model predicts. In other words, working “backward” from the omics data, can we predict what time point a particular observation of omics data comes from? This is once again done by cross-validation, where we purposefully hide some data from the model fitting process, fit the model to the remaining data, and predict what groups the held-out data points belong to.
We see that the error rate is relatively low with one component (~4% overall) and the prediction performance improves even more when adding the second component.
perf.diablo <- perf(final.diablo.model, validation = 'Mfold',
folds = 3, nrepeat = 10,
dist = 'centroids.dist', seed = 321)
perf.diablo$WeightedVote.error.rate## $centroids.dist
## comp1 comp2
## 0h 0.00000000 0.000000000
## 4h 0.05000000 0.000000000
## 24h 0.06666667 0.016666667
## Overall.ER 0.03888889 0.005555556
## Overall.BER 0.03888889 0.005555556It may not make a lot of sense why we want to check this, because we already know the time point for all the data that we have, and in practice you would not have a sample of omics data without knowing what time point it came from. After all, our goal here is more exploring patterns in existing omics data, than training a model so that we can predict outcomes for future omics data we collect. But even if our goal is inference and understanding rather than prediction, it is reassuring if our model has good predictive performance. We can say that our model has learned something from the data, which means that we can look at outputs from the model and use them to understand the data. The data are from the real world and the model is not, it’s only an abstraction of reality. So it’s only useful if we can tie it back to reality in some way.
On that inspiring note, let’s end our workshop on multiomics integration!