diff --git a/DESCRIPTION b/DESCRIPTION
index 4b5ba23cbcea4498c818378f23a4eba2e9ff8660..cff3bef59a083a1064a9a0bcd04121c0fd1ba037 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: phoenics
Type: Package
Title: Pathways Longitudinal and Differential Analysis in Metabolomics
Version: 0.4
-Date: 2024-09-23
+Date: 2024-11-29
Depends: R (>= 4.0.0)
Imports: lme4, blme, tibble, FactoMineR, factoextra, tidyr, stats
Suggests: KEGGREST
diff --git a/NAMESPACE b/NAMESPACE
index 68a007efd95b14c6563a7aa31014d056acce4978..026f39e151b1492f34a8ac387c8f4e17b203aed0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,9 +13,12 @@ export(head)
export(overlap_coefficient)
export(pathway_search)
export(test_pathway)
+importFrom(FactoMineR,MFA)
importFrom(FactoMineR,PCA)
importFrom(blme,blmer)
importFrom(factoextra,fviz_eig)
+importFrom(factoextra,fviz_mfa_ind)
+importFrom(factoextra,fviz_mfa_var)
importFrom(factoextra,fviz_pca_ind)
importFrom(factoextra,fviz_pca_var)
importFrom(lme4,lmer)
@@ -26,5 +29,6 @@ importFrom(stats,p.adjust.methods)
importFrom(tibble,column_to_rownames)
importFrom(tibble,remove_rownames)
importFrom(tibble,rownames_to_column)
+importFrom(tidyr,pivot_wider)
importFrom(tidyr,separate)
importFrom(utils,read.table)
diff --git a/NEWS.md b/NEWS.md
index a8c313b5a98d842b94d32dbd91a84d0cebceff14..608aa340845aa997881564c241a3f3064c28f78a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,8 @@
-# phoenics 0.4 [2024-09-23]
+# phoenics 0.4 [2024-11-26]
+
+## New features
+
+* added MFA to calculate pathway scores in `pathway_test`
## Improvement:
diff --git a/R/pathwayRes-class.R b/R/pathwayRes-class.R
index f0264458ab8ff46ba3370cb8d2db600ae0188f0d..b79d8fe7197e510c5d3013c5d8d883c92ef8c8cd 100644
--- a/R/pathwayRes-class.R
+++ b/R/pathwayRes-class.R
@@ -16,8 +16,12 @@ NULL
#' quantif <- from_ASICS_to_PHOENICS(quantif)
#' out_test <- test_pathway(quantif, design, pathways,
#' fixed = c("Age", "Treatment"), random = "Mouse",
-#' npc = 2, model = "blmer")
+#' analysis = "PCA", npc = 2, model = "blmer")
#' summary(out_test)
+#' out_test2 <- test_pathway(quantif, design, pathways,
+#' fixed = c("Age", "Treatment"), random = "Mouse",
+#' analysis = "MFA", npc = 2, model = "blmer")
+#' summary(out_test2)
#' @exportS3Method
summary.pathwayRes <- function(object, ...) {
cat("Number of pathways: ", length(object), "\n")
@@ -36,29 +40,39 @@ print.pathwayRes <- function(x, ...) {
cat("$", names(path_res)[2], ": ", "Pathway code \n", sep = "")
cat("$", names(path_res)[3], ": ", "Metabolites in the pathway \n", sep = "")
cat("$", names(path_res)[4], ": ", "PCA results \n", sep = "")
- cat("$", names(path_res)[5], ": ", "Model results \n", sep = "")
- cat("$", names(path_res)[6], ": ", "Pathway test results \n", sep = "")
+ cat("$", names(path_res)[5], ": ", "MFA results \n", sep = "")
+ cat("$", names(path_res)[6], ": ", "Model results \n", sep = "")
+ cat("$", names(path_res)[7], ": ", "Pathway test results \n", sep = "")
}
#' @rdname pathwayRes
#' @aliases plot.pathwayRes
#' @param pathway_id a character string or vector of pathway codes or names.
#' Default to \code{NULL} in which case the first pathway is displayed
+#' @param analysis a character string indicating the type of analysis to display
#' @param plot a character string indicating the type of plot to return. Default
-#' to \code{"eig"} (the screegraph of the PCA is displayed)
+#' to \code{"eig"} (the screegraph of the PCA or MFA is displayed)
#' @param habillage a character string indicating the column of the design used
#' to color the individuals. Only used when \code{plot = "ind"}. Default to
#' \code{"none"} (no color)
-#' @importFrom factoextra fviz_pca_ind fviz_eig fviz_pca_var
+#' @importFrom factoextra fviz_pca_ind fviz_eig fviz_pca_var fviz_mfa_ind fviz_mfa_var
#' @examples
-#' plot(out_test)
-#' plot(out_test, pathway_id = "mmu00052", plot = "var")
-#' plot(out_test, pathway_id = "mmu00052", plot = "ind", habillage = "Age")
+#' plot(out_test, analysis = "PCA")
+#' plot(out_test, pathway_id = "mmu00052", analysis = "PCA", plot = "var")
+#' plot(out_test, pathway_id = "mmu00052", analysis = "PCA", plot = "ind",
+#' habillage = "Age")
+#' plot(out_test2, pathway_id = "mmu00562", analysis = "MFA", plot = "eig")
+#' plot(out_test2, pathway_id = "mmu00562", analysis = "MFA", plot = "var")
+#' plot(out_test2, pathway_id = "mmu00562", analysis = "MFA", plot = "ind")
+#' plot(out_test2, pathway_id = "mmu00562", analysis = "MFA", plot = "group")
#' @exportS3Method
plot.pathwayRes <- function(x, ..., pathway_id = NULL,
- plot = c("eig", "var", "ind"), habillage = "none") {
+ analysis = c("PCA", "MFA"),
+ plot = c("eig", "var", "ind", "group"),
+ habillage = "none") {
plot <- match.arg(plot)
+ analysis <- match.arg(analysis)
if (is.null(pathway_id)) {
message("`pathway_id` not specified... defaulting to the first pathway")
@@ -66,20 +80,45 @@ plot.pathwayRes <- function(x, ..., pathway_id = NULL,
}
object_sub <- extract.pathwayRes(x, pathway_id)
+
lapply(object_sub, function (path_res) {
- if (plot == "var") {
- p <- factoextra::fviz_pca_var(path_res$PCA, title = path_res$pathway_name,
- repel = TRUE)
- }
- if (plot == "ind") {
- p <- factoextra::fviz_pca_ind(path_res$PCA, habillage = habillage,
- title = path_res$pathway_name,
- geom = "point", pointsize = 2,
- mean.point = FALSE)
+ if (analysis == "PCA") {
+ if (plot == "var") {
+ p <- factoextra::fviz_pca_var(path_res$PCA, title = path_res$pathway_name,
+ repel = TRUE)
+ }
+ if (plot == "ind") {
+ p <- factoextra::fviz_pca_ind(path_res$PCA, habillage = habillage,
+ title = path_res$pathway_name,
+ geom = "point", pointsize = 2,
+ mean.point = FALSE)
+ }
+ if (plot == "eig") {
+ p <- factoextra::fviz_eig(path_res$PCA, addlabels = TRUE,
+ title = path_res$pathway_name)
+ }
}
- if (plot == "eig") {
- p <- factoextra::fviz_eig(path_res$PCA, addlabels = TRUE,
- title = path_res$pathway_name)
+ if (analysis == "MFA") {
+ if (plot == "group") {
+ p <- factoextra::fviz_mfa_var(path_res$MFA, choice = "group",
+ title = path_res$pathway_name,
+ repel = TRUE)
+ }
+ if (plot == "var") {
+ p <- factoextra::fviz_mfa_var(path_res$MFA, choice = "quanti.var",
+ title = path_res$pathway_name,
+ repel = TRUE)
+ }
+ if (plot == "ind") {
+ p <- factoextra::fviz_mfa_ind(path_res$MFA, partial = "all",
+ habillage = habillage,
+ title = path_res$pathway_name,
+ pointsize = 2)
+ }
+ if (plot == "eig") {
+ p <- factoextra::fviz_eig(path_res$MFA, addlabels = TRUE,
+ title = path_res$pathway_name)
+ }
}
return(p)
})
diff --git a/R/test_pathway.R b/R/test_pathway.R
index fedac69f6ed656af64f2a080139437dfbd65695d..7231e885724dd1d7d0f580f56aaa5ff61724fe65 100644
--- a/R/test_pathway.R
+++ b/R/test_pathway.R
@@ -2,7 +2,7 @@
#'
#' @description Perform a differential analysis at pathway level based on
#' metabolite quantifications and information on pathway metabolite composition.
-#' The method relies on a PCA step.
+#' The method relies on a PCA or a MFA step.
#'
#' @param quantif data.frame or matrix of the metabolite quantification with
#' samples in rows (sample identifiers must be row names) and metabolites in
@@ -12,18 +12,21 @@
#' in columns. Column names must be provided and are used as arguments for
#' \code{fixed} and \code{random}
#' @param pathways data.frame or matrix with metabolites in rows (all
-#' metabolites in columns of \code{quantif} must have a row in this input) and the
-#' following information in columns: \itemize{
+#' metabolites in columns of \code{quantif} must have a row in this input) and
+#' the following information in columns: \itemize{
#' \item \code{metabolite_code} metabolite code
#' \item \code{metabolite_name} metabolite name
#' \item \code{pathway_code} pathway code (identifier)
#' \item \code{pathway_name} name of the pathway
#' }
#' @param npc number of PCs for the PCA step
+#' @param analysis character string indicating if the pathway scores are
+#' obtained using \link[FactoMineR]{PCA} or \link[FactoMineR]{MFA}
#' @param model a character string indicating if the model has to be fitted
#' using \link[lme4]{lmer} or \link[blme]{blmer}. Default to \code{"lmer"}
#' @param fixed character vector of fixed effects to be included in the model.
-#' They must correspond to column names of \code{design}
+#' They must correspond to column names of \code{design}. If \code{analysis} is
+#' "MFA", the first fixed effect must correspond to the time effect
#' @param random character vector of random effects to be included in the model.
#' They must correspond to column names of \code{design}
#' @param organism organism code in KEGG database. Required if
@@ -38,32 +41,42 @@
#' \item{metabolites}{a data.frame with the names and codes of the quantified
#' metabolites in the pathway}
#' \item{PCA}{the result of the pathway PCA (a \code{PCA} object as obtained
-#' from \link[FactoMineR]{PCA})}
+#' from \link[FactoMineR]{PCA}), is equal to \code{NULL} if
+#' \code{analysis = "MFA"}}
+#' \item{MFA}{the result of the pathway MFA (a \code{MFA} object as obtained
+#' from \link[FactoMineR]{MFA}), is equal to \code{NULL} if
+#' \code{analysis = "PCA"}}
#' \item{model}{the output of the mixed model fit}
#' \item{test_pathway}{a data.frame with the p-values for each tested fixed
#' effect}
#'
-#' @importFrom FactoMineR PCA
+#' @importFrom FactoMineR PCA MFA
#' @importFrom lme4 lmer
#' @importFrom blme blmer
#' @importFrom tibble column_to_rownames rownames_to_column
#' @importFrom stats as.formula anova
-#' @importFrom tidyr separate
+#' @importFrom tidyr separate pivot_wider
#'
#' @examples
#' data("MTBLS422")
#' quantif <- from_ASICS_to_PHOENICS(quantif)
#' out_test <- test_pathway(quantif, design, pathways,
#' fixed = c("Age", "Treatment"), random = "Mouse",
-#' npc = 2, model = "blmer")
+#' npc = 2, analysis = "PCA", model = "blmer")
#' out_test
#'
+#' out_test2 <- test_pathway(quantif, design, pathways,
+#' fixed = c("Age", "Treatment"), random = "Mouse",
+#' npc = 2, analysis = "MFA", model = "blmer")
+#' out_test2
+#'
#' if (requireNamespace("KEGGREST", quietly = TRUE)) {
#' \donttest{
-#' out_test2 <- test_pathway(quantif, design, pathways = "auto",
+#' out_test3 <- test_pathway(quantif, design, pathways = "auto",
#' fixed = c("Age", "Treatment"), random = "Mouse",
-#' npc = 2, model = "blmer", organism = "mmu")
-#' out_test2
+#' npc = 2, analysis = "PCA", model = "blmer",
+#' organism = "mmu")
+#' out_test3
#' }
#' }
#'
@@ -85,8 +98,9 @@
#' @export
test_pathway <- function(quantif, design, pathways = "auto", fixed, random,
- npc = 1L, model = c("lmer", "blmer"),
- organism = NULL, min_size = 2) {
+ npc = 1L, analysis = c("PCA", "MFA"),
+ model = c("lmer", "blmer"), organism = NULL,
+ min_size = 2) {
# Inputs tests
if (!is.character(fixed)) {
@@ -105,7 +119,6 @@ test_pathway <- function(quantif, design, pathways = "auto", fixed, random,
stop("Random effect must be in `design` column names.")
}
-
if (!identical(sort(rownames(quantif)), sort(rownames(design)))) {
stop("`quantif` and `design` row names must be identical.")
}
@@ -118,6 +131,7 @@ test_pathway <- function(quantif, design, pathways = "auto", fixed, random,
}
model <- match.arg(model)
+ analysis <- match.arg(analysis)
if (round(npc) != npc) {
stop("`npc` must be an integer.")
@@ -188,16 +202,48 @@ test_pathway <- function(quantif, design, pathways = "auto", fixed, random,
quantif_path <- merge(quantif_path, design, by = 0)
quantif_path <- column_to_rownames(quantif_path, "Row.names")
- res_pca <- PCA(quantif_path, graph = FALSE, ncp = npc,
- quali.sup = colnames(design))
+ if (analysis == "PCA") {
+ res_pca <- PCA(quantif_path, graph = FALSE, ncp = npc,
+ quali.sup = colnames(design))
+
+ score <- data.frame(res_pca$ind$coord[, 1:npc, drop = FALSE])
+ colnames(score) <- paste0(path, "_PC", 1:ncol(score))
+
+ score_path <- merge(score, design, by = 0)
+ score_path <- column_to_rownames(score_path, "Row.names")
+ }
- score <- data.frame(res_pca$ind$coord[, 1:npc, drop = FALSE])
- colnames(score) <- paste0(path, "_PC", 1:ncol(score))
+ if (analysis == "MFA") {
+ nb_metab <- length(metab_path)
+ timepoints <- unique(quantif_path[, fixed[1]])
+
+ quantif_path_wide <- pivot_wider(quantif_path,
+ id_cols = c(fixed[-1], random),
+ names_from = fixed[1],
+ values_from = metab_path,
+ names_vary = "slowest")
+ quantif_path_wide <- as.data.frame(quantif_path_wide)
+ rownames(quantif_path_wide) <- quantif_path_wide[, random[1]]
+
+ res_mfa <- MFA(quantif_path_wide,
+ group = c(length(c(fixed[-1], random)),
+ rep(nb_metab, length(timepoints))),
+ type = c("n", rep("s", length(timepoints))),
+ name.group = c("supplementary",
+ paste0("date_", timepoints)),
+ num.group.sup = 1, graph = FALSE)
+
+ score <- data.frame(res_mfa$ind$coord.partiel[, 1:npc, drop = FALSE])
+ colnames(score) <- paste0(path, "_PC", 1:ncol(score))
+
+ score_path <- rownames_to_column(score, "ind")
+ score_path <- separate(score_path, col = "ind",
+ into = c(random[1], fixed[1]), sep = ".date_")
+ score_path <- merge(unique(design[, c(fixed[-1], random)]), score_path,
+ by = c(random[1]))
+ }
# Mixed model
- score_path <- merge(score, design, by = 0)
- score_path <- column_to_rownames(score_path, "Row.names")
-
form <- paste0("~ ", paste(fixed, collapse = " + "), " + ",
paste(paste0("(1|", random, ")"), collapse = " + "))
fmla <- sapply(paste0(colnames(score), form), as.formula)
@@ -262,9 +308,17 @@ test_pathway <- function(quantif, design, pathways = "auto", fixed, random,
path_code <- unique(pathways_sub[, "pathway_code"])
out_met <- unique(pathways_sub[, c("metabolite_code", "metabolite_name")])
- res <- list("pathway_name" = path_name, "pathway_code" = path_code,
- "metabolites" = out_met, "PCA" = res_pca, "model" = res_lmm,
- "test_pathway" = res_test)
+ if (analysis == "PCA") {
+ res <- list("pathway_name" = path_name, "pathway_code" = path_code,
+ "metabolites" = out_met, "PCA" = res_pca, "MFA" = NULL,
+ "model" = res_lmm, "test_pathway" = res_test)
+ }
+ if (analysis == "MFA") {
+ res <- list("pathway_name" = path_name, "pathway_code" = path_code,
+ "metabolites" = out_met, "PCA" = NULL, "MFA" = res_mfa,
+ "model" = res_lmm, "test_pathway" = res_test)
+ }
+
return(res)
})
diff --git a/man/pathwayRes.Rd b/man/pathwayRes.Rd
index b57c12868d606c23b260a020e539cb7ce4eafdd7..3e2089c1e2f9f288660d515f071e4ee8f1a4af43 100644
--- a/man/pathwayRes.Rd
+++ b/man/pathwayRes.Rd
@@ -21,7 +21,8 @@
x,
...,
pathway_id = NULL,
- plot = c("eig", "var", "ind"),
+ analysis = c("PCA", "MFA"),
+ plot = c("eig", "var", "ind", "group"),
habillage = "none"
)
@@ -38,8 +39,10 @@ adjust_pval(object, method = p.adjust.methods, n = length(object))
\item{pathway_id}{a character string or vector of pathway codes or names}
+\item{analysis}{a character string indicating the type of analysis to display}
+
\item{plot}{a character string indicating the type of plot to return. Default
-to \code{"eig"} (the screegraph of the PCA is displayed)}
+to \code{"eig"} (the screegraph of the PCA or MFA is displayed)}
\item{habillage}{a character string indicating the column of the design used
to color the individuals. Only used when \code{plot = "ind"}. Default to
@@ -75,12 +78,21 @@ data("MTBLS422")
quantif <- from_ASICS_to_PHOENICS(quantif)
out_test <- test_pathway(quantif, design, pathways,
fixed = c("Age", "Treatment"), random = "Mouse",
- npc = 2, model = "blmer")
+ analysis = "PCA", npc = 2, model = "blmer")
summary(out_test)
+out_test2 <- test_pathway(quantif, design, pathways,
+ fixed = c("Age", "Treatment"), random = "Mouse",
+ analysis = "MFA", npc = 2, model = "blmer")
+summary(out_test2)
print(out_test)
-plot(out_test)
-plot(out_test, pathway_id = "mmu00052", plot = "var")
-plot(out_test, pathway_id = "mmu00052", plot = "ind", habillage = "Age")
+plot(out_test, analysis = "PCA")
+plot(out_test, pathway_id = "mmu00052", analysis = "PCA", plot = "var")
+plot(out_test, pathway_id = "mmu00052", analysis = "PCA", plot = "ind",
+ habillage = "Age")
+plot(out_test2, pathway_id = "mmu00562", analysis = "MFA", plot = "eig")
+plot(out_test2, pathway_id = "mmu00562", analysis = "MFA", plot = "var")
+plot(out_test2, pathway_id = "mmu00562", analysis = "MFA", plot = "ind")
+plot(out_test2, pathway_id = "mmu00562", analysis = "MFA", plot = "group")
head(out_test)
extract(out_test, "mmu00562")
adj_pval <- adjust_pval(out_test)
diff --git a/man/test_pathway.Rd b/man/test_pathway.Rd
index ea5f024087292bc55b9c196277b496af6d2c1f61..6efa704c0e3ffea9074e44a6cfcdca7a09ec7b9b 100644
--- a/man/test_pathway.Rd
+++ b/man/test_pathway.Rd
@@ -11,6 +11,7 @@ test_pathway(
fixed,
random,
npc = 1L,
+ analysis = c("PCA", "MFA"),
model = c("lmer", "blmer"),
organism = NULL,
min_size = 2
@@ -27,8 +28,8 @@ in columns. Column names must be provided and are used as arguments for
\code{fixed} and \code{random}}
\item{pathways}{data.frame or matrix with metabolites in rows (all
-metabolites in columns of \code{quantif} must have a row in this input) and the
-following information in columns: \itemize{
+metabolites in columns of \code{quantif} must have a row in this input) and
+the following information in columns: \itemize{
\item \code{metabolite_code} metabolite code
\item \code{metabolite_name} metabolite name
\item \code{pathway_code} pathway code (identifier)
@@ -36,13 +37,17 @@ following information in columns: \itemize{
}}
\item{fixed}{character vector of fixed effects to be included in the model.
-They must correspond to column names of \code{design}}
+They must correspond to column names of \code{design}. If \code{analysis} is
+"MFA", the first fixed effect must correspond to the time effect}
\item{random}{character vector of random effects to be included in the model.
They must correspond to column names of \code{design}}
\item{npc}{number of PCs for the PCA step}
+\item{analysis}{character string indicating if the pathway scores are
+obtained using \link[FactoMineR]{PCA} or \link[FactoMineR]{MFA}}
+
\item{model}{a character string indicating if the model has to be fitted
using \link[lme4]{lmer} or \link[blme]{blmer}. Default to \code{"lmer"}}
@@ -60,7 +65,11 @@ results. Each element of the list contain the following entries:
\item{metabolites}{a data.frame with the names and codes of the quantified
metabolites in the pathway}
\item{PCA}{the result of the pathway PCA (a \code{PCA} object as obtained
-from \link[FactoMineR]{PCA})}
+from \link[FactoMineR]{PCA}), is equal to \code{NULL} if
+\code{analysis = "MFA"}}
+\item{MFA}{the result of the pathway MFA (a \code{MFA} object as obtained
+from \link[FactoMineR]{MFA}), is equal to \code{NULL} if
+\code{analysis = "PCA"}}
\item{model}{the output of the mixed model fit}
\item{test_pathway}{a data.frame with the p-values for each tested fixed
effect}
@@ -68,7 +77,7 @@ effect}
\description{
Perform a differential analysis at pathway level based on
metabolite quantifications and information on pathway metabolite composition.
-The method relies on a PCA step.
+The method relies on a PCA or a MFA step.
}
\details{
If \code{pathways = "auto"}, information on pathways in which metabolites are
@@ -82,15 +91,21 @@ data("MTBLS422")
quantif <- from_ASICS_to_PHOENICS(quantif)
out_test <- test_pathway(quantif, design, pathways,
fixed = c("Age", "Treatment"), random = "Mouse",
- npc = 2, model = "blmer")
+ npc = 2, analysis = "PCA", model = "blmer")
out_test
+out_test2 <- test_pathway(quantif, design, pathways,
+ fixed = c("Age", "Treatment"), random = "Mouse",
+ npc = 2, analysis = "MFA", model = "blmer")
+out_test2
+
if (requireNamespace("KEGGREST", quietly = TRUE)) {
\donttest{
- out_test2 <- test_pathway(quantif, design, pathways = "auto",
+ out_test3 <- test_pathway(quantif, design, pathways = "auto",
fixed = c("Age", "Treatment"), random = "Mouse",
- npc = 2, model = "blmer", organism = "mmu")
- out_test2
+ npc = 2, analysis = "PCA", model = "blmer",
+ organism = "mmu")
+ out_test3
}
}