add MFA to `test_pathway` function

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

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)

NEWS.md

-# phoenics 0.4 [2024-09-23]
+# phoenics 0.4 [2024-11-26]
+
+## New features
+
+* added MFA to calculate pathway scores in `pathway_test`
 
 ## Improvement:
 

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)
   })

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)
   })
   

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)

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
  }
 }