long

Author: cafferychen777

.DS_Store

@@ -11,6 +11,8 @@ export(generate_alpha_spaghettiplot_long)
 export(generate_alpha_test_long)
 export(generate_alpha_test_pair)
 export(generate_alpha_test_single)
+export(generate_alpha_trend_test_long)
+export(generate_alpha_volatility_test_long)
 export(generate_beta_change_boxplot_pair)
 export(generate_beta_change_spaghettiplot_long)
 export(generate_beta_change_test_pair)
@@ -19,9 +21,13 @@ export(generate_beta_ordination_pair)
 export(generate_beta_ordination_single)
 export(generate_beta_pc_boxplot_long)
 export(generate_beta_pc_change_boxplot_pair)
+export(generate_beta_pc_trend_test_long)
+export(generate_beta_pc_volatility_test_long)
 export(generate_beta_test_long)
 export(generate_beta_test_pair)
 export(generate_beta_test_single)
+export(generate_beta_trend_test_long)
+export(generate_beta_volatility_test_long)
 export(generate_taxa_areaplot_long)
 export(generate_taxa_barplot_long)
 export(generate_taxa_barplot_pair)
@@ -48,6 +54,8 @@ export(generate_taxa_spaghettiplot_long)
 export(generate_taxa_test_long)
 export(generate_taxa_test_pair)
 export(generate_taxa_test_single)
+export(generate_taxa_trend_test_long)
+export(generate_taxa_volatility_test_long)
 export(is_categorical)
 export(is_continuous_numeric)
 export(is_rarefied)
@@ -56,6 +64,7 @@ export(linda.plot)
 export(mStat_aggregate_by_taxonomy)
 export(mStat_aggregate_data)
 export(mStat_calculate_PC)
+export(mStat_calculate_adjusted_distance)
 export(mStat_calculate_alpha_diversity)
 export(mStat_calculate_beta_diversity)
 export(mStat_combine_data)

NAMESPACE

@@ -13,6 +13,7 @@
 #' @param palette An optional color palette for the plot. If not provided, a default color palette will be used. The palette should be a vector of color codes in a format accepted by ggplot2 (e.g., hexadecimal color codes). The number of colors in the palette should be at least as large as the number of groups being plotted.
 #' @param group.var An optional variable in the metadata table that represents the grouping factor.
 #' @param strata.var An optional variable in the metadata table that represents the stratification factor.
+#' @param adj.vars A character vector of variable names to be used for adjustment.
 #' @param base.size The base font size for the plot.
 #' @param theme.choice A character string indicating the choice of pre-defined ggplot2 theme for the plot. Supported choices include "prism" (default), "classic", "gray", and "bw".
 #' @param custom.theme An optional custom ggplot2 theme. If provided, this theme will be used instead of the pre-defined themes.
@@ -220,7 +221,7 @@ generate_alpha_boxplot_long <- function (data.obj,
 
       # 使用mutate和residuals来添加残差
       alpha_df <- alpha_df %>%
-        mutate(!!sym(index) := residuals(lm(formula_obj, data = alpha_df)))
+        dplyr::mutate(!!sym(index) := stats::residuals(lm(formula_obj, data = alpha_df)))
 
       message("Alpha diversity has been adjusted for the following covariates: ", paste(adj.vars, collapse = ", "), ".")
     }

R/generate_alpha_boxplot_long.R

@@ -10,6 +10,7 @@
 #' @param t.level The specific time point to be analyzed.
 #' @param group.var An optional variable in the metadata table that represents the grouping factor.
 #' @param strata.var An optional variable in the metadata table that represents the stratification factor.
+#' @param adj.vars A character vector of variable names to be used for adjustment.
 #' @param base.size The base font size for the plot. Default is 16.
 #' @param theme.choice A character string indicating the choice of pre-defined ggplot2 theme for the plot. Supported choices include "prism" (default), "classic", "gray", and "bw".
 #' @param custom.theme An optional custom ggplot2 theme. If provided, this theme will be used instead of the pre-defined themes.
