The TARA Oceans expedition facilitated the study of plankton communities by providing ocean metagenomic data combined with environmental measures to the scientific community. This study focuses on 139 prokaryotic-enriched samples collected from 68 stations and spread across three depth layers: the surface (SRF), the deep chlorophyll maximum (DCM) layer and the mesopelagic (MES) zones. Samples were located in 8 different oceans or seas: Indian Ocean (IO), Mediterranean Sea (MS), North Atlantic Ocean (NAO), North Pacific Ocean (NPO), Red Sea (RS), South Atlantic Ocean (SAO), South Pacific Ocean (SPO) and South Ocean (SO).
In this vignette, we consider a subset of the original data analyzed in the article (Mariette & Villa-Vialaneix, 2018) and illustrate the usefulness of mixKernel to 1/ perform an integrative exploratory analysis as in (Mariette & Villa-Vialaneix, 2018) and to 2/ select relevant variables for unsupervised analysis.
The data include 1% of the 35,650 prokaryotic OTUs and of the 39,246 bacterial genes that were randomly selected. The aim is to integrate prokaryotic abundances and functional processes to environmental measure with an unsupervised method.
Install and load the mixOmics and mixKernel packages:
The (previously normalized) datasets are provided as matrices with matching sample names (rownames):
data(TARAoceans)
# more details with: ?TARAOceans
# we check the dimension of the data:
lapply(list("phychem" = TARAoceans$phychem, "pro.phylo" = TARAoceans$pro.phylo,
"pro.NOGs" = TARAoceans$pro.NOGs), dim)
## $phychem
## [1] 139 22
##
## $pro.phylo
## [1] 139 356
##
## $pro.NOGs
## [1] 139 638
For each input dataset, a kernel is computed using the function
compute.kernel
with the choice of linear, phylogenic or
abundance kernels. A user defined function can also be provided as
input(argument kernel.func
, see
?compute.kernel
).
The results are lists with a ‘kernel’ entry that stores the kernel matrix. The resulting kernels are symmetric matrices with a size equal to the number of observations (rows) in the input datasets.
phychem.kernel <- compute.kernel(TARAoceans$phychem, kernel.func = "linear")
pro.phylo.kernel <- compute.kernel(TARAoceans$pro.phylo, kernel.func = "abundance")
pro.NOGs.kernel <- compute.kernel(TARAoceans$pro.NOGs, kernel.func = "abundance")
# check dimensions
dim(pro.NOGs.kernel$kernel)
## [1] 139 139
A general overview of the correlation structure between datasets is
obtained as described in Mariette and Villa-Vialaneix (2018) and
displayed using the function cim.kernel
:
cim.kernel(phychem = phychem.kernel,
pro.phylo = pro.phylo.kernel,
pro.NOGs = pro.NOGs.kernel,
method = "square")
The figure shows that pro.phylo
and
pro.NOGs
is the most correlated pair of kernels. This
result is expected as both kernels provide a summary of prokaryotic
communities.
The function combine.kernels
implements 3 different
methods for combining kernels: STATIS-UMKL, sparse-UMKL and full-UMKL
(see more details in Mariette and Villa-Vialaneix, 2018). It returns a
meta-kernel that can be used as an input for the function
kernel.pca
(kernel PCA). The three methods bring
complementary information and must be chosen according to the research
question.
The STATIS-UMKL
approach gives an overview on the common
information between the different datasets. The full-UMKL
computes a kernel that minimizes the distortion between all input
kernels. The sparse-UMKL
is a sparse variant of
full-UMKL
that selects the most relevant kernels in
addition to distortion minimization.
A kernel PCA can be performed from the combined kernel with the
function kernel.pca``. The argument
ncomp` allows to choose
how many components to extract from KPCA.
Sample plots using the plotIndiv
function from
mixOmics
:
all.depths <- levels(factor(TARAoceans$sample$depth))
depth.pch <- c(20, 17, 4, 3)[match(TARAoceans$sample$depth, all.depths)]
plotIndiv(kernel.pca.result,
comp = c(1, 2),
ind.names = FALSE,
legend = TRUE,
group = as.vector(TARAoceans$sample$ocean),
col.per.group = c("#f99943", "#44a7c4", "#05b052", "#2f6395",
"#bb5352", "#87c242", "#07080a", "#92bbdb"),
pch = depth.pch,
pch.levels = TARAoceans$sample$depth,
legend.title = "Ocean / Sea",
title = "Projection of TARA Oceans stations",
size.title = 10,
legend.title.pch = "Depth")
The explained variance supported by each axis of KPCA is displayed
with the plot
function, and can help choosing the number of
components in KPCA.
The first axis summarizes ~ 20% of the total variance.
Here we focus on the information summarized on the first component.
Variable values are randomly permuted with the function
permute.kernel.pca
.
In the following example, physical variable are permuted at the
variable level (kernel phychem
), OTU abundances from
pro.phylo
kernel are permuted at the phylum level (OTU
phyla are stored in the second column, named Phylum
, of the
taxonomy annotation provided in TARAoceans
object in the
entry taxonomy
) and gene abundances from
pro.NOGs
are permuted at the GO level (GO are provided in
the entry GO
of the dataset):
## [1] Proteobacteria Proteobacteria Proteobacteria Proteobacteria Proteobacteria
## [6] Cyanobacteria Proteobacteria Proteobacteria Chloroflexi Proteobacteria
## 56 Levels: Acidobacteria Actinobacteria aquifer1 Aquificae ... WCHB1-60
## [1] NA NA "K" NA NA "S" "S" "S" NA "S"
# here we set a seed for reproducible results with this tutorial
set.seed(17051753)
kernel.pca.result <- kernel.pca.permute(kernel.pca.result, ncomp = 1,
phychem = colnames(TARAoceans$phychem),
pro.phylo = TARAoceans$taxonomy[, "Phylum"],
pro.NOGs = TARAoceans$GO)
Results are displayed with the function
plotVar.kernel.pca
. The argument ndisplay
indicates the number of variables to display for each kernel:
Proteobacteria
is the most important variable for the
pro.phylo
kernel.
The relative abundance of `Proteobacteria`` is then extracted in each of our 139 samples, and each sample is colored according to the value of this variable in the KPCA projection plot:
selected <- which(TARAoceans$taxonomy[, "Phylum"] == "Proteobacteria")
proteobacteria.per.sample <- apply(TARAoceans$pro.phylo[, selected], 1, sum) /
apply(TARAoceans$pro.phylo, 1, sum)
colfunc <- colorRampPalette(c("royalblue", "red"))
col.proteo <- colfunc(length(proteobacteria.per.sample))
col.proteo <- col.proteo[rank(proteobacteria.per.sample, ties = "first")]
plotIndiv(kernel.pca.result,
comp = c(1, 2),
ind.names = FALSE,
legend = FALSE,
col = col.proteo,
pch = depth.pch,
pch.levels = TARAoceans$sample$depth,
legend.title = "Ocean / Sea",
title = "Representation of Proteobacteria abundance",
legend.title.pch = "Depth")
Similarly, the temperature is the most important variable for the
phychem
kernel. The temperature values can be displayed on
the kernel PCA projection as follows:
col.temp <- colfunc(length(TARAoceans$phychem[, 4]))
col.temp <- col.temp[rank(TARAoceans$phychem[, 4], ties = "first")]
plotIndiv(kernel.pca.result,
comp = c(1, 2),
ind.names = FALSE,
legend = FALSE,
col = col.temp,
pch = depth.pch,
pch.levels = TARAoceans$sample$depth,
legend.title = "Ocean / Sea",
title = "Representation of mean temperature",
legend.title.pch = "Depth")