Refactor NewInterface.R to unified API design

Rewrote NewInterface.R to provide a simplified, unified API for microbiome data analysis, consolidating multiple plotting and testing functions into core unified functions (plot_taxa, plot_alpha, plot_beta, test_taxa, test_beta, test_alpha). The new interface automatically detects study design and routes to the appropriate underlying function, improving usability and maintainability. Added detailed documentation and parameter validation, and removed legacy and redundant code.

Author: cafferychen777

R/NewInterface.R

@@ -7,7 +7,6 @@
 #' @param alpha.name The alpha diversity index to be plotted. Supported indices include "shannon", "simpson", "observed_species", "chao1", "ace", "pielou", and "faith_pd".
 #' @param depth An integer specifying the sequencing depth for the "Rarefy" and "Rarefy-TSS" methods.
 #' If NULL, no rarefaction is performed.
-#' @param subject.var The variable in the metadata table that represents the subject.
 #' @param time.var The variable in the metadata table that represents the time.
 #' @param t.level Character string specifying the time level/value to subset data to,
 #' if a time variable is provided. Default NULL does not subset data.
@@ -76,7 +75,6 @@
 #'   data.obj     = peerj32.obj,
 #'   alpha.obj    = NULL,
 #'   alpha.name   = c("simpson"),
-#'   subject.var  = "subject",
 #'   time.var     = "time",
 #'   t.level      = "2",
 #'   group.var    = "group",
@@ -96,7 +94,6 @@
 #'   data.obj     = peerj32.obj,
 #'   alpha.obj    = alpha.obj,
 #'   alpha.name   = c("simpson"),
-#'   subject.var  = "subject",
 #'   time.var     = "time",
 #'   t.level      = "2",
 #'   group.var    = "group",
@@ -119,7 +116,6 @@
 #'   data.obj     = subset_T2D.obj,
 #'   alpha.obj    = NULL,
 #'   alpha.name   = c("shannon"),
-#'   subject.var  = "subject_id",
 #'   time.var     = "visit_number",
 #'   t.level      = "   3",
 #'   group.var    = "subject_race",
@@ -146,7 +142,6 @@
 #'   data.obj     = peerj32.obj,
 #'   alpha.obj    = NULL,
 #'   alpha.name   = c("simpson"),
-#'   subject.var  = "subject",
 #'   time.var     = "time",
 #'   t.level      = "2",
 #'   group.var    = "group",