@@ -185,7 +186,7 @@ generate_alpha_boxplot_single <- function (data.obj,
 
       # 使用mutate和residuals来添加残差
       alpha_df <- alpha_df %>%
-        mutate(!!sym(index) := residuals(lm(formula_obj, data = alpha_df)))
+        dplyr::mutate(!!sym(index) := stats::residuals(lm(formula_obj, data = alpha_df)))
 
       message("Alpha diversity has been adjusted for the following covariates: ", paste(adj.vars, collapse = ", "), ".")
     }

R/generate_alpha_boxplot_single.R

@@ -15,6 +15,7 @@
 #' @param time.var The variable in the metadata table that represents the time.
 #' @param group.var (Optional) The variable in the metadata table that represents the grouping factor.
 #' @param strata.var (Optional) The variable in the metadata table that represents the stratification factor.
+#' @param adj.vars A character vector of variable names to be used for adjustment.
 #' @param change.base The base time for calculating the change in alpha diversity.
 #' @param change.func (Optional) A function for calculating the change in alpha diversity; can be either "log" (default) or "absolute."
 #' @param pdf (Optional) A boolean indicating whether to save the output as a PDF file.
@@ -239,7 +240,7 @@ generate_alpha_change_boxplot_pair <-
 
         # 使用mutate和residuals来添加残差
         combined_alpha <- combined_alpha %>%
-          mutate(!!sym(paste0(index, "_diff")) := residuals(lm(formula_obj, data = combined_alpha)))
+          dplyr::mutate(!!sym(paste0(index, "_diff")) := stats::residuals(lm(formula_obj, data = combined_alpha)))
 
         message("Alpha diversity Change has been adjusted for the following covariates: ", paste(adj.vars, collapse = ", "), ".")
       }

R/generate_alpha_change_boxplot_pair.R

@@ -18,6 +18,7 @@
 #' @param time.var The name of the time variable.
 #' @param group.var The name of the group variable.
 #' @param strata.var The name of the strata variable (default is NULL).
+#' @param adj.vars A character vector of variable names to be used for adjustment.
 #' @param pdf Logical, whether to save the plot as a PDF (default is TRUE).
 #' @param file.ann The annotation to be added to the PDF file name (default is NULL).
 #' @param ... Additional arguments passed to ggplot().
@@ -182,7 +183,7 @@ generate_alpha_spaghettiplot_long <-
 
         # 使用mutate和residuals来添加残差
         sub_alpha.df <- sub_alpha.df %>%
-          mutate(!!sym(index) := residuals(lm(formula_obj, data = sub_alpha.df)))
+          dplyr::mutate(!!sym(index) := stats::residuals(lm(formula_obj, data = sub_alpha.df)))
 
         message("Alpha diversity has been adjusted for the following covariates: ", paste(adj.vars, collapse = ", "), ".")
       }

R/generate_alpha_spaghettiplot_long.R

@@ -46,6 +46,10 @@ construct_formula <- function(index, group.var, time.var, subject.var, adj.vars)
 #' @param data.obj A list object that includes feature.tab (an OTU table with
 #' taxa as rows and samples as columns) and meta.dat (a metadata table with
 #' samples as rows and variables as columns).
+#' @param alpha.obj An object containing precomputed alpha diversity indices,
+#' typically calculated by the function `mStat_calculate_alpha_diversity`.
+#' If NULL (the default), alpha diversity will be calculated from the data.obj using
+#' `mStat_calculate_alpha_diversity`.
 #' @param alpha.name A string with the name of the alpha diversity index to compute.
 #' Options could include: "shannon", "simpson", "observed_species", "chao1", "ace", and "pielou".
 #' @param time.var A string representing the time variable's name in the
@@ -104,12 +108,14 @@ generate_alpha_trend_test_long <- function(data.obj,
     }
   }
 