@@ -167,7 +162,6 @@ generate_alpha_boxplot_single <- function (data.obj,
                                            alpha.name = c("shannon",
                                                           "observed_species"),
                                            depth = NULL,
-                                           subject.var = NULL,
                                            time.var = NULL,
                                            t.level = NULL,
                                            group.var = NULL,

R/generate_alpha_boxplot_single.R

@@ -23,12 +23,12 @@
 #' @param alpha.name The alpha diversity index to be plotted. Supported indices include "shannon", "simpson", "observed_species", "chao1", "ace", "pielou", and "faith_pd". Previously named as `alpha.index`.
 #' @param depth An integer specifying the sequencing depth for the "Rarefy" and "Rarefy-TSS" methods.
 #' If NULL, no rarefaction is performed.
-#' @param time.var Character string specifying the column name in metadata containing
-#'                time values for each sample. Required to identify pairs of time
-#'                points to calculate changes between.
 #' @param subject.var Character string specifying the column name in metadata containing
 #'                    unique subject IDs. Required to pair samples from the same subject
 #'                    across time points.
+#' @param time.var Character string specifying the column name in metadata containing
+#'                time values for each sample. Required to identify pairs of time
+#'                points to calculate changes between.
 #' @param group.var Character string specifying the column name in metadata containing
 #'                 grouping categories. Used as a predictor in the models to test for
 #'                 differences in changes between groups. Optional, can be NULL.
@@ -83,9 +83,9 @@
 #' generate_alpha_change_test_pair(
 #'   data.obj = peerj32.obj,
 #'   alpha.obj = NULL,
-#'   time.var = "time",
 #'   alpha.name = c("shannon"),
 #'   subject.var = "subject",
+#'   time.var = "time",
 #'   group.var = "sex",
 #'   adj.vars = NULL,
 #'   change.base = "2",
@@ -100,9 +100,9 @@
 #' generate_alpha_change_test_pair(
 #'   data.obj = peerj32.obj,
 #'   alpha.obj = NULL,
-#'   time.var = "Day",
 #'   alpha.name = c("shannon"),
 #'   subject.var = "subject",
+#'   time.var = "Day",
 #'   group.var = "sex",
 #'   adj.vars = c("group"),
 #'   change.base = "2",
@@ -115,9 +115,9 @@
 #' generate_alpha_change_test_pair(
 #'   data.obj = subset_pairs.obj,
 #'   alpha.obj = NULL,
-#'   time.var = "Antibiotic",
 #'   alpha.name = c("shannon"),
 #'   subject.var = "MouseID",
+#'   time.var = "Antibiotic",
 #'   group.var = "Sex",
 #'   adj.vars = NULL,
 #'   change.base = "Baseline",
@@ -131,8 +131,8 @@ generate_alpha_change_test_pair <-
            alpha.obj = NULL,
            alpha.name = NULL,
            depth = NULL,
-           time.var,
            subject.var,
+           time.var,
            group.var,
            adj.vars = NULL,
            change.base,

R/generate_alpha_change_test_pair.R

@@ -68,11 +68,11 @@ extract_coef <- function(model) {
 #' "observed_species", "chao1", "ace", "pielou", and "faith_pd".
 #' @param depth An integer specifying the sequencing depth for the "Rarefy" and "Rarefy-TSS" methods.
 #' If NULL, no rarefaction is performed.
-#' @param group.var A string representing the group variable's name in the
-#' metadata.
+#' @param subject.var A string specifying the subject variable column in the metadata.
 #' @param time.var A string representing the time variable's name in the
 #' metadata. The default is NULL.
-#' @param subject.var A string specifying the subject variable column in the metadata.
+#' @param group.var A string representing the group variable's name in the
+#' metadata.
 #' @param adj.vars A character vector with the names of adjustment variables in
 #' the metadata.
 #' @param change.base A value indicating the base level for the time variable.
@@ -87,37 +87,37 @@ extract_coef <- function(model) {
 #' generate_alpha_test_pair(
 #' data.obj = peerj32.obj,
 #' alpha.obj = NULL,
-#' time.var = "time",
 #' alpha.name = c("shannon", "simpson", "ace"),
 #' subject.var = "subject",
+#' time.var = "time",
 #' group.var = NULL
 #' )
 #'
 #' generate_alpha_test_pair(
 #' data.obj = peerj32.obj,
 #' alpha.obj = NULL,
-#' time.var = "time",
 #' alpha.name = c("shannon", "simpson", "ace"),
 #' subject.var = "subject",
+#' time.var = "time",
 #' group.var = NULL,
 #' change.base = "2"
 #' )
 #'
 #' generate_alpha_test_pair(
 #' data.obj = peerj32.obj,
 #' alpha.obj = NULL,
-#' time.var = "time",
 #' alpha.name = c("shannon", "simpson", "ace"),
 #' subject.var = "subject",
+#' time.var = "time",
 #' group.var = "group"
 #' )
 #'
 #' generate_alpha_test_pair(
 #' data.obj = peerj32.obj,
 #' alpha.obj = NULL,
-#' time.var = "time",
 #' alpha.name = c("shannon", "simpson", "ace"),
 #' subject.var = "subject",
+#' time.var = "time",
 #' group.var = "group",
 #' adj.vars = "sex"
 #' )
@@ -126,17 +126,17 @@ extract_coef <- function(model) {
 #' generate_alpha_test_pair(
 #' data.obj = subset_pairs.obj,
 #' alpha.obj = NULL,
-#' time.var = "Antibiotic",
 #' alpha.name = c("shannon", "simpson", "ace"),
 #' subject.var = "MouseID",
+#' time.var = "Antibiotic",
 #' group.var = "Sex"
 #' )
 #' generate_alpha_test_pair(
 #' data.obj = subset_pairs.obj,
 #' alpha.obj = NULL,
-#' time.var = "Antibiotic",
 #' alpha.name = c("shannon", "simpson", "ace"),
 #' subject.var = "MouseID",
+#' time.var = "Antibiotic",
 #' group.var = "Sex",
 #' change.base = "Week 2"
 #' )
@@ -146,8 +146,8 @@ generate_alpha_test_pair <-
            alpha.obj = NULL,
            alpha.name = NULL,
            depth = NULL,
-           time.var,
            subject.var,
+           time.var,
            group.var,
            adj.vars = NULL,
            change.base = NULL) {

R/generate_alpha_test_pair.R

@@ -15,12 +15,12 @@
 #' \code{\link[MicrobiomeStat]{mStat_calculate_beta_diversity}} function.
 #' If NULL, beta diversity will be automatically computed from \code{data.obj}
 #' using \code{mStat_calculate_beta_diversity}.
-#' @param time.var The name of the column in metadata containing the time variable.
-#'                This should be a column with time points for each sample. Required
-#'                to identify pairs of samples for the same subject across time.
 #' @param subject.var The name of the column in metadata containing the subject IDs.
 #'                    This should uniquely identify each subject in the study. Required
 #'                    to identify samples that belong to the same subject.
+#' @param time.var The name of the column in metadata containing the time variable.
+#'                This should be a column with time points for each sample. Required
+#'                to identify pairs of samples for the same subject across time.
 #' @param group.var The name of the column in metadata containing the grouping variable
 #'                 to use in linear modeling. This grouping variable will be used as a
 #'                 predictor in the linear models for beta diversity change. Optional.
@@ -43,8 +43,8 @@
 #' generate_beta_change_test_pair(
 #'   data.obj = peerj32.obj,
 #'   dist.obj = NULL,
-#'   time.var = "time",
 #'   subject.var = "subject",
+#'   time.var = "time",
 #'   group.var = "group",
 #'   adj.vars = NULL,
 #'   change.base = "1",
@@ -53,8 +53,8 @@
 #' generate_beta_change_test_pair(
 #'   data.obj = peerj32.obj,
 #'   dist.obj = NULL,
-#'   time.var = "time",
 #'   subject.var = "subject",
+#'   time.var = "time",
 #'   group.var = "group",
 #'   adj.vars = "sex",
 #'   change.base = "1",
@@ -65,8 +65,8 @@
 #' generate_beta_change_test_pair(
 #'   data.obj = subset_pairs.obj,
 #'   dist.obj = NULL,
-#'   time.var = "Antibiotic",
 #'   subject.var = "MouseID",
+#'   time.var = "Antibiotic",
 #'   group.var = "Sex",
 #'   adj.vars = NULL,
 #'   change.base = "Baseline",
@@ -89,8 +89,8 @@
 generate_beta_change_test_pair <-
   function(data.obj,
            dist.obj = NULL,
-           time.var = NULL,
            subject.var,
+           time.var = NULL,
            group.var,
            adj.vars = NULL,
            change.base = NULL,

R/generate_beta_change_test_pair.R

@@ -21,7 +21,6 @@
 #'   \item{Other metadata}{like method, dist.name, etc.}
 #' }
 #' See \code{\link[MicrobiomeStat]{mStat_calculate_PC}} function for details on output format.
-#' @param subject.var String. Variable to be used as subject.
 #' @param time.var String. Variable to be used for time. Default is NULL.
 #' @param t.level Character string specifying the time level/value to subset data to,
 #' if a time variable is provided. Default NULL does not subset data.
@@ -95,7 +94,6 @@
 #'   data.obj = peerj32.obj,
 #'   dist.obj = NULL,
 #'   pc.obj = NULL,
-#'   subject.var = "subject",
 #'   time.var = "time",
 #'   t.level = "2",
 #'   group.var = "group",
@@ -120,7 +118,6 @@
 #'   data.obj = subset_T2D.obj,
 #'   dist.obj = dist.obj,
 #'   pc.obj = pc.obj,
-#'   subject.var = "subject_id",
 #'   time.var = "visit_number_num",
 #'   t.level = NULL,
 #'   group.var = "subject_race",
@@ -141,7 +138,6 @@
 #' @export
 generate_beta_ordination_single <-
   function(data.obj,
-           subject.var,
            time.var = NULL,
            t.level = NULL,
            group.var = NULL,
@@ -174,12 +170,12 @@ generate_beta_ordination_single <-
           # Subset data to specific time point if t.level is provided
           condition <- paste(time.var, "== '", t.level, "'", sep = "")
           data.obj <- mStat_subset_data(data.obj, condition = condition)
-          meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(subject.var,group.var,strata.var,time.var)))
+          meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(group.var,strata.var,time.var)))
           dist.obj <-
             mStat_calculate_beta_diversity(data.obj = data.obj, dist.name = dist.name)
         } else {
           # If no specific time point is provided, use all time points but warn if multiple exist
-          meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(subject.var,group.var,strata.var,time.var)))
+          meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(group.var,strata.var,time.var)))
           if (length(levels(as.factor(meta_tab[,time.var]))) != 1){
             message("Multiple time points detected in your dataset. It is recommended to either set t.level or utilize functions for longitudinal data analysis.")
           }
@@ -188,7 +184,7 @@ generate_beta_ordination_single <-
         }
       } else {
         # If no time variable is provided, proceed with all data
-        meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(subject.var,group.var,strata.var,time.var)))
+        meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(group.var,strata.var,time.var)))
         dist.obj <-
           mStat_calculate_beta_diversity(data.obj = data.obj, dist.name = dist.name)
       }
@@ -203,19 +199,19 @@ generate_beta_ordination_single <-
           if (!is.null(t.level)){
             condition <- paste(time.var, "== '", t.level, "'", sep = "")
             data.obj <- mStat_subset_data(data.obj, condition = condition)
-            meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(subject.var,group.var,strata.var,time.var)))
+            meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(group.var,strata.var,time.var)))
           } else {
-            meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(subject.var,group.var,strata.var,time.var)))
+            meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(group.var,strata.var,time.var)))
             if (length(levels(as.factor(meta_tab[,time.var]))) != 1){
               message("Multiple time points detected in your dataset. It is recommended to either set t.level or utilize functions for longitudinal data analysis.")
             }
           }
         } else {
-          meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(subject.var,group.var,strata.var,time.var)))
+          meta_tab <- data.obj$meta.dat %>% dplyr::select(all_of(c(group.var,strata.var,time.var)))
         }
       }
       if (!is.null(attr(dist.obj[[dist.name[1]]], "labels"))){
-        meta_tab <- attr(dist.obj[[dist.name[1]]], "labels")  %>% dplyr::select(all_of(c(subject.var,group.var,strata.var,time.var)))
+        meta_tab <- attr(dist.obj[[dist.name[1]]], "labels")  %>% dplyr::select(all_of(c(group.var,strata.var,time.var)))
       }
     }
 