-  message("The volatility calculation in generate_alpha_volatility_test_long relies on a numeric time variable.
-         Please check that your time variable is coded as numeric.
-         If the time variable is not numeric, it may cause issues in computing the results of the volatility test.
-         You can ensure the time variable is numeric by mutating it in the metadata.")
+  message(
+    "The trend test in 'generate_alpha_trend_test_long' relies on a numeric time variable.\n",
+    "Please ensure that your time variable is coded as numeric.\n",
+    "If the time variable is not numeric, it may cause issues in computing the results of the trend test.\n",
+    "The time variable will be converted to numeric within the function if needed."
+  )
 
-  data.obj$meta.dat <- data.obj$meta.dat %>% mutate(!!sym(time.var) := as.numeric(!!sym(time.var)))
+  data.obj$meta.dat <- data.obj$meta.dat %>% dplyr::mutate(!!sym(time.var) := as.numeric(!!sym(time.var)))
 
   meta_tab <-
     load_data_obj_metadata(data.obj) %>% as.data.frame() %>% select(all_of(c(
@@ -123,8 +129,6 @@ generate_alpha_trend_test_long <- function(data.obj,
 
   test.list <- lapply(alpha.name, function(index) {
 
-    # Calculate the volatility for each subject
-
     formula <- construct_formula(index, group.var, time.var, subject.var, adj.vars)
 
     model <- lmer(formula, data = alpha_df)

R/generate_alpha_trend_test_long.R

@@ -16,6 +16,10 @@
 #' @param data.obj A list object that includes feature.tab (an OTU table with
 #' taxa as rows and samples as columns) and meta.dat (a metadata table with
 #' samples as rows and variables as columns).
+#' @param alpha.obj An object containing precomputed alpha diversity indices,
+#' typically calculated by the function `mStat_calculate_alpha_diversity`.
+#' If NULL (the default), alpha diversity will be calculated from the data.obj using
+#' `mStat_calculate_alpha_diversity`.
 #' @param alpha.name A string with the name of the alpha diversity index to compute.
 #' Options could include: "shannon", "simpson", "observed_species", "chao1", "ace", and "pielou".
 #' @param time.var A string representing the time variable's name in the
@@ -35,11 +39,12 @@
 #' subset_T2D.obj2$meta.dat$visit_number <- as.numeric(subset_T2D.obj2$meta.dat$visit_number)
 #' generate_alpha_volatility_test_long(
 #' data.obj = subset_T2D.obj2,
+#' alpha.obj = NULL,
 #' alpha.name = c("shannon","simpson"),
 #' time.var = "visit_number",
 #' subject.var = "subject_id",
 #' group.var = "subject_race",
-#' adj.vars = "subject_gender"
+#' adj.vars = "sample_body_site"
 #' )
 #' }
 #' @export
@@ -74,12 +79,14 @@ generate_alpha_volatility_test_long <- function(data.obj,
     }
   }
 
-  message("The volatility calculation in generate_alpha_volatility_test_long relies on a numeric time variable.
-         Please check that your time variable is coded as numeric.
-         If the time variable is not numeric, it may cause issues in computing the results of the volatility test.
-         You can ensure the time variable is numeric by mutating it in the metadata.")
+  message(
+    "The volatility calculation in generate_alpha_volatility_test_long relies on a numeric time variable.\n",
+    "Please check that your time variable is coded as numeric.\n",
+    "If the time variable is not numeric, it may cause issues in computing the results of the volatility test.\n",
+    "You can ensure the time variable is numeric by mutating it in the metadata."
+  )
 
-  data.obj$meta.dat <- data.obj$meta.dat %>% mutate(!!sym(time.var) := as.numeric(!!sym(time.var)))
+  data.obj$meta.dat <- data.obj$meta.dat %>% dplyr::mutate(!!sym(time.var) := as.numeric(!!sym(time.var)))
 
   meta_tab <-
     load_data_obj_metadata(data.obj) %>% as.data.frame() %>% select(all_of(c(
@@ -95,7 +102,9 @@ generate_alpha_volatility_test_long <- function(data.obj,
 
   if (!is.null(adj.vars)){
     # Obtain the residuals by fitting lm(alpha ~ adj.vars)
-    residuals_lm <- lm(as.formula(paste0(index, "~", paste(adj.vars, collapse = "+"))), data = alpha_df)$residuals
+    residuals_lm <- stats::lm(as.formula(paste0(index, "~", paste(adj.vars, collapse = "+"))), data = alpha_df)$residuals
+  } else {
+    residuals_lm <- alpha_df %>% dplyr::select(all_of(c(index))) %>% dplyr::pull()
   }
 
   # Calculate the volatility for each subject
@@ -104,23 +113,24 @@ generate_alpha_volatility_test_long <- function(data.obj,
 
   # Group data by subject and calculate volatility
   volatility_df <- alpha_df %>%
-    arrange(!!sym(subject.var), !!sym(time.var)) %>%
-    group_by(!!sym(subject.var)) %>%
-    mutate(diff_residuals = abs(residuals - lag(residuals)),
-           diff_time = !!sym(time.var) - lag(!!sym(time.var))) %>%
-    filter(!is.na(diff_residuals), !is.na(diff_time)) %>%
-    filter(diff_time != 0) %>%
-    summarize(volatility = mean(diff_residuals / diff_time), .groups = 'drop')
+    dplyr::group_by(!!sym(subject.var)) %>%
+    dplyr::arrange(!!sym(time.var)) %>%
+    dplyr::mutate(diff_residuals = abs(residuals - lag(residuals)),
+           diff_time = !!sym(time.var) - dplyr::lag(!!sym(time.var))) %>%
+    dplyr::filter(!is.na(diff_residuals), !is.na(diff_time)) %>%
+    dplyr::filter(diff_time != 0) %>%
+    dplyr::summarize(volatility = mean(diff_residuals / diff_time), .groups = 'drop')
 
   test_df <- volatility_df %>%
-    left_join(meta_tab %>%
+    dplyr::left_join(meta_tab %>%
                 select(all_of(c(subject.var, group.var))) %>%
-                distinct(),
+                  dplyr::distinct(),
               by = subject.var)
 
   # Test the association between the volatility and the grp.var
   formula <- as.formula(paste("volatility ~", group.var))
-  test_result <- lm(formula, data = test_df)
+
+  test_result <- stats::lm(formula, data = test_df)
 
   coef.tab <- extract_coef(test_result)
 

R/generate_alpha_volatility_test_long.R

@@ -11,6 +11,7 @@
 #' @param time.var A string specifying the name of the time variable
 #' @param group.var A string specifying the name of the group variable or NULL (default)
 #' @param strata.var A string specifying the name of the strata variable or NULL (default)
+#' @param adj.vars A character vector containing the names of the columns in data.obj$meta.dat to include as covariates in the PERMANOVA analysis. If no covariates are needed, use NULL (default).
 #' @param change.base A string or numeric value specifying the first time point used for computing the beta diversity change
 #' @param dist.name A character vector specifying which beta diversity indices to calculate. Supported indices are "BC" (Bray-Curtis), "Jaccard", "UniFrac" (unweighted UniFrac), "GUniFrac" (generalized UniFrac), "WUniFrac" (weighted UniFrac), and "JS" (Jensen-Shannon divergence). If a name is provided but the corresponding object does not exist within dist.obj, it will be computed internally. If the specific index is not supported, an error message will be returned.
 #' @param base.size (Optional) Base font size for the plot (default is 16).
@@ -38,7 +39,7 @@
 #' dist.obj <- mStat_calculate_beta_diversity(peerj32.obj, dist.name = "BC")
 #' generate_beta_change_boxplot_pair(
 #'   data.obj = peerj32.obj,
-#'   dist.obj = dist.obj,
+#'   dist.obj = NULL,
 #'   subject.var = "subject",
 #'   time.var = "time",
 #'   group.var = "group",
@@ -81,6 +82,9 @@ generate_beta_change_boxplot_pair <-
       meta_tab <- load_data_obj_metadata(data.obj) %>% select(all_of(c(subject.var, time.var, group.var, strata.var, adj.vars)))
       dist.obj <-
         mStat_calculate_beta_diversity(data.obj = data.obj, dist.name = dist.name)
+      if (!is.null(adj.vars)){
+        dist.obj <- mStat_calculate_adjusted_distance(data.obj = data.obj, dist.obj = dist.obj, adj.vars = adj.vars, dist.name = dist.name)
+      }
     } else {
       if (!is.null(data.obj) & !is.null(data.obj$meta.dat)){
         meta_tab <- load_data_obj_metadata(data.obj) %>% select(all_of(c(subject.var, time.var, group.var, strata.var, adj.vars)))
@@ -93,10 +97,6 @@ generate_beta_change_boxplot_pair <-
 
     meta_tab <- meta_tab %>% rownames_to_column("sample")
 
-    if (!is.null(adj.vars)){
-      dist.obj <- mStat_calculate_adjusted_distance(data.obj = data.obj, dist.obj = dist.obj, adj.vars = adj.vars, dist.name = dist.name)
-    }
-
     if (is.null(change.base)){
       change.base <- unique(meta_tab %>% select(all_of(c(time.var))))[1,]
       message("The 'change.base' variable was NULL. It has been set to the first unique value in the 'time.var' column of the 'alpha_df' data frame: ", change.base)

R/generate_beta_change_boxplot_pair.R

@@ -9,6 +9,7 @@
 #' @param t0.level (Optional) The baseline time point level.
 #' @param ts.levels (Optional) Time point levels.
 #' @param group.var (Optional) The variable in the metadata table that represents the grouping factor.
+#' @param adj.vars A character vector containing the names of the columns in data.obj$meta.dat to include as covariates in the PERMANOVA analysis. If no covariates are needed, use NULL (default).
 #' @param strata.var (Optional) The variable in the metadata table that represents the stratification factor.
 #' @param dist.name A character vector specifying which beta diversity indices to calculate. Supported indices are "BC" (Bray-Curtis), "Jaccard", "UniFrac" (unweighted UniFrac), "GUniFrac" (generalized UniFrac), "WUniFrac" (weighted UniFrac), and "JS" (Jensen-Shannon divergence). If a name is provided but the corresponding object does not exist within dist.obj, it will be computed internally. If the specific index is not supported, an error message will be returned.
 #' @param base.size (Optional) Base font size for the plot (default is 16).
@@ -104,7 +105,7 @@ generate_beta_change_spaghettiplot_long <-
                              classic = theme_classic(),
                              gray = theme_gray(),
                              bw = theme_bw(),
-                             ggprism::theme_prism()) 
+                             ggprism::theme_prism())
 
     theme_to_use <- if (!is.null(custom.theme)) custom.theme else theme_function
 
@@ -122,6 +123,9 @@ generate_beta_change_spaghettiplot_long <-
         meta_tab <- load_data_obj_metadata(data.obj) %>% select(all_of(c(subject.var, time.var, group.var, strata.var)))
         dist.obj <-
           mStat_calculate_beta_diversity(data.obj = data.obj, dist.name = dist.name)
+        if (!is.null(adj.vars)){
+          dist.obj <- mStat_calculate_adjusted_distance(data.obj = data.obj, dist.obj = dist.obj, adj.vars = adj.vars, dist.name = dist.name)
+        }
       } else {
         if (!is.null(data.obj) & !is.null(data.obj$meta.dat)){
           data.obj <-
@@ -135,10 +139,6 @@ generate_beta_change_spaghettiplot_long <-
         }
       }
 
-      if (!is.null(adj.vars)){
-        dist.obj <- mStat_calculate_adjusted_distance(data.obj = data.obj, dist.obj = dist.obj, adj.vars = adj.vars, dist.name = dist.name)
-      }
-
       if (is.null(dist.obj[[dist.name]])) {
         message(paste("dist.obj does not contain", dist.name))
       } else {