@@ -263,7 +259,7 @@ generate_beta_ordination_single <-
         as.data.frame() %>%
         rownames_to_column("sample") %>%
         dplyr::left_join(meta_tab %>%
-                           dplyr::select(all_of(c(subject.var, time.var, group.var, strata.var))) %>%
+                           dplyr::select(all_of(c(time.var, group.var, strata.var))) %>%
                            rownames_to_column("sample"), by = "sample")
 
       # Filter out NA values in time variable if it exists
@@ -362,17 +358,10 @@ generate_beta_ordination_single <-
 
       # Save the plot as a PDF file if requested
       if (pdf) {
-        pdf_name <- paste0(
-          "beta_ordination_single_",
-          "subject_",
-          subject.var,
-          "_",
-          "time_",
-          time.var,
-          "_",
-          "dist.name_",
-          dist.name
-        )
+        pdf_name <- paste0("beta_ordination_single_", "dist.name_", dist.name)
+        if (!is.null(time.var)) {
+          pdf_name <- paste0(pdf_name, "_", "time_", time.var)
+        }
         if (!is.null(group.var)) {
           pdf_name <- paste0(pdf_name, "_", "group_", group.var)
         }

R/generate_beta_ordination_single.R

@@ -163,6 +163,12 @@ generate_beta_test_single <- function(data.obj,
 
   formula <- as.formula(formula_str)
 
+  # Ensure .Random.seed exists (PermanovaG2 requires it)
+  # This initializes the RNG state if it hasn't been used yet
+  if (!exists(".Random.seed", envir = globalenv())) {
+    set.seed(NULL)  # Initialize RNG with system seed
+  }
+
   # Run PermanovaG2 for all distance matrices
   # PERMANOVA is a non-parametric multivariate statistical test used to assess the significance of compositional differences among groups
   message("Running PermanovaG2 for all distances...")

R/generate_beta_test_single.R

@@ -298,9 +298,9 @@ generate_taxa_areaplot_long <-
         otu_tax_agg <- data.obj$feature.tab
       }
 
-      # Subset features if specified
+      # Subset features if specified (using %in% for robustness against NA or non-existent features)
       if (!is.null(features.plot)){
-        otu_tax_agg <- otu_tax_agg[na.omit(features.plot),]
+        otu_tax_agg <- otu_tax_agg[rownames(otu_tax_agg) %in% features.plot,]
       }
 
       # Prepare the feature table for analysis

R/generate_taxa_areaplot_long.R

@@ -317,9 +317,9 @@ generate_taxa_barplot_long <-
         otu_tax_agg <- data.obj$feature.tab
       }
 
-      # Subset features if specified
+      # Subset features if specified (using %in% for robustness against NA or non-existent features)
       if (!is.null(features.plot)){
-        otu_tax_agg <- otu_tax_agg[na.omit(features.plot),]
+        otu_tax_agg <- otu_tax_agg[rownames(otu_tax_agg) %in% features.plot,]
       }
 
       # Prepare the feature table for analysis