From 234bc474fd4efe9f60574a7ec84d750739e257e2 Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Mon, 26 Aug 2024 09:34:17 +0300
Subject: [PATCH 01/16] short term change

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 R/shortTermChange.R | 82 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 82 insertions(+)
 create mode 100644 R/shortTermChange.R

diff --git a/R/shortTermChange.R b/R/shortTermChange.R
new file mode 100644
index 0000000..c7744a5
--- /dev/null
+++ b/R/shortTermChange.R
@@ -0,0 +1,82 @@
+#' @rdname shortTermChange
+#' @export
+setGeneric("shortTermChange", signature = c("x"),
+    function(x,assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
+        depth = NULL)
+        standardGeneric("shortTermChange"))
+
+#' @rdname shortTermChange
+#' @export
+setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
+    function(x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
+        depth = NULL){
+        ############################## Input check #############################
+        # Check validity of object
+        if(nrow(x) == 0L){
+            stop("No data available in `x` ('x' has nrow(x) == 0L.)",
+                 call. = FALSE)
+        }
+        # Check assay.type
+        .check_assay_present(assay.type, x)
+        if(!.is_a_bool(rarefy)){
+            stop("'rarefy' must be TRUE or FALSE.", call. = FALSE)
+        }
+        if(!.is_a_bool(compositional)){
+            stop("'compositional' must be TRUE or FALSE.", call. = FALSE)
+        }
+        # Ensure that the provided depth is valid
+        if (!is.null(depth)) {
+            if (depth > min(assay(x, assay.type))) {
+                stop("Depth cannot be greater than the minimum number of counts in your data")
+            }
+        } else {
+            stop("Depth cannot be NULL")
+        }
+        
+        ############################ Input check end ###########################
+        
+        ############################ Data Preparation #########################
+        # Initialize the filtered object based on rarefy and compositional arguments
+        if (rarefy == TRUE && compositional == FALSE) {
+            message("rarefy is set to TRUE, calculating short term change using counts")
+            x <- rarefyAssay(x, assay.type = assay.type, depth = depth)
+        } else if (rarefy == FALSE && compositional == FALSE) {
+            message("rarefy is set to FALSE, compositional==FALSE, using raw counts")
+            x <- x
+        } else if (rarefy == FALSE && compositional == TRUE) {
+            message("rarefy is set to FALSE, compositional==TRUE, using relative abundances")
+            x <- transformCounts(x, method = "relative", assay.type = assay.type)
+        } else if (rarefy == TRUE && compositional == TRUE) {
+            stop("Both rarefy and compositional cannot be TRUE simultaneously")
+        }
+        
+        assay.data <- assay(x, assay.type)
+        time <- colData(x)[, "Time.hr"]
+        colnames(assay.data) <- time
+        
+        assay.data <- reshape2::melt(assay.data)
+        
+        ########################### Growth Metrics ############################
+        grwt <- as_tibble(assay.data) %>%
+            arrange(Var2) %>%
+            group_by(Var1) %>%
+            mutate(
+                time_lag = Var2 - lag(Var2), 
+                growth_diff = value - lag(value),
+                growth_rate = value - lag(value) / lag(value),
+                var_abund = value - lag(value) / time_lag
+            )
+        
+        colnames(grwt)[colnames(grwt) == "variable"] <- "OTU"
+        
+        maxgrs <- grwt %>%
+            summarize(max.growth = max(growth_diff, na.rm = T))
+        colnames(maxgrs)[colnames(maxgrs) == "variable"] <- "OTU"
+        grs.all <- merge(gwrt, maxgrs)
+        grs.all <- grs.all %>%
+            mutate(ismax = ifelse(growth_diff == max.growth, T, F))
+
+        grs.all$OTU <- gsub("_", " ", grs.all$OTU)
+    }
+)
+        
\ No newline at end of file

From 4e9568a70c56068f43c29cc4e2b01cb8df315fcc Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Tue, 27 Aug 2024 15:24:12 +0300
Subject: [PATCH 02/16] up

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 R/shortTermChange.R | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/R/shortTermChange.R b/R/shortTermChange.R
index c7744a5..014dfe3 100644
--- a/R/shortTermChange.R
+++ b/R/shortTermChange.R
@@ -67,16 +67,17 @@ setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
                 var_abund = value - lag(value) / time_lag
             )
         
-        colnames(grwt)[colnames(grwt) == "variable"] <- "OTU"
+        colnames(grwt)[colnames(grwt) == "Var1"] <- "OTU"
+        colnames(grwt)[colnames(grwt) == "Var2"] <- "OTU"
         
         maxgrs <- grwt %>%
             summarize(max.growth = max(growth_diff, na.rm = T))
-        colnames(maxgrs)[colnames(maxgrs) == "variable"] <- "OTU"
         grs.all <- merge(gwrt, maxgrs)
         grs.all <- grs.all %>%
             mutate(ismax = ifelse(growth_diff == max.growth, T, F))
 
         grs.all$OTU <- gsub("_", " ", grs.all$OTU)
+        return(grs.all)
     }
 )
         
\ No newline at end of file

From 19937c77995008e9fc78659e8f4b8db366c2ec5d Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Wed, 4 Sep 2024 14:05:16 +0300
Subject: [PATCH 03/16] add plot option

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 NAMESPACE              |   2 +
 R/shortTermChange.R    | 156 ++++++++++++++++++++++++++++++++---------
 man/shortTermChange.Rd |  59 ++++++++++++++++
 3 files changed, 183 insertions(+), 34 deletions(-)
 create mode 100644 man/shortTermChange.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 0dae9a6..95c39ca 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,7 +3,9 @@
 export(getBaselineDivergence)
 export(getStepwiseDivergence)
 export(getTimeDivergence)
+export(shortTermChange)
 exportMethods(getTimeDivergence)
+exportMethods(shortTermChange)
 importFrom(S4Vectors,DataFrame)
 importFrom(SingleCellExperiment,altExp)
 importFrom(SummarizedExperiment,"colData<-")
diff --git a/R/shortTermChange.R b/R/shortTermChange.R
index 014dfe3..0bc758e 100644
--- a/R/shortTermChange.R
+++ b/R/shortTermChange.R
@@ -1,15 +1,57 @@
+
+#' @title Short term Changes in Abundance
+#'
+#' @description Calculates short term changes in abundance of taxa
+#'              using temporal Abundance data.
+#'  
+#' @param x a \code{\link{SummarizedExperiment}} object.
+#' 
+#' @param assay.type \code{Character scalar}. Specifies the name of assay 
+#'   used in calculation. (Default: \code{"counts"})
+#'   
+#' @param rarefy \code{Logical scalar}. Whether to rarefy counts.
+#' (Default: \code{FALSE})
+#' 
+#' @param compositional \code{Logical scalar}. Whether to transform counts.
+#' (Default: \code{FALSE})
+#' 
+#' @param depth \code{Integer scalar}. Specifies the depth used in rarefying. 
+#' (Default: \code{min(assay(x, assay.type)}))
+#' 
+#' @param ... additional arguments.
+#' 
+#' 
+#' @return \code{dataframe} with \code{short term change} calculations.
+#' 
+#' @details This approach is used by Wisnoski NI and colleagues
+#'          \url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
+#'          the following calculation log(present abundance/past abundance).
+#'          Also a compositional version using relative abundance similar to
+#'          Brian WJi, Sheth R et al
+#'          \url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
+#'          This approach is useful for identifying short term growth behaviors of taxa.
+#'          
+#' @name shortTermChange
+#' 
+#' 
+#' @example
+#' data(minimalgut)
+#' tse <- minimalgut
+
+
+
 #' @rdname shortTermChange
 #' @export
 setGeneric("shortTermChange", signature = c("x"),
     function(x,assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
-        depth = NULL)
+        depth = min(assay(x, assay.type)), ...)
         standardGeneric("shortTermChange"))
 
 #' @rdname shortTermChange
 #' @export
 setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
     function(x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
-        depth = NULL){
+        depth = min(assay(x, assay.type)), plot = FALSE, ...){
         ############################## Input check #############################
         # Check validity of object
         if(nrow(x) == 0L){
@@ -17,22 +59,18 @@ setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
                  call. = FALSE)
         }
         # Check assay.type
-        .check_assay_present(assay.type, x)
-        if(!.is_a_bool(rarefy)){
+        mia:::.check_assay_present(assay.type, x)
+        if(!mia:::.is_a_bool(rarefy)){
             stop("'rarefy' must be TRUE or FALSE.", call. = FALSE)
         }
-        if(!.is_a_bool(compositional)){
+        if(!mia:::.is_a_bool(compositional)){
             stop("'compositional' must be TRUE or FALSE.", call. = FALSE)
         }
         # Ensure that the provided depth is valid
-        if (!is.null(depth)) {
-            if (depth > min(assay(x, assay.type))) {
-                stop("Depth cannot be greater than the minimum number of counts in your data")
-            }
-        } else {
-            stop("Depth cannot be NULL")
-        }
+        if ( !is.null(depth) && depth > min(assay(x, assay.type)) ) {
+            stop("Depth cannot be greater than the minimum number of counts in your data")
         
+        }
         ############################ Input check end ###########################
         
         ############################ Data Preparation #########################
@@ -40,44 +78,94 @@ setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
         if (rarefy == TRUE && compositional == FALSE) {
             message("rarefy is set to TRUE, calculating short term change using counts")
             x <- rarefyAssay(x, assay.type = assay.type, depth = depth)
+            assay.type <- "subsampled"
         } else if (rarefy == FALSE && compositional == FALSE) {
             message("rarefy is set to FALSE, compositional==FALSE, using raw counts")
             x <- x
         } else if (rarefy == FALSE && compositional == TRUE) {
             message("rarefy is set to FALSE, compositional==TRUE, using relative abundances")
-            x <- transformCounts(x, method = "relative", assay.type = assay.type)
+            x <- transformAssay(x, method = "relabundance", assay.type = assay.type)
+            assay.type <- "relabundance"
         } else if (rarefy == TRUE && compositional == TRUE) {
             stop("Both rarefy and compositional cannot be TRUE simultaneously")
         }
-        
-        assay.data <- assay(x, assay.type)
-        time <- colData(x)[, "Time.hr"]
-        colnames(assay.data) <- time
-        
-        assay.data <- reshape2::melt(assay.data)
-        
         ########################### Growth Metrics ############################
-        grwt <- as_tibble(assay.data) %>%
-            arrange(Var2) %>%
-            group_by(Var1) %>%
-            mutate(
-                time_lag = Var2 - lag(Var2), 
-                growth_diff = value - lag(value),
-                growth_rate = value - lag(value) / lag(value),
-                var_abund = value - lag(value) / time_lag
-            )
-        
-        colnames(grwt)[colnames(grwt) == "Var1"] <- "OTU"
-        colnames(grwt)[colnames(grwt) == "Var2"] <- "OTU"
-        
+        grwt <- .calculate_growth_metrics(x, assay.type = assay.type, ...)
+        #max growth
         maxgrs <- grwt %>%
             summarize(max.growth = max(growth_diff, na.rm = T))
-        grs.all <- merge(gwrt, maxgrs)
+        grs.all <- merge(grwt, maxgrs, by = "OTU")
         grs.all <- grs.all %>%
             mutate(ismax = ifelse(growth_diff == max.growth, T, F))
 
         grs.all$OTU <- gsub("_", " ", grs.all$OTU)
+        
+        grs.all$OTUabb <- toupper(abbreviate(grs.all$OTU,
+                                             minlength = 3,
+                                             method = "both.sides"
+        ))
+        grs.all$otu.time <- paste0(grs.all$OTUabb, " ", grs.all$time, "h")
+        if(plot){
+            grs.all <- ggplot(grs.all, aes(x = Time, group = OTU, col = OTU)) +
+                geom_line(aes(y = growth_diff), alpha = 0.6, size = 1) +
+                geom_point(
+                    data = subset(grs.all, ismax == T),
+                    aes(y = max.growth), alpha = 0.8, size = 3
+                ) +
+                # coord_flip() +
+                geom_ribbon(aes(group = NULL, col = NULL, ymax = 0, ymin = min(growth_diff)),
+                            fill = "white", col = "#edf3f5", alpha = 0.7
+                ) +
+                theme(
+                    legend.position = "top", legend.text = element_text(size = 9),
+                    panel.background = element_rect(fill = "white", color = "black"),
+                    panel.grid.major = element_blank(),
+                    panel.grid.minor = element_blank(),
+                    legend.key = element_blank()
+                ) +
+                geom_text_repel(
+                    data = subset(grs.all, ismax == T),
+                    aes(y = max.growth, label = otu.time),
+                    nudge_y = 0.2, box.padding = .5, max.iter = 10000,
+                    size = 3, color = "black", segment.alpha = .5,
+                    fontface = "italic", direction = "both"
+                ) +
+                geom_hline(yintercept = 0) +
+                labs(
+                    x = "Time (hr)",
+                    y = expression(paste("Change in abundance ", " \U00B5 = ", 
+                                         abund[t + 1] / abund[t]))
+                )
+        }
         return(grs.all)
+    
     }
 )
+        ########################################################################
+# wrapper to calculate growth matrix
+.calculate_growth_metrics <- function(x, assay.type, ...) {
+    assay_data <- assay(x, assay.type)
+    time <- colData(x)[, "Time.hr"]
+    colnames(assay_data) <- time
+    
+    # Reshape data
+    assay_data <- reshape2::melt(assay_data)
+    assay_data <- as_tibble(assay_data) %>%
+        arrange(Var2) %>%
+        group_by(Var1) %>%
+        mutate(
+            time_lag = Var2 - lag(Var2), 
+            growth_diff = value - lag(value),
+            growth_rate = (value - lag(value)) / lag(value),
+            var_abund = (value - lag(value)) / time_lag
+        )
+    
+    colnames(assay_data)[colnames(assay_data) == "Var1"] <- "OTU"
+    colnames(assay_data)[colnames(assay_data) == "Var2"] <- "Time"
+    
+    return(assay_data)
+}
+
+
+
         
\ No newline at end of file
diff --git a/man/shortTermChange.Rd b/man/shortTermChange.Rd
new file mode 100644
index 0000000..d2bb6b7
--- /dev/null
+++ b/man/shortTermChange.Rd
@@ -0,0 +1,59 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/shortTermChange.R
+\name{shortTermChange}
+\alias{shortTermChange}
+\alias{shortTermChange,SummarizedExperiment-method}
+\title{Short term Changes in Abundance}
+\usage{
+shortTermChange(
+  x,
+  assay.type = "counts",
+  rarefy = FALSE,
+  compositional = FALSE,
+  depth = min(assay(x, assay.type)),
+  ...
+)
+
+\S4method{shortTermChange}{SummarizedExperiment}(
+  x,
+  assay.type = "counts",
+  rarefy = FALSE,
+  compositional = FALSE,
+  depth = min(assay(x, assay.type)),
+  plot = FALSE,
+  ...
+)
+}
+\arguments{
+\item{x}{a \code{\link{SummarizedExperiment}} object.}
+
+\item{assay.type}{\code{Character scalar}. Specifies the name of assay
+used in calculation. (Default: \code{"counts"})}
+
+\item{rarefy}{\code{Logical scalar}. Whether to rarefy counts.
+(Default: \code{FALSE})}
+
+\item{compositional}{\code{Logical scalar}. Whether to transform counts.
+(Default: \code{FALSE})}
+
+\item{depth}{\code{Integer scalar}. Specifies the depth used in rarefying.
+(Default: \code{min(assay(x, assay.type)}))}
+
+\item{...}{additional arguments.}
+}
+\value{
+\code{dataframe} with \code{short term change} calculations.
+}
+\description{
+Calculates short term changes in abundance of taxa
+using temporal Abundance data.
+}
+\details{
+This approach is used by Wisnoski NI and colleagues
+\url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
+the following calculation log(present abundance/past abundance).
+Also a compositional version using relative abundance similar to
+Brian WJi, Sheth R et al
+\url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
+This approach is useful for identifying short term growth behaviors of taxa.
+}

From e451d8659df7fe96cb1cada81816398e7cec6ebc Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Fri, 6 Sep 2024 14:38:03 +0300
Subject: [PATCH 04/16] update fxn

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 DESCRIPTION            |  8 ++++---
 NAMESPACE              |  5 ++++
 R/shortTermChange.R    | 53 ++++++++++++++++++++++++++++--------------
 man/shortTermChange.Rd | 21 ++++++++++++++++-
 4 files changed, 65 insertions(+), 22 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index ad2f372..0af1c5a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -25,12 +25,14 @@ Imports:
     SummarizedExperiment,
     SingleCellExperiment,
     vegan,
-    scater
+    scater,
+    ggplot2,
+    ggrepel,
+    reshape2
 Suggests: 
     TreeSummarizedExperiment,
     tidySingleCellExperiment,
     tidySummarizedExperiment,
-    ggplot2,
     miaViz,
     lubridate,
     rmarkdown,
@@ -42,7 +44,7 @@ Suggests:
 Encoding: UTF-8
 URL: https://github.com/microbiome/miaTime
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr
 LazyData: false
 Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index 95c39ca..393d27c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,11 +12,16 @@ importFrom(SummarizedExperiment,"colData<-")
 importFrom(SummarizedExperiment,assay)
 importFrom(SummarizedExperiment,colData)
 importFrom(dplyr,"%>%")
+importFrom(dplyr,arrange)
+importFrom(dplyr,as_tibble)
 importFrom(dplyr,filter)
 importFrom(dplyr,group_by)
 importFrom(dplyr,mutate)
 importFrom(dplyr,select)
+importFrom(ggplot2,ggplot)
+importFrom(ggrepel,geom_text_repel)
 importFrom(methods,is)
 importFrom(mia,estimateDivergence)
 importFrom(mia,mergeSEs)
+importFrom(mia,rarefyAssay)
 importFrom(vegan,vegdist)
diff --git a/R/shortTermChange.R b/R/shortTermChange.R
index 0bc758e..381f05b 100644
--- a/R/shortTermChange.R
+++ b/R/shortTermChange.R
@@ -1,4 +1,3 @@
-
 #' @title Short term Changes in Abundance
 #'
 #' @description Calculates short term changes in abundance of taxa
@@ -18,10 +17,14 @@
 #' @param depth \code{Integer scalar}. Specifies the depth used in rarefying. 
 #' (Default: \code{min(assay(x, assay.type)}))
 #' 
+#' @param plot \code{Logical scalar}. Whether to plot short term change.
+#' (Default: \code{FALSE})
+#' 
 #' @param ... additional arguments.
 #' 
 #' 
-#' @return \code{dataframe} with \code{short term change} calculations.
+#' @return \code{dataframe} or \code{plot} with \code{short term change} 
+#' calculations.
 #' 
 #' @details This approach is used by Wisnoski NI and colleagues
 #'          \url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
@@ -34,10 +37,20 @@
 #' @name shortTermChange
 #' 
 #' 
-#' @example
+#' @examples
+#' 
+#' # Load time series data
 #' data(minimalgut)
 #' tse <- minimalgut
-
+#' 
+#' short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
+#' 
+#' # Subset samples by Time_lable and StudyIdentifier
+#' tse <- tse[, !(colData(tse)$Time_label %in% short_time_labels)]
+#' tse <- tse[, (colData(tse)$StudyIdentifier == "Bioreactor A")]
+#' 
+#' # Plot short term change in abundance
+#' shortTermChange(tse, rarefy = TRUE, plot = TRUE) + ggtitle("Bioreactor A")
 
 
 #' @rdname shortTermChange
@@ -49,6 +62,10 @@ setGeneric("shortTermChange", signature = c("x"),
 
 #' @rdname shortTermChange
 #' @export
+#' @importFrom dplyr arrange as_tibble
+#' @importFrom ggplot2 ggplot
+#' @importFrom ggrepel geom_text_repel
+#' @importFrom mia rarefyAssay
 setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
     function(x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
         depth = min(assay(x, assay.type)), plot = FALSE, ...){
@@ -112,15 +129,20 @@ setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
                     data = subset(grs.all, ismax == T),
                     aes(y = max.growth), alpha = 0.8, size = 3
                 ) +
-                # coord_flip() +
-                geom_ribbon(aes(group = NULL, col = NULL, ymax = 0, ymin = min(growth_diff)),
-                            fill = "white", col = "#edf3f5", alpha = 0.7
+                geom_ribbon(aes(group = NULL, col = NULL, ymax = 0, ymin = -9),
+                            fill = "#edf3f5", col = "#edf3f5", alpha = 0.5
                 ) +
                 theme(
-                    legend.position = "top", legend.text = element_text(size = 9),
-                    panel.background = element_rect(fill = "white", color = "black"),
-                    panel.grid.major = element_blank(),
-                    panel.grid.minor = element_blank(),
+                    legend.position = "right", legend.text = element_text(size = 9),
+                    panel.background = element_rect(fill = "white"),
+                    panel.grid.major = element_line(
+                        size = 0.5, linetype = "solid",
+                        colour = "#CCD1D1"
+                    ),
+                    panel.grid.minor = element_line(
+                        size = 0.5, linetype = "solid",
+                        colour = "#CCD1D1"
+                    ),
                     legend.key = element_blank()
                 ) +
                 geom_text_repel(
@@ -133,8 +155,7 @@ setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
                 geom_hline(yintercept = 0) +
                 labs(
                     x = "Time (hr)",
-                    y = expression(paste("Change in abundance ", " \U00B5 = ", 
-                                         abund[t + 1] / abund[t]))
+                    y = expression(paste("Change in abundance ", " \U00B5 = ", abund[t + 1] / abund[t]))
                 )
         }
         return(grs.all)
@@ -164,8 +185,4 @@ setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
     colnames(assay_data)[colnames(assay_data) == "Var2"] <- "Time"
     
     return(assay_data)
-}
-
-
-
-        
\ No newline at end of file
+}
\ No newline at end of file
diff --git a/man/shortTermChange.Rd b/man/shortTermChange.Rd
index d2bb6b7..0a2c002 100644
--- a/man/shortTermChange.Rd
+++ b/man/shortTermChange.Rd
@@ -40,9 +40,13 @@ used in calculation. (Default: \code{"counts"})}
 (Default: \code{min(assay(x, assay.type)}))}
 
 \item{...}{additional arguments.}
+
+\item{plot}{\code{Logical scalar}. Whether to plot short term change.
+(Default: \code{FALSE})}
 }
 \value{
-\code{dataframe} with \code{short term change} calculations.
+\code{dataframe} or \code{plot} with \code{short term change}
+calculations.
 }
 \description{
 Calculates short term changes in abundance of taxa
@@ -57,3 +61,18 @@ Brian WJi, Sheth R et al
 \url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
 This approach is useful for identifying short term growth behaviors of taxa.
 }
+\examples{
+
+# Load time series data
+data(minimalgut)
+tse <- minimalgut
+
+short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
+
+# Subset samples by Time_lable and StudyIdentifier
+tse <- tse[, !(colData(tse)$Time_label \%in\% short_time_labels)]
+tse <- tse[, (colData(tse)$StudyIdentifier == "Bioreactor A")]
+
+# Plot short term change in abundance
+shortTermChange(tse, rarefy = TRUE, plot = TRUE) + ggtitle("Bioreactor A")
+}

From 07e33b488016d725e65ffd24500b204e24202afe Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Wed, 9 Oct 2024 09:23:19 +0300
Subject: [PATCH 05/16] tests for short term change

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 NAMESPACE                             |  2 +
 R/shortTermChange.R                   |  4 +-
 R/utils.R                             |  1 -
 tests/testthat/test-shortTermChange.R | 80 +++++++++++++++++++++++++++
 4 files changed, 84 insertions(+), 3 deletions(-)
 create mode 100644 tests/testthat/test-shortTermChange.R

diff --git a/NAMESPACE b/NAMESPACE
index 393d27c..dfc1273 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -18,6 +18,8 @@ importFrom(dplyr,filter)
 importFrom(dplyr,group_by)
 importFrom(dplyr,mutate)
 importFrom(dplyr,select)
+importFrom(dplyr,summarize)
+importFrom(ggplot2,aes)
 importFrom(ggplot2,ggplot)
 importFrom(ggrepel,geom_text_repel)
 importFrom(methods,is)
diff --git a/R/shortTermChange.R b/R/shortTermChange.R
index 381f05b..41a695e 100644
--- a/R/shortTermChange.R
+++ b/R/shortTermChange.R
@@ -62,8 +62,8 @@ setGeneric("shortTermChange", signature = c("x"),
 
 #' @rdname shortTermChange
 #' @export
-#' @importFrom dplyr arrange as_tibble
-#' @importFrom ggplot2 ggplot
+#' @importFrom dplyr arrange as_tibble summarize
+#' @importFrom ggplot2 ggplot aes
 #' @importFrom ggrepel geom_text_repel
 #' @importFrom mia rarefyAssay
 setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
diff --git a/R/utils.R b/R/utils.R
index 2b5fdf7..f3a0077 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -2,7 +2,6 @@
 # internal methods loaded from other packages
 
 .check_altExp_present <- mia:::.check_altExp_present
-.calc_reference_dist <- mia:::.calc_reference_dist
 .get_mat_from_sce <- scater:::.get_mat_from_sce
 
 ################################################################################
diff --git a/tests/testthat/test-shortTermChange.R b/tests/testthat/test-shortTermChange.R
new file mode 100644
index 0000000..a76ee92
--- /dev/null
+++ b/tests/testthat/test-shortTermChange.R
@@ -0,0 +1,80 @@
+test_that("shortTermChange", {
+    library(SummarizedExperiment)
+    # Load dataset
+    data(minimalgut)
+    tse <- minimalgut
+    # Check if the function handles empty input
+    empty_se <- SummarizedExperiment()
+    expect_error(shortTermChange(empty_se), 
+                 "No data available in `x`")
+    # Check if assay.type argument works
+    # tse_invalid <- tse
+    # expect_error(
+    #     shortTermChange(tse_invalid, assay.type = "invalid_assay"),
+    #     "'assay.type' must be a valid name of assays(x)"
+    # )
+    # Check that rarefy and compositional cannot both be TRUE
+    expect_error(shortTermChange(tse, rarefy = TRUE, compositional = TRUE), 
+                 "Both rarefy and compositional cannot be TRUE simultaneously")
+    # Check if the depth argument is greater than minimum counts
+    min_depth <- min(assay(tse, "counts"))
+    expect_error(shortTermChange(tse, depth = min_depth + 1),
+                 "Depth cannot be greater than the minimum number of counts in your data")
+    # Check if rarefy = TRUE works
+    result <- shortTermChange(tse, rarefy = TRUE)
+    expect_true(is.data.frame(result))
+    # Check if compositional = TRUE works
+    result <- shortTermChange(tse, compositional = TRUE)
+    expect_true(is.data.frame(result))  
+    # Expected output should be a ggplot object
+    result <- shortTermChange(tse, plot = TRUE)
+    expect_s3_class(result, "gg")  
+    # Should still return a dataframe
+    result <- shortTermChange(tse, rarefy = TRUE, 
+                              compositional = FALSE, 
+                              additional_arg = "value")
+    expect_true(is.data.frame(result))  
+    
+    short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
+    # Subset samples by Time_label and StudyIdentifier 
+    tse_filtered <- tse[, !(colData(tse)$Time_label %in% short_time_labels)]
+    tse_filtered <- tse_filtered[, (colData(tse_filtered)$StudyIdentifier == "Bioreactor A")]
+    
+    expect_true(all(!(colData(tse_filtered)$Time_label %in% short_time_labels)))
+    # Expected output should be a ggplot object
+    result <- shortTermChange(tse_filtered, 
+                              rarefy = TRUE, plot = TRUE)
+    expect_s3_class(result, "gg")  
+    
+    result <- shortTermChange(tse_filtered, 
+                              rarefy = FALSE, plot = FALSE)
+    # Expected output is a dataframe
+    expect_true(is.data.frame(result))  
+    # Ensure dataframe has expected columns
+    expect_true("OTU" %in% colnames(result))  
+    expect_true("growth_diff" %in% colnames(result))
+    
+    result <- shortTermChange(tse_filtered, 
+                              rarefy = TRUE, plot = FALSE)
+    
+    result <- shortTermChange(tse_filtered, 
+                              compositional = TRUE, plot = FALSE)
+    # Expected output is a dataframe
+    expect_true(is.data.frame(result))  
+    
+    result <- shortTermChange(tse_filtered, rarefy = FALSE, 
+                              compositional = FALSE, 
+                              plot = FALSE)
+    expect_true("growth_diff" %in% colnames(result))
+    # Test some expected properties (e.g., that growth_diff isn't all NAs)
+    expect_false(all(is.na(result$growth_diff)))
+    
+    min_depth <- min(assay(tse_filtered, "counts"))
+    result <- shortTermChange(tse_filtered, rarefy = TRUE, depth = min_depth)
+    expect_true(is.data.frame(result)) 
+    expect_error(shortTermChange(tse_filtered, 
+                                 rarefy = TRUE, 
+                                 depth = min_depth + 1),
+                 "Depth cannot be greater than the minimum number of counts in your data")
+    
+})
\ No newline at end of file

From 452a709867af9f5aa08a030ba3a0cf3c66d9b28a Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Thu, 10 Oct 2024 11:46:58 +0300
Subject: [PATCH 06/16] import functions

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 NAMESPACE           | 1 +
 R/shortTermChange.R | 3 ++-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/NAMESPACE b/NAMESPACE
index dfc1273..51335d8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,4 +26,5 @@ importFrom(methods,is)
 importFrom(mia,estimateDivergence)
 importFrom(mia,mergeSEs)
 importFrom(mia,rarefyAssay)
+importFrom(mia,transformAssay)
 importFrom(vegan,vegdist)
diff --git a/R/shortTermChange.R b/R/shortTermChange.R
index 41a695e..b7967c0 100644
--- a/R/shortTermChange.R
+++ b/R/shortTermChange.R
@@ -65,7 +65,8 @@ setGeneric("shortTermChange", signature = c("x"),
 #' @importFrom dplyr arrange as_tibble summarize
 #' @importFrom ggplot2 ggplot aes
 #' @importFrom ggrepel geom_text_repel
-#' @importFrom mia rarefyAssay
+#' @importFrom mia rarefyAssay transformAssay
+#' @importFrom SummarizedExperiment colData
 setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
     function(x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
         depth = min(assay(x, assay.type)), plot = FALSE, ...){

From 6755ca97dd8d596febcb84336214eb72b0e9c4ee Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Thu, 10 Oct 2024 12:20:01 +0300
Subject: [PATCH 07/16] update dependencies

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 DESCRIPTION | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 0af1c5a..7f56d4b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -16,17 +16,17 @@ Description:
 biocViews: Microbiome, Software, Sequencing, Coverage
 License: Artistic-2.0 | file LICENSE
 Depends:
-    R (>= 4.0)
+    R (>= 4.0),
+    ggplot2,
+    mia,
+    SummarizedExperiment,
 Imports:
     dplyr,
     methods,    
-    mia,
     S4Vectors,
-    SummarizedExperiment,
     SingleCellExperiment,
     vegan,
     scater,
-    ggplot2,
     ggrepel,
     reshape2
 Suggests: 

From d0ac31fbdc85800f86e783f0f429eccd227f8c7e Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Fri, 8 Nov 2024 13:43:18 +0200
Subject: [PATCH 08/16] update short term change

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 DESCRIPTION                                   |   3 +-
 NAMESPACE                                     |  10 +
 R/getShortTermChange.R                        | 146 ++++++++++++++
 R/shortTermChange.R                           | 189 ------------------
 ...ortTermChange.Rd => getShortTermChange.Rd} |  26 +--
 man/miaTime-package.Rd                        |   4 +-
 tests/testthat/test-getShortTermChange.R      |  64 ++++++
 tests/testthat/test-shortTermChange.R         |  80 --------
 8 files changed, 234 insertions(+), 288 deletions(-)
 create mode 100644 R/getShortTermChange.R
 delete mode 100644 R/shortTermChange.R
 rename man/{shortTermChange.Rd => getShortTermChange.Rd} (71%)
 create mode 100644 tests/testthat/test-getShortTermChange.R
 delete mode 100644 tests/testthat/test-shortTermChange.R

diff --git a/DESCRIPTION b/DESCRIPTION
index e502f7f..b941c2e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -25,8 +25,7 @@ Depends:
     R (>= 4.4.0),
     ggplot2,
     mia,
-    SummarizedExperiment,,
-    mia
+    SummarizedExperiment
 Imports:
     dplyr,
     methods,
diff --git a/NAMESPACE b/NAMESPACE
index 1ec2e56..0060ed1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,15 +3,25 @@
 export(addBaselineDivergence)
 export(addStepwiseDivergence)
 export(getBaselineDivergence)
+export(getShortTermChange)
 export(getStepwiseDivergence)
 exportMethods(addBaselineDivergence)
 exportMethods(addStepwiseDivergence)
 exportMethods(getBaselineDivergence)
+exportMethods(getShortTermChange)
 exportMethods(getStepwiseDivergence)
 import(TreeSummarizedExperiment)
 import(mia)
+importFrom(SummarizedExperiment,colData)
 importFrom(dplyr,arrange)
+importFrom(dplyr,as_tibble)
 importFrom(dplyr,group_by)
 importFrom(dplyr,lag)
 importFrom(dplyr,mutate)
+importFrom(dplyr,summarize)
 importFrom(dplyr,ungroup)
+importFrom(ggplot2,aes)
+importFrom(ggplot2,ggplot)
+importFrom(ggrepel,geom_text_repel)
+importFrom(mia,rarefyAssay)
+importFrom(mia,transformAssay)
diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
new file mode 100644
index 0000000..1aacccd
--- /dev/null
+++ b/R/getShortTermChange.R
@@ -0,0 +1,146 @@
+#' @title Short term Changes in Abundance
+#'
+#' @description Calculates short term changes in abundance of taxa
+#'              using temporal Abundance data.
+#'  
+#' @param x a \code{\link{SummarizedExperiment}} object.
+#' @param assay.type \code{Character scalar}. Specifies the name of assay 
+#'   used in calculation. (Default: \code{"counts"})
+#' @param rarefy \code{Logical scalar}. Whether to rarefy counts.
+#' (Default: \code{FALSE})
+#' @param compositional \code{Logical scalar}. Whether to transform counts.
+#' (Default: \code{FALSE})
+#' @param depth \code{Integer scalar}. Specifies the depth used in rarefying. 
+#' (Default: \code{min(assay(x, assay.type)}))
+#' @param ... additional arguments.
+#' 
+#' 
+#' @return \code{dataframe} with \code{short term change} 
+#' calculations.
+#' 
+#' @details This approach is used by Wisnoski NI and colleagues
+#'          \url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
+#'          the following calculation log(present abundance/past abundance).
+#'          Also a compositional version using relative abundance similar to
+#'          Brian WJi, Sheth R et al
+#'          \url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
+#'          This approach is useful for identifying short term growth behaviors of taxa.
+#'          
+#' @name getShortTermChange
+#' 
+#' 
+#' @examples
+#' 
+#' # Load time series data
+#' data(minimalgut)
+#' tse <- minimalgut
+#' 
+#' short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
+#' 
+#' # Subset samples by Time_lable and StudyIdentifier
+#' tse <- tse[, !(tse$Time_label %in% short_time_labels)]
+#' tse <- tse[, (tse$StudyIdentifier == "Bioreactor A")]
+#' 
+#' # Get short term change
+#' getShortTermChange(tse, rarefy = TRUE, time.col = "Time.hr")
+NULL
+
+#' @rdname getShortTermChange
+#' @export
+setGeneric("getShortTermChange", signature = c("x"),
+    function( x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
+        depth = min(assay(x, assay.type)), ...)
+        standardGeneric("getShortTermChange"))
+
+#' @rdname getShortTermChange
+#' @export
+#' @importFrom dplyr arrange as_tibble summarize
+#' @importFrom ggplot2 ggplot aes
+#' @importFrom ggrepel geom_text_repel
+#' @importFrom mia rarefyAssay transformAssay
+#' @importFrom SummarizedExperiment colData
+setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
+    function(x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
+        depth = min(assay(x, assay.type)), ...){
+        ############################## Input check #############################
+        # Check validity of object
+        if(nrow(x) == 0L){
+            stop("No data available in `x` ('x' has nrow(x) == 0L.)",
+                 call. = FALSE)
+        }
+        # Check assay.type
+        .check_assay_present(assay.type, x)
+        if(!.is_a_bool(rarefy)){
+            stop("'rarefy' must be TRUE or FALSE.", call. = FALSE)
+        }
+        if(!.is_a_bool(compositional)){
+            stop("'compositional' must be TRUE or FALSE.", call. = FALSE)
+        }
+        # Ensure that the provided depth is valid
+        if ( !is.null(depth) && depth > min(assay(x, assay.type), na.rm = TRUE) ) {
+            stop("Depth cannot be greater than the minimum number of counts in your data", call. = FALSE)
+        
+        }
+        ########################### Growth Metrics ############################
+        grwt <- .calculate_growth_metrics(x, assay.type = assay.type, 
+                                          rarefy = rarefy, 
+                                          compositional = compositional, 
+                                          depth = depth, ...)
+        # Clean and format growth matrics
+        grs.all <- .clean_growth_metrics(grwt, ...)
+        return(grs.all)
+    }
+)
+# wrapper to calculate growth matrix
+.calculate_growth_metrics <- function(x, assay.type, time.col = NULL, 
+                                      rarefy, compositional, depth, ...) {
+    ############################ Data Preparation ##############################
+    # Initialize the filtered object based on rarefy and compositional arguments
+    if (rarefy == TRUE && compositional == FALSE) {
+        message("rarefy is set to TRUE, calculating short term change using counts")
+        x <- rarefyAssay(x, assay.type = assay.type, depth = depth)
+        assay.type <- "subsampled"
+    } else if (rarefy == FALSE && compositional == FALSE) {
+        message("rarefy is set to FALSE, compositional==FALSE, using raw counts")
+        x <- x
+    } else if (rarefy == FALSE && compositional == TRUE) {
+        message("rarefy is set to FALSE, compositional==TRUE, using relative abundances")
+        x <- transformAssay(x, method = "relabundance", assay.type = assay.type)
+        assay.type <- "relabundance"
+    } else if (rarefy == TRUE && compositional == TRUE) {
+        stop("Both rarefy and compositional cannot be TRUE simultaneously", call. = FALSE)
+    }
+    # Reshape data and calcualte grwoth metrics
+    assay_data <- meltSE(x, assay.type = assay.type, add.col = time.col)
+    assay_data <- assay_data %>%
+        arrange( !!sym(time.col) ) %>%
+        group_by(FeatureID) %>%
+        mutate(
+            time_lag = !!sym(time.col) - lag( !!sym(time.col) ), 
+            growth_diff =!!sym(assay.type) - lag(!!sym(assay.type)),
+            growth_rate = (!!sym(assay.type) - lag(!!sym(assay.type))) / lag(!!sym(assay.type)),
+            var_abund = (!!sym(assay.type) - lag(!!sym(assay.type))) / time_lag
+        )
+    return(assay_data)
+}
+
+.clean_growth_metrics <- function(grwt, time.col = NULL, ...) {
+    # Calculate max growth
+    maxgrs <- grwt %>%
+        summarize(max.growth = max(growth_diff, na.rm = TRUE))
+    # Merge growth data with max growth
+    grs.all <- merge(grwt, maxgrs, by = "FeatureID")
+    # Add 'ismax' column indicating if the growth is the maximum
+    grs.all <- grs.all %>%
+        mutate(ismax = ifelse(growth_diff == max.growth, TRUE, FALSE))
+    # Clean and abbreviate FeatureID names
+    grs.all$FeatureID <- gsub("_", " ", grs.all$FeatureID)
+    grs.all$FeatureIDabb <- toupper(abbreviate(grs.all$FeatureID, 
+                                               minlength = 3, 
+                                               method = "both.sides"))
+    # Create 'Feature.time' column combining abbreviation and time information
+    grs.all$Feature.time <- paste0(grs.all$FeatureIDabb, " ", 
+                                   grs.all[[time.col]], "h")
+    
+    return(grs.all)
+}
diff --git a/R/shortTermChange.R b/R/shortTermChange.R
deleted file mode 100644
index b7967c0..0000000
--- a/R/shortTermChange.R
+++ /dev/null
@@ -1,189 +0,0 @@
-#' @title Short term Changes in Abundance
-#'
-#' @description Calculates short term changes in abundance of taxa
-#'              using temporal Abundance data.
-#'  
-#' @param x a \code{\link{SummarizedExperiment}} object.
-#' 
-#' @param assay.type \code{Character scalar}. Specifies the name of assay 
-#'   used in calculation. (Default: \code{"counts"})
-#'   
-#' @param rarefy \code{Logical scalar}. Whether to rarefy counts.
-#' (Default: \code{FALSE})
-#' 
-#' @param compositional \code{Logical scalar}. Whether to transform counts.
-#' (Default: \code{FALSE})
-#' 
-#' @param depth \code{Integer scalar}. Specifies the depth used in rarefying. 
-#' (Default: \code{min(assay(x, assay.type)}))
-#' 
-#' @param plot \code{Logical scalar}. Whether to plot short term change.
-#' (Default: \code{FALSE})
-#' 
-#' @param ... additional arguments.
-#' 
-#' 
-#' @return \code{dataframe} or \code{plot} with \code{short term change} 
-#' calculations.
-#' 
-#' @details This approach is used by Wisnoski NI and colleagues
-#'          \url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
-#'          the following calculation log(present abundance/past abundance).
-#'          Also a compositional version using relative abundance similar to
-#'          Brian WJi, Sheth R et al
-#'          \url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
-#'          This approach is useful for identifying short term growth behaviors of taxa.
-#'          
-#' @name shortTermChange
-#' 
-#' 
-#' @examples
-#' 
-#' # Load time series data
-#' data(minimalgut)
-#' tse <- minimalgut
-#' 
-#' short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
-#' 
-#' # Subset samples by Time_lable and StudyIdentifier
-#' tse <- tse[, !(colData(tse)$Time_label %in% short_time_labels)]
-#' tse <- tse[, (colData(tse)$StudyIdentifier == "Bioreactor A")]
-#' 
-#' # Plot short term change in abundance
-#' shortTermChange(tse, rarefy = TRUE, plot = TRUE) + ggtitle("Bioreactor A")
-
-
-#' @rdname shortTermChange
-#' @export
-setGeneric("shortTermChange", signature = c("x"),
-    function(x,assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
-        depth = min(assay(x, assay.type)), ...)
-        standardGeneric("shortTermChange"))
-
-#' @rdname shortTermChange
-#' @export
-#' @importFrom dplyr arrange as_tibble summarize
-#' @importFrom ggplot2 ggplot aes
-#' @importFrom ggrepel geom_text_repel
-#' @importFrom mia rarefyAssay transformAssay
-#' @importFrom SummarizedExperiment colData
-setMethod("shortTermChange", signature = c(x = "SummarizedExperiment"),
-    function(x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
-        depth = min(assay(x, assay.type)), plot = FALSE, ...){
-        ############################## Input check #############################
-        # Check validity of object
-        if(nrow(x) == 0L){
-            stop("No data available in `x` ('x' has nrow(x) == 0L.)",
-                 call. = FALSE)
-        }
-        # Check assay.type
-        mia:::.check_assay_present(assay.type, x)
-        if(!mia:::.is_a_bool(rarefy)){
-            stop("'rarefy' must be TRUE or FALSE.", call. = FALSE)
-        }
-        if(!mia:::.is_a_bool(compositional)){
-            stop("'compositional' must be TRUE or FALSE.", call. = FALSE)
-        }
-        # Ensure that the provided depth is valid
-        if ( !is.null(depth) && depth > min(assay(x, assay.type)) ) {
-            stop("Depth cannot be greater than the minimum number of counts in your data")
-        
-        }
-        ############################ Input check end ###########################
-        
-        ############################ Data Preparation #########################
-        # Initialize the filtered object based on rarefy and compositional arguments
-        if (rarefy == TRUE && compositional == FALSE) {
-            message("rarefy is set to TRUE, calculating short term change using counts")
-            x <- rarefyAssay(x, assay.type = assay.type, depth = depth)
-            assay.type <- "subsampled"
-        } else if (rarefy == FALSE && compositional == FALSE) {
-            message("rarefy is set to FALSE, compositional==FALSE, using raw counts")
-            x <- x
-        } else if (rarefy == FALSE && compositional == TRUE) {
-            message("rarefy is set to FALSE, compositional==TRUE, using relative abundances")
-            x <- transformAssay(x, method = "relabundance", assay.type = assay.type)
-            assay.type <- "relabundance"
-        } else if (rarefy == TRUE && compositional == TRUE) {
-            stop("Both rarefy and compositional cannot be TRUE simultaneously")
-        }
-        ########################### Growth Metrics ############################
-        grwt <- .calculate_growth_metrics(x, assay.type = assay.type, ...)
-        #max growth
-        maxgrs <- grwt %>%
-            summarize(max.growth = max(growth_diff, na.rm = T))
-        grs.all <- merge(grwt, maxgrs, by = "OTU")
-        grs.all <- grs.all %>%
-            mutate(ismax = ifelse(growth_diff == max.growth, T, F))
-
-        grs.all$OTU <- gsub("_", " ", grs.all$OTU)
-        
-        grs.all$OTUabb <- toupper(abbreviate(grs.all$OTU,
-                                             minlength = 3,
-                                             method = "both.sides"
-        ))
-        grs.all$otu.time <- paste0(grs.all$OTUabb, " ", grs.all$time, "h")
-        if(plot){
-            grs.all <- ggplot(grs.all, aes(x = Time, group = OTU, col = OTU)) +
-                geom_line(aes(y = growth_diff), alpha = 0.6, size = 1) +
-                geom_point(
-                    data = subset(grs.all, ismax == T),
-                    aes(y = max.growth), alpha = 0.8, size = 3
-                ) +
-                geom_ribbon(aes(group = NULL, col = NULL, ymax = 0, ymin = -9),
-                            fill = "#edf3f5", col = "#edf3f5", alpha = 0.5
-                ) +
-                theme(
-                    legend.position = "right", legend.text = element_text(size = 9),
-                    panel.background = element_rect(fill = "white"),
-                    panel.grid.major = element_line(
-                        size = 0.5, linetype = "solid",
-                        colour = "#CCD1D1"
-                    ),
-                    panel.grid.minor = element_line(
-                        size = 0.5, linetype = "solid",
-                        colour = "#CCD1D1"
-                    ),
-                    legend.key = element_blank()
-                ) +
-                geom_text_repel(
-                    data = subset(grs.all, ismax == T),
-                    aes(y = max.growth, label = otu.time),
-                    nudge_y = 0.2, box.padding = .5, max.iter = 10000,
-                    size = 3, color = "black", segment.alpha = .5,
-                    fontface = "italic", direction = "both"
-                ) +
-                geom_hline(yintercept = 0) +
-                labs(
-                    x = "Time (hr)",
-                    y = expression(paste("Change in abundance ", " \U00B5 = ", abund[t + 1] / abund[t]))
-                )
-        }
-        return(grs.all)
-    
-    }
-)
-        ########################################################################
-# wrapper to calculate growth matrix
-.calculate_growth_metrics <- function(x, assay.type, ...) {
-    assay_data <- assay(x, assay.type)
-    time <- colData(x)[, "Time.hr"]
-    colnames(assay_data) <- time
-    
-    # Reshape data
-    assay_data <- reshape2::melt(assay_data)
-    assay_data <- as_tibble(assay_data) %>%
-        arrange(Var2) %>%
-        group_by(Var1) %>%
-        mutate(
-            time_lag = Var2 - lag(Var2), 
-            growth_diff = value - lag(value),
-            growth_rate = (value - lag(value)) / lag(value),
-            var_abund = (value - lag(value)) / time_lag
-        )
-    
-    colnames(assay_data)[colnames(assay_data) == "Var1"] <- "OTU"
-    colnames(assay_data)[colnames(assay_data) == "Var2"] <- "Time"
-    
-    return(assay_data)
-}
\ No newline at end of file
diff --git a/man/shortTermChange.Rd b/man/getShortTermChange.Rd
similarity index 71%
rename from man/shortTermChange.Rd
rename to man/getShortTermChange.Rd
index 0a2c002..a7e1e23 100644
--- a/man/shortTermChange.Rd
+++ b/man/getShortTermChange.Rd
@@ -1,11 +1,11 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/shortTermChange.R
-\name{shortTermChange}
-\alias{shortTermChange}
-\alias{shortTermChange,SummarizedExperiment-method}
+% Please edit documentation in R/getShortTermChange.R
+\name{getShortTermChange}
+\alias{getShortTermChange}
+\alias{getShortTermChange,SummarizedExperiment-method}
 \title{Short term Changes in Abundance}
 \usage{
-shortTermChange(
+getShortTermChange(
   x,
   assay.type = "counts",
   rarefy = FALSE,
@@ -14,13 +14,12 @@ shortTermChange(
   ...
 )
 
-\S4method{shortTermChange}{SummarizedExperiment}(
+\S4method{getShortTermChange}{SummarizedExperiment}(
   x,
   assay.type = "counts",
   rarefy = FALSE,
   compositional = FALSE,
   depth = min(assay(x, assay.type)),
-  plot = FALSE,
   ...
 )
 }
@@ -40,12 +39,9 @@ used in calculation. (Default: \code{"counts"})}
 (Default: \code{min(assay(x, assay.type)}))}
 
 \item{...}{additional arguments.}
-
-\item{plot}{\code{Logical scalar}. Whether to plot short term change.
-(Default: \code{FALSE})}
 }
 \value{
-\code{dataframe} or \code{plot} with \code{short term change}
+\code{dataframe} with \code{short term change}
 calculations.
 }
 \description{
@@ -70,9 +66,9 @@ tse <- minimalgut
 short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
 
 # Subset samples by Time_lable and StudyIdentifier
-tse <- tse[, !(colData(tse)$Time_label \%in\% short_time_labels)]
-tse <- tse[, (colData(tse)$StudyIdentifier == "Bioreactor A")]
+tse <- tse[, !(tse$Time_label \%in\% short_time_labels)]
+tse <- tse[, (tse$StudyIdentifier == "Bioreactor A")]
 
-# Plot short term change in abundance
-shortTermChange(tse, rarefy = TRUE, plot = TRUE) + ggtitle("Bioreactor A")
+# Get short term change
+getShortTermChange(tse, rarefy = TRUE, time.col = "Time.hr")
 }
diff --git a/man/miaTime-package.Rd b/man/miaTime-package.Rd
index 41e150e..a2676ce 100644
--- a/man/miaTime-package.Rd
+++ b/man/miaTime-package.Rd
@@ -12,18 +12,18 @@
 \link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment} class
 }
 \author{
-\strong{Maintainer}: Yagmur Simsek \email{yagmur.simsek@hsrw.org}
+\strong{Maintainer}: Tuomas Borman \email{tuomas.v.borman@utu.fi} (\href{https://orcid.org/0000-0002-8563-8884}{ORCID})
 
 Authors:
 \itemize{
   \item Leo Lahti \email{leo.lahti@iki.fi} (\href{https://orcid.org/0000-0001-5537-637X}{ORCID})
+  \item Yagmur Simsek \email{yagmur.simsek@hsrw.org}
 }
 
 Other contributors:
 \itemize{
   \item Sudarshan Shetty \email{sudarshanshetty9@gmail.com} [contributor]
   \item Chouaib Benchraka [contributor]
-  \item Tuomas Borman \email{tuomas.v.borman@utu.fi} (\href{https://orcid.org/0000-0002-8563-8884}{ORCID}) [contributor]
   \item Muluh Muluh [contributor]
 }
 
diff --git a/tests/testthat/test-getShortTermChange.R b/tests/testthat/test-getShortTermChange.R
new file mode 100644
index 0000000..bab1915
--- /dev/null
+++ b/tests/testthat/test-getShortTermChange.R
@@ -0,0 +1,64 @@
+test_that("getShortTermChange", {
+    library(SummarizedExperiment)
+    # Load dataset
+    data(minimalgut)
+    tse <- minimalgut
+    # Check if the function handles empty input
+    empty_se <- SummarizedExperiment()
+    expect_error(getShortTermChange(empty_se), 
+                 "No data available in `x`")
+    # Check if assay.type argument works
+    # tse_invalid <- tse
+    # expect_error(
+    #     getShortTermChange(tse_invalid, assay.type = "invalid_assay"),
+    #     "'assay.type' must be a valid name of assays(x)"
+    # )
+    # Check that rarefy and compositional cannot both be TRUE
+    expect_error(getShortTermChange(tse, rarefy = TRUE, compositional = TRUE, 
+                                    time.col = "Time.hr"), 
+                 "Both rarefy and compositional cannot be TRUE simultaneously")
+    # Check if the depth argument is greater than minimum counts
+    min_depth <- min(assay(tse, "counts"))
+    expect_error(getShortTermChange(tse, depth = min_depth + 1, 
+                                    time.col = "Time.hr"),
+                 "Depth cannot be greater than the minimum number of counts in your data")
+    # Check if rarefy = TRUE works
+    result <- getShortTermChange(tse, rarefy = TRUE, time.col = "Time.hr")
+    expect_true(is.data.frame(result))
+    # Check if compositional = TRUE works
+    result <- getShortTermChange(tse, compositional = TRUE, time.col = "Time.hr")
+    expect_true(is.data.frame(result))  
+    # Should still return a dataframe
+    result <- getShortTermChange(tse, rarefy = TRUE, 
+                              compositional = FALSE, 
+                              time.col = "Time.hr")
+    expect_true(is.data.frame(result))  
+    
+    short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
+    # Subset samples by Time_label and StudyIdentifier 
+    tse_filtered <- tse[, !(tse$Time_label %in% short_time_labels)]
+    tse_filtered <- tse_filtered[, (tse_filtered$StudyIdentifier == "Bioreactor A")]
+    
+    expect_true(all(!(tse_filtered$Time_label %in% short_time_labels)))
+    
+    result <- getShortTermChange(tse_filtered, 
+                              rarefy = TRUE, time.col = "Time.hr")
+    
+    result <- getShortTermChange(tse_filtered, 
+                              compositional = TRUE, time.col = "Time.hr")
+    # Expected output is a dataframe
+    expect_true(is.data.frame(result))  
+    expect_true("growth_diff" %in% colnames(result))
+    # Test some expected properties (e.g., that growth_diff isn't all NAs)
+    expect_false(all(is.na(result$growth_diff)))
+    
+    min_depth <- min(assay(tse_filtered, "counts"))
+    result <- getShortTermChange(tse_filtered, rarefy = TRUE, 
+                                 depth = min_depth, time.col = "Time.hr")
+    expect_true(is.data.frame(result)) 
+    expect_error(getShortTermChange(tse_filtered, 
+                                 rarefy = TRUE, 
+                                 depth = min_depth + 1, time.col = "Time.hr"),
+                 "Depth cannot be greater than the minimum number of counts in your data")
+    
+})
diff --git a/tests/testthat/test-shortTermChange.R b/tests/testthat/test-shortTermChange.R
deleted file mode 100644
index a76ee92..0000000
--- a/tests/testthat/test-shortTermChange.R
+++ /dev/null
@@ -1,80 +0,0 @@
-test_that("shortTermChange", {
-    library(SummarizedExperiment)
-    # Load dataset
-    data(minimalgut)
-    tse <- minimalgut
-    # Check if the function handles empty input
-    empty_se <- SummarizedExperiment()
-    expect_error(shortTermChange(empty_se), 
-                 "No data available in `x`")
-    # Check if assay.type argument works
-    # tse_invalid <- tse
-    # expect_error(
-    #     shortTermChange(tse_invalid, assay.type = "invalid_assay"),
-    #     "'assay.type' must be a valid name of assays(x)"
-    # )
-    # Check that rarefy and compositional cannot both be TRUE
-    expect_error(shortTermChange(tse, rarefy = TRUE, compositional = TRUE), 
-                 "Both rarefy and compositional cannot be TRUE simultaneously")
-    # Check if the depth argument is greater than minimum counts
-    min_depth <- min(assay(tse, "counts"))
-    expect_error(shortTermChange(tse, depth = min_depth + 1),
-                 "Depth cannot be greater than the minimum number of counts in your data")
-    # Check if rarefy = TRUE works
-    result <- shortTermChange(tse, rarefy = TRUE)
-    expect_true(is.data.frame(result))
-    # Check if compositional = TRUE works
-    result <- shortTermChange(tse, compositional = TRUE)
-    expect_true(is.data.frame(result))  
-    # Expected output should be a ggplot object
-    result <- shortTermChange(tse, plot = TRUE)
-    expect_s3_class(result, "gg")  
-    # Should still return a dataframe
-    result <- shortTermChange(tse, rarefy = TRUE, 
-                              compositional = FALSE, 
-                              additional_arg = "value")
-    expect_true(is.data.frame(result))  
-    
-    short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
-    # Subset samples by Time_label and StudyIdentifier 
-    tse_filtered <- tse[, !(colData(tse)$Time_label %in% short_time_labels)]
-    tse_filtered <- tse_filtered[, (colData(tse_filtered)$StudyIdentifier == "Bioreactor A")]
-    
-    expect_true(all(!(colData(tse_filtered)$Time_label %in% short_time_labels)))
-    # Expected output should be a ggplot object
-    result <- shortTermChange(tse_filtered, 
-                              rarefy = TRUE, plot = TRUE)
-    expect_s3_class(result, "gg")  
-    
-    result <- shortTermChange(tse_filtered, 
-                              rarefy = FALSE, plot = FALSE)
-    # Expected output is a dataframe
-    expect_true(is.data.frame(result))  
-    # Ensure dataframe has expected columns
-    expect_true("OTU" %in% colnames(result))  
-    expect_true("growth_diff" %in% colnames(result))
-    
-    result <- shortTermChange(tse_filtered, 
-                              rarefy = TRUE, plot = FALSE)
-    
-    result <- shortTermChange(tse_filtered, 
-                              compositional = TRUE, plot = FALSE)
-    # Expected output is a dataframe
-    expect_true(is.data.frame(result))  
-    
-    result <- shortTermChange(tse_filtered, rarefy = FALSE, 
-                              compositional = FALSE, 
-                              plot = FALSE)
-    expect_true("growth_diff" %in% colnames(result))
-    # Test some expected properties (e.g., that growth_diff isn't all NAs)
-    expect_false(all(is.na(result$growth_diff)))
-    
-    min_depth <- min(assay(tse_filtered, "counts"))
-    result <- shortTermChange(tse_filtered, rarefy = TRUE, depth = min_depth)
-    expect_true(is.data.frame(result)) 
-    expect_error(shortTermChange(tse_filtered, 
-                                 rarefy = TRUE, 
-                                 depth = min_depth + 1),
-                 "Depth cannot be greater than the minimum number of counts in your data")
-    
-})
\ No newline at end of file

From 99b98e32a837af1c77ccc10a8eacd4e9743becab Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Fri, 8 Nov 2024 13:54:20 +0200
Subject: [PATCH 09/16] update

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 NAMESPACE              | 1 +
 R/getShortTermChange.R | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/NAMESPACE b/NAMESPACE
index 0060ed1..e576a75 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,6 +13,7 @@ exportMethods(getStepwiseDivergence)
 import(TreeSummarizedExperiment)
 import(mia)
 importFrom(SummarizedExperiment,colData)
+importFrom(dplyr,"%>%")
 importFrom(dplyr,arrange)
 importFrom(dplyr,as_tibble)
 importFrom(dplyr,group_by)
diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
index 1aacccd..1c87174 100644
--- a/R/getShortTermChange.R
+++ b/R/getShortTermChange.R
@@ -54,7 +54,7 @@ setGeneric("getShortTermChange", signature = c("x"),
 
 #' @rdname getShortTermChange
 #' @export
-#' @importFrom dplyr arrange as_tibble summarize
+#' @importFrom dplyr arrange as_tibble summarize "%>%"
 #' @importFrom ggplot2 ggplot aes
 #' @importFrom ggrepel geom_text_repel
 #' @importFrom mia rarefyAssay transformAssay

From 9529271a1659f782d8dd71e78ab0234d15fb124e Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Wed, 13 Nov 2024 12:47:20 +0200
Subject: [PATCH 10/16] update

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 NAMESPACE                                     |   2 +
 R/AllGenerics.R                               |  12 ++
 R/getShortTermChange.R                        | 112 ++++++++----------
 R/utils.R                                     |   1 +
 ...ortTermChange.Rd => addShortTermChange.Rd} |  59 +++++----
 tests/testthat/test-getShortTermChange.R      |  43 +------
 6 files changed, 91 insertions(+), 138 deletions(-)
 rename man/{getShortTermChange.Rd => addShortTermChange.Rd} (50%)

diff --git a/NAMESPACE b/NAMESPACE
index e576a75..d750737 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,11 +1,13 @@
 # Generated by roxygen2: do not edit by hand
 
 export(addBaselineDivergence)
+export(addShortTermChange)
 export(addStepwiseDivergence)
 export(getBaselineDivergence)
 export(getShortTermChange)
 export(getStepwiseDivergence)
 exportMethods(addBaselineDivergence)
+exportMethods(addShortTermChange)
 exportMethods(addStepwiseDivergence)
 exportMethods(getBaselineDivergence)
 exportMethods(getShortTermChange)
diff --git a/R/AllGenerics.R b/R/AllGenerics.R
index 55b5d16..5b57da5 100644
--- a/R/AllGenerics.R
+++ b/R/AllGenerics.R
@@ -20,3 +20,15 @@ setGeneric("getStepwiseDivergence", signature = c("x"), function(x, ...)
 #' @export
 setGeneric("addStepwiseDivergence", signature = "x", function(x, ...)
     standardGeneric("addStepwiseDivergence"))
+
+#' @rdname addShortTermChange
+#' @export
+setGeneric("addShortTermChange", signature = "x", function(x, ...)
+    standardGeneric("addShortTermChange"))
+
+#' @rdname addShortTermChange
+#' @export
+setGeneric("getShortTermChange", signature = "x", function(x, ...)
+    standardGeneric("getShortTermChange"))
+
+
diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
index 1c87174..106896f 100644
--- a/R/getShortTermChange.R
+++ b/R/getShortTermChange.R
@@ -6,17 +6,16 @@
 #' @param x a \code{\link{SummarizedExperiment}} object.
 #' @param assay.type \code{Character scalar}. Specifies the name of assay 
 #'   used in calculation. (Default: \code{"counts"})
-#' @param rarefy \code{Logical scalar}. Whether to rarefy counts.
-#' (Default: \code{FALSE})
-#' @param compositional \code{Logical scalar}. Whether to transform counts.
-#' (Default: \code{FALSE})
-#' @param depth \code{Integer scalar}. Specifies the depth used in rarefying. 
-#' (Default: \code{min(assay(x, assay.type)}))
+#' @param name \code{Character scalar}. Specifies a name for storing
+#' short term results. (Default: \code{"short_term_change"})
 #' @param ... additional arguments.
 #' 
 #' 
-#' @return \code{dataframe} with \code{short term change} 
-#' calculations.
+#' @return code{getShortTermChange} returns \code{DataFrame} object containing 
+#' the short term change in abundance over time for a microbe. 
+#' \code{addShortTermChange}, on the other hand, returns a
+#' \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}
+#' object with these results in its \code{metadata}.
 #' 
 #' @details This approach is used by Wisnoski NI and colleagues
 #'          \url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
@@ -26,7 +25,7 @@
 #'          \url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
 #'          This approach is useful for identifying short term growth behaviors of taxa.
 #'          
-#' @name getShortTermChange
+#' @name addShortTermChange
 #' 
 #' 
 #' @examples
@@ -37,22 +36,21 @@
 #' 
 #' short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
 #' 
-#' # Subset samples by Time_lable and StudyIdentifier
+#' # Subset samples by Time_label and StudyIdentifier
 #' tse <- tse[, !(tse$Time_label %in% short_time_labels)]
 #' tse <- tse[, (tse$StudyIdentifier == "Bioreactor A")]
 #' 
 #' # Get short term change
-#' getShortTermChange(tse, rarefy = TRUE, time.col = "Time.hr")
+#' # Case of rarefying counts
+#' tse <- transformAssay(tse, method = "relabundance")
+#' getShortTermChange(tse, assay.type = relabundance, time.col = "Time.hr")
+#' 
+#' # Case of transforming counts
+#' tse <- rarefyAssay(tse, assay.type = "counts")
+#' getShortTermChange(tse, assay.type = subsampled, time.col = "Time.hr")
 NULL
 
-#' @rdname getShortTermChange
-#' @export
-setGeneric("getShortTermChange", signature = c("x"),
-    function( x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
-        depth = min(assay(x, assay.type)), ...)
-        standardGeneric("getShortTermChange"))
-
-#' @rdname getShortTermChange
+#' @rdname addShortTermChange
 #' @export
 #' @importFrom dplyr arrange as_tibble summarize "%>%"
 #' @importFrom ggplot2 ggplot aes
@@ -60,61 +58,46 @@ setGeneric("getShortTermChange", signature = c("x"),
 #' @importFrom mia rarefyAssay transformAssay
 #' @importFrom SummarizedExperiment colData
 setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
-    function(x, assay.type = "counts", rarefy = FALSE, compositional = FALSE, 
-        depth = min(assay(x, assay.type)), ...){
+    function(x, assay.type = "counts", ...){
         ############################## Input check #############################
         # Check validity of object
-        if(nrow(x) == 0L){
+        if (nrow(x) == 0L){
             stop("No data available in `x` ('x' has nrow(x) == 0L.)",
                  call. = FALSE)
         }
         # Check assay.type
         .check_assay_present(assay.type, x)
-        if(!.is_a_bool(rarefy)){
-            stop("'rarefy' must be TRUE or FALSE.", call. = FALSE)
-        }
-        if(!.is_a_bool(compositional)){
-            stop("'compositional' must be TRUE or FALSE.", call. = FALSE)
-        }
-        # Ensure that the provided depth is valid
-        if ( !is.null(depth) && depth > min(assay(x, assay.type), na.rm = TRUE) ) {
-            stop("Depth cannot be greater than the minimum number of counts in your data", call. = FALSE)
-        
-        }
         ########################### Growth Metrics ############################
-        grwt <- .calculate_growth_metrics(x, assay.type = assay.type, 
-                                          rarefy = rarefy, 
-                                          compositional = compositional, 
-                                          depth = depth, ...)
-        # Clean and format growth matrics
+        grwt <- .calculate_growth_metrics(x, assay.type = assay.type, ...)
+        # Clean and format growth metrics
         grs.all <- .clean_growth_metrics(grwt, ...)
         return(grs.all)
     }
 )
+
+#' @rdname addShortTermChange
+#' @export
+setMethod("addShortTermChange", signature = c(x = "SummarizedExperiment"),
+          function(x, assay.type = "counts", name = "short_term_change", ...){
+              # Calculate short term change
+              res <- getShortTermChange(x, ...)
+              # Add to metadata
+              x <- .add_values_to_metadata(x, name, res, ...)
+              return(x)
+          }
+)
+
+################################ HELP FUNCTIONS ################################
+
 # wrapper to calculate growth matrix
-.calculate_growth_metrics <- function(x, assay.type, time.col = NULL, 
-                                      rarefy, compositional, depth, ...) {
-    ############################ Data Preparation ##############################
-    # Initialize the filtered object based on rarefy and compositional arguments
-    if (rarefy == TRUE && compositional == FALSE) {
-        message("rarefy is set to TRUE, calculating short term change using counts")
-        x <- rarefyAssay(x, assay.type = assay.type, depth = depth)
-        assay.type <- "subsampled"
-    } else if (rarefy == FALSE && compositional == FALSE) {
-        message("rarefy is set to FALSE, compositional==FALSE, using raw counts")
-        x <- x
-    } else if (rarefy == FALSE && compositional == TRUE) {
-        message("rarefy is set to FALSE, compositional==TRUE, using relative abundances")
-        x <- transformAssay(x, method = "relabundance", assay.type = assay.type)
-        assay.type <- "relabundance"
-    } else if (rarefy == TRUE && compositional == TRUE) {
-        stop("Both rarefy and compositional cannot be TRUE simultaneously", call. = FALSE)
-    }
-    # Reshape data and calcualte grwoth metrics
-    assay_data <- meltSE(x, assay.type = assay.type, add.col = time.col)
+.calculate_growth_metrics <- function(x, assay.type, time.col = NULL, ...) {
+
+    # Reshape data and calculate growth metrics
+    assay_data <- meltSE(x, assay.type = assay.type, 
+                         add.col = time.col, row.name = "Feature_ID")
     assay_data <- assay_data %>%
         arrange( !!sym(time.col) ) %>%
-        group_by(FeatureID) %>%
+        group_by(Feature_ID) %>%
         mutate(
             time_lag = !!sym(time.col) - lag( !!sym(time.col) ), 
             growth_diff =!!sym(assay.type) - lag(!!sym(assay.type)),
@@ -127,19 +110,18 @@ setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
 .clean_growth_metrics <- function(grwt, time.col = NULL, ...) {
     # Calculate max growth
     maxgrs <- grwt %>%
-        summarize(max.growth = max(growth_diff, na.rm = TRUE))
+        summarize(max_growth = max(growth_diff, na.rm = TRUE))
     # Merge growth data with max growth
-    grs.all <- merge(grwt, maxgrs, by = "FeatureID")
+    grs.all <- merge(grwt, maxgrs, by = "Feature_ID")
     # Add 'ismax' column indicating if the growth is the maximum
     grs.all <- grs.all %>%
-        mutate(ismax = ifelse(growth_diff == max.growth, TRUE, FALSE))
+        mutate(ismax = ifelse(growth_diff == max_growth, TRUE, FALSE))
     # Clean and abbreviate FeatureID names
-    grs.all$FeatureID <- gsub("_", " ", grs.all$FeatureID)
-    grs.all$FeatureIDabb <- toupper(abbreviate(grs.all$FeatureID, 
+    grs.all$Feature_IDabb <- toupper(abbreviate(grs.all$Feature_ID, 
                                                minlength = 3, 
                                                method = "both.sides"))
     # Create 'Feature.time' column combining abbreviation and time information
-    grs.all$Feature.time <- paste0(grs.all$FeatureIDabb, " ", 
+    grs.all$Feature_time <- paste0(grs.all$Feature_IDabb, " ", 
                                    grs.all[[time.col]], "h")
     
     return(grs.all)
diff --git a/R/utils.R b/R/utils.R
index 6a4d25b..1cbd6bf 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -166,3 +166,4 @@
 .check_assay_present <- mia:::.check_assay_present
 .add_values_to_colData <- mia:::.add_values_to_colData
 .check_and_get_altExp <- mia:::.check_and_get_altExp
+.add_values_to_metadata <- mia:::.add_values_to_metadata
diff --git a/man/getShortTermChange.Rd b/man/addShortTermChange.Rd
similarity index 50%
rename from man/getShortTermChange.Rd
rename to man/addShortTermChange.Rd
index a7e1e23..9e08b6a 100644
--- a/man/getShortTermChange.Rd
+++ b/man/addShortTermChange.Rd
@@ -1,48 +1,37 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/getShortTermChange.R
-\name{getShortTermChange}
+% Please edit documentation in R/AllGenerics.R, R/getShortTermChange.R
+\name{addShortTermChange}
+\alias{addShortTermChange}
 \alias{getShortTermChange}
 \alias{getShortTermChange,SummarizedExperiment-method}
+\alias{addShortTermChange,SummarizedExperiment-method}
 \title{Short term Changes in Abundance}
 \usage{
-getShortTermChange(
-  x,
-  assay.type = "counts",
-  rarefy = FALSE,
-  compositional = FALSE,
-  depth = min(assay(x, assay.type)),
-  ...
-)
+addShortTermChange(x, ...)
 
-\S4method{getShortTermChange}{SummarizedExperiment}(
-  x,
-  assay.type = "counts",
-  rarefy = FALSE,
-  compositional = FALSE,
-  depth = min(assay(x, assay.type)),
-  ...
-)
+getShortTermChange(x, ...)
+
+\S4method{getShortTermChange}{SummarizedExperiment}(x, assay.type = "counts", ...)
+
+\S4method{addShortTermChange}{SummarizedExperiment}(x, assay.type = "counts", name = "short_term_change", ...)
 }
 \arguments{
 \item{x}{a \code{\link{SummarizedExperiment}} object.}
 
+\item{...}{additional arguments.}
+
 \item{assay.type}{\code{Character scalar}. Specifies the name of assay
 used in calculation. (Default: \code{"counts"})}
 
-\item{rarefy}{\code{Logical scalar}. Whether to rarefy counts.
-(Default: \code{FALSE})}
-
-\item{compositional}{\code{Logical scalar}. Whether to transform counts.
-(Default: \code{FALSE})}
-
-\item{depth}{\code{Integer scalar}. Specifies the depth used in rarefying.
-(Default: \code{min(assay(x, assay.type)}))}
-
-\item{...}{additional arguments.}
+\item{name}{\code{Character scalar}. Specifies a name for storing
+short term results. (Default: \code{"short_term_change"})}
 }
 \value{
-\code{dataframe} with \code{short term change}
-calculations.
+code{getShortTermChange} returns \code{DataFrame} object containing
+the short term change in abundance over time for a microbe.
+\code{addShortTermChange}, on the other hand, returns a
+\code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}
+object with these results in its \code{metadata}.
 }
 \description{
 Calculates short term changes in abundance of taxa
@@ -65,10 +54,16 @@ tse <- minimalgut
 
 short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
 
-# Subset samples by Time_lable and StudyIdentifier
+# Subset samples by Time_label and StudyIdentifier
 tse <- tse[, !(tse$Time_label \%in\% short_time_labels)]
 tse <- tse[, (tse$StudyIdentifier == "Bioreactor A")]
 
 # Get short term change
-getShortTermChange(tse, rarefy = TRUE, time.col = "Time.hr")
+# Case of rarefying counts
+tse <- transformAssay(tse, method = "relabundance")
+getShortTermChange(tse, assay.type = relabundance, time.col = "Time.hr")
+
+# Case of transforming counts
+tse <- rarefyAssay(tse, assay.type = "counts")
+getShortTermChange(tse, assay.type = subsampled, time.col = "Time.hr")
 }
diff --git a/tests/testthat/test-getShortTermChange.R b/tests/testthat/test-getShortTermChange.R
index bab1915..26ce8c8 100644
--- a/tests/testthat/test-getShortTermChange.R
+++ b/tests/testthat/test-getShortTermChange.R
@@ -7,33 +7,8 @@ test_that("getShortTermChange", {
     empty_se <- SummarizedExperiment()
     expect_error(getShortTermChange(empty_se), 
                  "No data available in `x`")
-    # Check if assay.type argument works
-    # tse_invalid <- tse
-    # expect_error(
-    #     getShortTermChange(tse_invalid, assay.type = "invalid_assay"),
-    #     "'assay.type' must be a valid name of assays(x)"
-    # )
-    # Check that rarefy and compositional cannot both be TRUE
-    expect_error(getShortTermChange(tse, rarefy = TRUE, compositional = TRUE, 
-                                    time.col = "Time.hr"), 
-                 "Both rarefy and compositional cannot be TRUE simultaneously")
-    # Check if the depth argument is greater than minimum counts
-    min_depth <- min(assay(tse, "counts"))
-    expect_error(getShortTermChange(tse, depth = min_depth + 1, 
-                                    time.col = "Time.hr"),
-                 "Depth cannot be greater than the minimum number of counts in your data")
-    # Check if rarefy = TRUE works
-    result <- getShortTermChange(tse, rarefy = TRUE, time.col = "Time.hr")
-    expect_true(is.data.frame(result))
-    # Check if compositional = TRUE works
-    result <- getShortTermChange(tse, compositional = TRUE, time.col = "Time.hr")
-    expect_true(is.data.frame(result))  
+
     # Should still return a dataframe
-    result <- getShortTermChange(tse, rarefy = TRUE, 
-                              compositional = FALSE, 
-                              time.col = "Time.hr")
-    expect_true(is.data.frame(result))  
-    
     short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
     # Subset samples by Time_label and StudyIdentifier 
     tse_filtered <- tse[, !(tse$Time_label %in% short_time_labels)]
@@ -41,24 +16,10 @@ test_that("getShortTermChange", {
     
     expect_true(all(!(tse_filtered$Time_label %in% short_time_labels)))
     
-    result <- getShortTermChange(tse_filtered, 
-                              rarefy = TRUE, time.col = "Time.hr")
-    
-    result <- getShortTermChange(tse_filtered, 
-                              compositional = TRUE, time.col = "Time.hr")
+    result <- getShortTermChange(tse_filtered, time.col = "Time.hr")
     # Expected output is a dataframe
     expect_true(is.data.frame(result))  
     expect_true("growth_diff" %in% colnames(result))
     # Test some expected properties (e.g., that growth_diff isn't all NAs)
     expect_false(all(is.na(result$growth_diff)))
-    
-    min_depth <- min(assay(tse_filtered, "counts"))
-    result <- getShortTermChange(tse_filtered, rarefy = TRUE, 
-                                 depth = min_depth, time.col = "Time.hr")
-    expect_true(is.data.frame(result)) 
-    expect_error(getShortTermChange(tse_filtered, 
-                                 rarefy = TRUE, 
-                                 depth = min_depth + 1, time.col = "Time.hr"),
-                 "Depth cannot be greater than the minimum number of counts in your data")
-    
 })

From b2afe4bfb9108e05ed7332ee5408cd34024d77b9 Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Wed, 13 Nov 2024 13:36:47 +0200
Subject: [PATCH 11/16] update

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 R/getShortTermChange.R    | 2 +-
 man/addShortTermChange.Rd | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
index 106896f..7cd4ff1 100644
--- a/R/getShortTermChange.R
+++ b/R/getShortTermChange.R
@@ -11,7 +11,7 @@
 #' @param ... additional arguments.
 #' 
 #' 
-#' @return code{getShortTermChange} returns \code{DataFrame} object containing 
+#' @return \code{getShortTermChange} returns \code{DataFrame} object containing 
 #' the short term change in abundance over time for a microbe. 
 #' \code{addShortTermChange}, on the other hand, returns a
 #' \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}
diff --git a/man/addShortTermChange.Rd b/man/addShortTermChange.Rd
index 9e08b6a..1362e20 100644
--- a/man/addShortTermChange.Rd
+++ b/man/addShortTermChange.Rd
@@ -27,7 +27,7 @@ used in calculation. (Default: \code{"counts"})}
 short term results. (Default: \code{"short_term_change"})}
 }
 \value{
-code{getShortTermChange} returns \code{DataFrame} object containing
+\code{getShortTermChange} returns \code{DataFrame} object containing
 the short term change in abundance over time for a microbe.
 \code{addShortTermChange}, on the other hand, returns a
 \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}

From c34b0d7cc8d914b6c7e6ccbcbfd0dc69c4448a44 Mon Sep 17 00:00:00 2001
From: Daena Rys <rysdaena8@gmail.com>
Date: Wed, 13 Nov 2024 14:03:25 +0200
Subject: [PATCH 12/16] update

Signed-off-by: Daena Rys <rysdaena8@gmail.com>
---
 R/getShortTermChange.R    | 4 ++--
 man/addShortTermChange.Rd | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
index 7cd4ff1..72aab05 100644
--- a/R/getShortTermChange.R
+++ b/R/getShortTermChange.R
@@ -43,11 +43,11 @@
 #' # Get short term change
 #' # Case of rarefying counts
 #' tse <- transformAssay(tse, method = "relabundance")
-#' getShortTermChange(tse, assay.type = relabundance, time.col = "Time.hr")
+#' getShortTermChange(tse, assay.type = "relabundance", time.col = "Time.hr")
 #' 
 #' # Case of transforming counts
 #' tse <- rarefyAssay(tse, assay.type = "counts")
-#' getShortTermChange(tse, assay.type = subsampled, time.col = "Time.hr")
+#' getShortTermChange(tse, assay.type = "subsampled", time.col = "Time.hr")
 NULL
 
 #' @rdname addShortTermChange
diff --git a/man/addShortTermChange.Rd b/man/addShortTermChange.Rd
index 1362e20..4653f23 100644
--- a/man/addShortTermChange.Rd
+++ b/man/addShortTermChange.Rd
@@ -61,9 +61,9 @@ tse <- tse[, (tse$StudyIdentifier == "Bioreactor A")]
 # Get short term change
 # Case of rarefying counts
 tse <- transformAssay(tse, method = "relabundance")
-getShortTermChange(tse, assay.type = relabundance, time.col = "Time.hr")
+getShortTermChange(tse, assay.type = "relabundance", time.col = "Time.hr")
 
 # Case of transforming counts
 tse <- rarefyAssay(tse, assay.type = "counts")
-getShortTermChange(tse, assay.type = subsampled, time.col = "Time.hr")
+getShortTermChange(tse, assay.type = "subsampled", time.col = "Time.hr")
 }

From f79bce66d2ddc9c9a7467462d7c574407900a676 Mon Sep 17 00:00:00 2001
From: TuomasBorman <tvborm@utu.fi>
Date: Fri, 24 Jan 2025 17:03:07 +0200
Subject: [PATCH 13/16] up

---
 R/getShortTermChange.R | 225 ++++++++++++++++++++++-------------------
 1 file changed, 122 insertions(+), 103 deletions(-)

diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
index 72aab05..9a331d6 100644
--- a/R/getShortTermChange.R
+++ b/R/getShortTermChange.R
@@ -1,128 +1,147 @@
-#' @title Short term Changes in Abundance
+#' @name
+#' getShortTermChange
+#'
+#' @export
+#'
+#' @title
+#' Short term changes in abundance
+#'
+#' @description
+#' Calculates short term changes in abundance of taxa using temporal
+#' abundance data.
+#'
+#' @details
+#' This approach is used by Wisnoski NI and colleagues
+#' \url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
+#' the following calculation log(present abundance/past abundance).
+#' Also a compositional version using relative abundance similar to
+#' Brian WJi, Sheth R et al
+#' \url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
+#' This approach is useful for identifying short term growth behaviors of taxa.
+#'
+#' @references
+#'
+#' Ji, B.W., et al. (2020) Macroecological dynamics of gut microbiota.
+#' Nat Microbiol 5, 768–775 . doi: https://doi.org/10.1038/s41564-020-0685-1
+#'
+#' @return
+#' \code{getShortTermChange} returns \code{DataFrame} object containing
+#' the short term change in abundance over time for a microbe.
+#' \code{addShortTermChange}, on the other hand, returns a
+#' \code{\link[SummarizedExperiment:SummarizedExperiment]{SummarizedExperiment}}
+#' object with these results in its \code{metadata}.
+#'
+#' @inheritParams addBaselineDivergence
 #'
-#' @description Calculates short term changes in abundance of taxa
-#'              using temporal Abundance data.
-#'  
-#' @param x a \code{\link{SummarizedExperiment}} object.
-#' @param assay.type \code{Character scalar}. Specifies the name of assay 
-#'   used in calculation. (Default: \code{"counts"})
 #' @param name \code{Character scalar}. Specifies a name for storing
 #' short term results. (Default: \code{"short_term_change"})
+#'
 #' @param ... additional arguments.
-#' 
-#' 
-#' @return \code{getShortTermChange} returns \code{DataFrame} object containing 
-#' the short term change in abundance over time for a microbe. 
-#' \code{addShortTermChange}, on the other hand, returns a
-#' \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}
-#' object with these results in its \code{metadata}.
-#' 
-#' @details This approach is used by Wisnoski NI and colleagues
-#'          \url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
-#'          the following calculation log(present abundance/past abundance).
-#'          Also a compositional version using relative abundance similar to
-#'          Brian WJi, Sheth R et al
-#'          \url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
-#'          This approach is useful for identifying short term growth behaviors of taxa.
-#'          
-#' @name addShortTermChange
-#' 
-#' 
+#
 #' @examples
-#' 
+#' library(miaTime)
+#'
 #' # Load time series data
 #' data(minimalgut)
 #' tse <- minimalgut
-#' 
-#' short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
-#' 
-#' # Subset samples by Time_label and StudyIdentifier
-#' tse <- tse[, !(tse$Time_label %in% short_time_labels)]
-#' tse <- tse[, (tse$StudyIdentifier == "Bioreactor A")]
-#' 
-#' # Get short term change
-#' # Case of rarefying counts
+#'
+#' # Get relative abundances
 #' tse <- transformAssay(tse, method = "relabundance")
-#' getShortTermChange(tse, assay.type = "relabundance", time.col = "Time.hr")
-#' 
-#' # Case of transforming counts
-#' tse <- rarefyAssay(tse, assay.type = "counts")
-#' getShortTermChange(tse, assay.type = "subsampled", time.col = "Time.hr")
+#' # Calculate short term changes
+#' df <- getShortTermChange(tse, assay.type = "relabundance", time.col = "Time.hr", group = "StudyIdentifier")
+#'
+#' # Select certain bacteria
+#' select <- df[["FeatureID"]] %in% c(
+#'     "Ruminococcus_bromii", "Coprococcus_catus", "Akkermansia_muciniphila")
+#' df <- df[ select, ]
+#'
+#' # Plot results
+#' library(ggplot2)
+#' p <- ggplot(df, aes(x = Time.hr, y = growth_rate, colour = FeatureID)) +
+#'     geom_smooth() +
+#'     facet_grid(. ~ StudyIdentifier, scales = "free") +
+#'     scale_y_log10()
+#' p
+#'
+#' @seealso
+#' \code{\link[mia:getStepwiseDivergence]{mia::getStepwiseDivergence()}}
 NULL
 
-#' @rdname addShortTermChange
+#' @rdname getShortTermChange
+#' @export
+setMethod("addShortTermChange", signature = c(x = "SummarizedExperiment"),
+    function(x, name = "short_term_change", ...){
+        x <- .check_and_get_altExp(x, ...)
+        args <- c(list(x = x), list(...)[!names(list(...)) %in% c("altexp")])
+        # Calculate short term change
+        res <- do.call(getShortTermChange, args)
+        # Add to metadata
+        x <- .add_values_to_metadata(x, name, res, ...)
+        return(x)
+    }
+)
+
+#' @rdname getShortTermChange
 #' @export
-#' @importFrom dplyr arrange as_tibble summarize "%>%"
-#' @importFrom ggplot2 ggplot aes
-#' @importFrom ggrepel geom_text_repel
-#' @importFrom mia rarefyAssay transformAssay
-#' @importFrom SummarizedExperiment colData
 setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
-    function(x, assay.type = "counts", ...){
+    function(
+        x, time.col, assay.type = "counts", time.interval = 1L, group = NULL,
+        ...){
         ############################## Input check #############################
-        # Check validity of object
-        if (nrow(x) == 0L){
-            stop("No data available in `x` ('x' has nrow(x) == 0L.)",
+        x <- .check_and_get_altExp(x, ...)
+        temp <- .check_input(
+            time.col, list("character scalar"), colnames(colData(x)))
+        if( !is.numeric(x[[time.col]]) ){
+            stop("'time.col' must specify numeric column from colData(x)",
                  call. = FALSE)
         }
-        # Check assay.type
+        #
         .check_assay_present(assay.type, x)
-        ########################### Growth Metrics ############################
-        grwt <- .calculate_growth_metrics(x, assay.type = assay.type, ...)
-        # Clean and format growth metrics
-        grs.all <- .clean_growth_metrics(grwt, ...)
-        return(grs.all)
+        #
+        temp <- .check_input(
+            group, list(NULL, "character scalar"), colnames(colData(x)))
+        ########################### Input check end ############################
+        res <- .calculate_growth_metrics(
+            x, assay.type, time.col, group, time.interval, ...)
+        return(res)
     }
 )
 
-#' @rdname addShortTermChange
-#' @export
-setMethod("addShortTermChange", signature = c(x = "SummarizedExperiment"),
-          function(x, assay.type = "counts", name = "short_term_change", ...){
-              # Calculate short term change
-              res <- getShortTermChange(x, ...)
-              # Add to metadata
-              x <- .add_values_to_metadata(x, name, res, ...)
-              return(x)
-          }
-)
-
 ################################ HELP FUNCTIONS ################################
 
 # wrapper to calculate growth matrix
-.calculate_growth_metrics <- function(x, assay.type, time.col = NULL, ...) {
-
-    # Reshape data and calculate growth metrics
-    assay_data <- meltSE(x, assay.type = assay.type, 
-                         add.col = time.col, row.name = "Feature_ID")
-    assay_data <- assay_data %>%
-        arrange( !!sym(time.col) ) %>%
-        group_by(Feature_ID) %>%
+#' @importFrom dplyr arrange group_by summarise mutate select
+.calculate_growth_metrics <- function(
+        x, assay.type, time.col, group, time.interval, ...) {
+    # Get data in long format
+    df <- meltSE(x, assay.type, add.col = c(time.col, group))
+    # If there are replicated samples, give warning that average is calculated
+    if( anyDuplicated(df[, c("FeatureID", group, time.col)]) ){
+        warning("The dataset contains replicated samples. The average is ",
+                "calculated for each time point.", call. = FALSE)
+    }
+    # Sort data based on time
+    df <- df |> arrange( !!sym(time.col) )
+    # Group based on features and time points. If group was specified, take that
+    # also into account
+    if( !is.null(group) ){
+        df <- df |> group_by(!!sym(group), FeatureID, !!sym(time.col))
+    } else{
+        df <- df |> group_by(FeatureID, !!sym(time.col))
+    }
+    df <- df |>
+        # Summarize duplicated samples by taking an average
+        summarise(value = mean(.data[[assay.type]], na.rm = TRUE)) |>
+        # For each feature in a sample group, calculate growth metrics
         mutate(
-            time_lag = !!sym(time.col) - lag( !!sym(time.col) ), 
-            growth_diff =!!sym(assay.type) - lag(!!sym(assay.type)),
-            growth_rate = (!!sym(assay.type) - lag(!!sym(assay.type))) / lag(!!sym(assay.type)),
-            var_abund = (!!sym(assay.type) - lag(!!sym(assay.type))) / time_lag
-        )
-    return(assay_data)
-}
-
-.clean_growth_metrics <- function(grwt, time.col = NULL, ...) {
-    # Calculate max growth
-    maxgrs <- grwt %>%
-        summarize(max_growth = max(growth_diff, na.rm = TRUE))
-    # Merge growth data with max growth
-    grs.all <- merge(grwt, maxgrs, by = "Feature_ID")
-    # Add 'ismax' column indicating if the growth is the maximum
-    grs.all <- grs.all %>%
-        mutate(ismax = ifelse(growth_diff == max_growth, TRUE, FALSE))
-    # Clean and abbreviate FeatureID names
-    grs.all$Feature_IDabb <- toupper(abbreviate(grs.all$Feature_ID, 
-                                               minlength = 3, 
-                                               method = "both.sides"))
-    # Create 'Feature.time' column combining abbreviation and time information
-    grs.all$Feature_time <- paste0(grs.all$Feature_IDabb, " ", 
-                                   grs.all[[time.col]], "h")
-    
-    return(grs.all)
+            time_lag = !!sym(time.col) -
+                lag(!!sym(time.col), n = time.interval),
+            growth_diff = value - lag(value, n = time.interval),
+            growth_rate = (value - lag(value, n = time.interval)) /
+                lag(value, n = time.interval),
+            var_abund = (value - lag(value, n = time.interval)) / time_lag
+        ) |>
+        # Remove value column that includes average abundances
+        select(-value)
+    return(df)
 }

From 5ed978442357ef198b8d990204c1671bff102909 Mon Sep 17 00:00:00 2001
From: TuomasBorman <tvborm@utu.fi>
Date: Sat, 25 Jan 2025 11:22:06 +0200
Subject: [PATCH 14/16] up

---
 NAMESPACE                                     |   7 +-
 R/AllGenerics.R                               |  12 +-
 R/getBaselineDivergence.R                     |   8 +-
 R/getShortTermChange.R                        |  61 ++++++----
 R/getStepwiseDivergence.R                     |   8 +-
 man/addShortTermChange.Rd                     |  69 -----------
 ...Divergence.Rd => getBaselineDivergence.Rd} |   0
 man/getShortTermChange.Rd                     | 108 ++++++++++++++++++
 ...Divergence.Rd => getStepwiseDivergence.Rd} |   0
 9 files changed, 168 insertions(+), 105 deletions(-)
 delete mode 100644 man/addShortTermChange.Rd
 rename man/{addBaselineDivergence.Rd => getBaselineDivergence.Rd} (100%)
 create mode 100644 man/getShortTermChange.Rd
 rename man/{addStepwiseDivergence.Rd => getStepwiseDivergence.Rd} (100%)

diff --git a/NAMESPACE b/NAMESPACE
index be4be9e..fb32bc4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,15 +2,19 @@
 
 export(addBaselineDivergence)
 export(addBimodality)
+export(addShortTermChange)
 export(addStepwiseDivergence)
 export(getBaselineDivergence)
 export(getBimodality)
+export(getShortTermChange)
 export(getStepwiseDivergence)
 exportMethods(addBaselineDivergence)
 exportMethods(addBimodality)
+exportMethods(addShortTermChange)
 exportMethods(addStepwiseDivergence)
 exportMethods(getBaselineDivergence)
 exportMethods(getBimodality)
+exportMethods(getShortTermChange)
 exportMethods(getStepwiseDivergence)
 import(S4Vectors)
 import(SingleCellExperiment)
@@ -21,10 +25,11 @@ import(mia)
 importFrom(dplyr,across)
 importFrom(dplyr,all_of)
 importFrom(dplyr,arrange)
-importFrom(dplyr,as_tibble)
 importFrom(dplyr,group_by)
 importFrom(dplyr,lag)
 importFrom(dplyr,mutate)
+importFrom(dplyr,select)
+importFrom(dplyr,summarise)
 importFrom(dplyr,summarize)
 importFrom(dplyr,ungroup)
 importFrom(tidyr,pivot_wider)
diff --git a/R/AllGenerics.R b/R/AllGenerics.R
index 7a35e9c..3e4f15e 100644
--- a/R/AllGenerics.R
+++ b/R/AllGenerics.R
@@ -1,22 +1,22 @@
 # All generic methods are listed here
 
-#' @rdname addBaselineDivergence
+#' @rdname getBaselineDivergence
 #' @export
 setGeneric("getBaselineDivergence", signature = "x", function(x, ...)
     standardGeneric("getBaselineDivergence"))
 
-#' @rdname addBaselineDivergence
+#' @rdname getBaselineDivergence
 #' @export
 setGeneric("addBaselineDivergence", signature = "x", function(x, ...)
     standardGeneric("addBaselineDivergence"))
 
-#' @rdname addStepwiseDivergence
+#' @rdname getStepwiseDivergence
 #' @export
 #'
 setGeneric("getStepwiseDivergence", signature = c("x"), function(x, ...)
     standardGeneric("getStepwiseDivergence"))
 
-#' @rdname addStepwiseDivergence
+#' @rdname getStepwiseDivergence
 #' @export
 setGeneric("addStepwiseDivergence", signature = "x", function(x, ...)
     standardGeneric("addStepwiseDivergence"))
@@ -31,12 +31,12 @@ setGeneric("getBimodality", signature = "x", function(x, ...)
 setGeneric("addBimodality", signature = "x", function(x, ...)
     standardGeneric("addBimodality"))
 
-#' @rdname addShortTermChange
+#' @rdname getShortTermChange
 #' @export
 setGeneric("addShortTermChange", signature = "x", function(x, ...)
     standardGeneric("addShortTermChange"))
 
-#' @rdname addShortTermChange
+#' @rdname getShortTermChange
 #' @export
 setGeneric("getShortTermChange", signature = "x", function(x, ...)
     standardGeneric("getShortTermChange"))
diff --git a/R/getBaselineDivergence.R b/R/getBaselineDivergence.R
index a3c3fa6..644832e 100644
--- a/R/getBaselineDivergence.R
+++ b/R/getBaselineDivergence.R
@@ -1,4 +1,6 @@
-#' @name addBaselineDivergence
+#' @name
+#' getBaselineDivergence
+#'
 #' @export
 #'
 #' @title
@@ -94,7 +96,7 @@
 #'
 NULL
 
-#' @rdname addBaselineDivergence
+#' @rdname getBaselineDivergence
 #' @export
 setMethod("getBaselineDivergence", signature = c(x = "SummarizedExperiment"),
     function(
@@ -155,7 +157,7 @@ setMethod("getBaselineDivergence", signature = c(x = "SummarizedExperiment"),
     }
 )
 
-#' @rdname addBaselineDivergence
+#' @rdname getBaselineDivergence
 #' @export
 setMethod("addBaselineDivergence", signature = c(x = "SummarizedExperiment"),
     function(
diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
index 9a331d6..1df5650 100644
--- a/R/getShortTermChange.R
+++ b/R/getShortTermChange.R
@@ -11,16 +11,19 @@
 #' abundance data.
 #'
 #' @details
-#' This approach is used by Wisnoski NI and colleagues
-#' \url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
-#' the following calculation log(present abundance/past abundance).
-#' Also a compositional version using relative abundance similar to
-#' Brian WJi, Sheth R et al
-#' \url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
-#' This approach is useful for identifying short term growth behaviors of taxa.
+#' These functions can be utilized to calculate growth metrics for short term
+#' change. In specific, the functions calculate the metrics with the following
+#' equations:
 #'
-#' @references
+#' \deqn{time\_diff = time_{t} - time_{t-1}}
+#'
+#' \deqn{abundance\_diff = abundance_{t} - abundance_{t-1}}
+#'
+#' \deqn{growth\_rate = abundance\_diff - abundance_{t-1}}
 #'
+#' \deqn{rate\_of\_change = abundance\_diff - time\_diff}
+#'
+#' @references
 #' Ji, B.W., et al. (2020) Macroecological dynamics of gut microbiota.
 #' Nat Microbiol 5, 768–775 . doi: https://doi.org/10.1038/s41564-020-0685-1
 #'
@@ -37,6 +40,12 @@
 #' short term results. (Default: \code{"short_term_change"})
 #'
 #' @param ... additional arguments.
+#' \itemize{
+#'   \item \code{time.interval}: \code{Integer scalar}. Indicates the increment
+#'   between time steps. By default, the function compares each sample to the
+#'   previous one. If you need to take every second, every third, or so, time
+#'   step, then increase this accordingly. (Default: \code{1L})
+#' }
 #
 #' @examples
 #' library(miaTime)
@@ -48,7 +57,15 @@
 #' # Get relative abundances
 #' tse <- transformAssay(tse, method = "relabundance")
 #' # Calculate short term changes
-#' df <- getShortTermChange(tse, assay.type = "relabundance", time.col = "Time.hr", group = "StudyIdentifier")
+#' df <- getShortTermChange(
+#'   tse, assay.type = "relabundance", time.col = "Time.hr",
+#'   group = "StudyIdentifier")
+#'
+#' # Calculate the logarithm of the ratio described in Ji, B.W., et al. (2020)
+#' tse <- transformAssay(
+#'     tse, assay.type = "relabundance", method = "log10", pseudocount = TRUE)
+#' df <- getShortTermChange(
+#'   tse, assay.type = "log10", time.col = "Time.hr", group = "StudyIdentifier")
 #'
 #' # Select certain bacteria
 #' select <- df[["FeatureID"]] %in% c(
@@ -84,9 +101,7 @@ setMethod("addShortTermChange", signature = c(x = "SummarizedExperiment"),
 #' @rdname getShortTermChange
 #' @export
 setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
-    function(
-        x, time.col, assay.type = "counts", time.interval = 1L, group = NULL,
-        ...){
+    function(x, time.col, assay.type = "counts", group = NULL, ...){
         ############################## Input check #############################
         x <- .check_and_get_altExp(x, ...)
         temp <- .check_input(
@@ -101,20 +116,21 @@ setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
         temp <- .check_input(
             group, list(NULL, "character scalar"), colnames(colData(x)))
         ########################### Input check end ############################
-        res <- .calculate_growth_metrics(
-            x, assay.type, time.col, group, time.interval, ...)
+        # Get data in long format
+        df <- meltSE(x, assay.type, add.col = c(time.col, group))
+        # Calculate metrics
+        res <- .calculate_growth_metrics(df, assay.type, time.col, group, ...)
         return(res)
     }
 )
 
 ################################ HELP FUNCTIONS ################################
 
-# wrapper to calculate growth matrix
+# This function calculates the growth metrics from the data.frame.
 #' @importFrom dplyr arrange group_by summarise mutate select
 .calculate_growth_metrics <- function(
-        x, assay.type, time.col, group, time.interval, ...) {
-    # Get data in long format
-    df <- meltSE(x, assay.type, add.col = c(time.col, group))
+        df, assay.type, time.col, group, time.interval = 1L, ...) {
+    temp <- .check_input(time.interval, list("numeric scalar"))
     # If there are replicated samples, give warning that average is calculated
     if( anyDuplicated(df[, c("FeatureID", group, time.col)]) ){
         warning("The dataset contains replicated samples. The average is ",
@@ -134,12 +150,11 @@ setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
         summarise(value = mean(.data[[assay.type]], na.rm = TRUE)) |>
         # For each feature in a sample group, calculate growth metrics
         mutate(
-            time_lag = !!sym(time.col) -
+            time_diff = !!sym(time.col) -
                 lag(!!sym(time.col), n = time.interval),
-            growth_diff = value - lag(value, n = time.interval),
-            growth_rate = (value - lag(value, n = time.interval)) /
-                lag(value, n = time.interval),
-            var_abund = (value - lag(value, n = time.interval)) / time_lag
+            abundance_diff = value - lag(value, n = time.interval),
+            growth_rate = abundance_diff / lag(value, n = time.interval),
+            rate_of_change = abundance_diff / time_diff
         ) |>
         # Remove value column that includes average abundances
         select(-value)
diff --git a/R/getStepwiseDivergence.R b/R/getStepwiseDivergence.R
index d2bef69..cd735d5 100644
--- a/R/getStepwiseDivergence.R
+++ b/R/getStepwiseDivergence.R
@@ -1,4 +1,6 @@
-#' @name addStepwiseDivergence
+#' @name
+#' getStepwiseDivergence
+#'
 #' @export
 #'
 #' @title
@@ -50,7 +52,7 @@
 #'
 NULL
 
-#' @rdname addStepwiseDivergence
+#' @rdname getStepwiseDivergence
 #' @export
 setMethod("getStepwiseDivergence", signature = c(x = "ANY"),
     function(
@@ -107,7 +109,7 @@ setMethod("getStepwiseDivergence", signature = c(x = "ANY"),
     }
 )
 
-#' @rdname addStepwiseDivergence
+#' @rdname getStepwiseDivergence
 #' @export
 setMethod("addStepwiseDivergence", signature = c(x = "SummarizedExperiment"),
     function(
diff --git a/man/addShortTermChange.Rd b/man/addShortTermChange.Rd
deleted file mode 100644
index 4653f23..0000000
--- a/man/addShortTermChange.Rd
+++ /dev/null
@@ -1,69 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/AllGenerics.R, R/getShortTermChange.R
-\name{addShortTermChange}
-\alias{addShortTermChange}
-\alias{getShortTermChange}
-\alias{getShortTermChange,SummarizedExperiment-method}
-\alias{addShortTermChange,SummarizedExperiment-method}
-\title{Short term Changes in Abundance}
-\usage{
-addShortTermChange(x, ...)
-
-getShortTermChange(x, ...)
-
-\S4method{getShortTermChange}{SummarizedExperiment}(x, assay.type = "counts", ...)
-
-\S4method{addShortTermChange}{SummarizedExperiment}(x, assay.type = "counts", name = "short_term_change", ...)
-}
-\arguments{
-\item{x}{a \code{\link{SummarizedExperiment}} object.}
-
-\item{...}{additional arguments.}
-
-\item{assay.type}{\code{Character scalar}. Specifies the name of assay
-used in calculation. (Default: \code{"counts"})}
-
-\item{name}{\code{Character scalar}. Specifies a name for storing
-short term results. (Default: \code{"short_term_change"})}
-}
-\value{
-\code{getShortTermChange} returns \code{DataFrame} object containing
-the short term change in abundance over time for a microbe.
-\code{addShortTermChange}, on the other hand, returns a
-\code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}
-object with these results in its \code{metadata}.
-}
-\description{
-Calculates short term changes in abundance of taxa
-using temporal Abundance data.
-}
-\details{
-This approach is used by Wisnoski NI and colleagues
-\url{https://github.com/nwisnoski/ul-seedbank}. Their approach is based on
-the following calculation log(present abundance/past abundance).
-Also a compositional version using relative abundance similar to
-Brian WJi, Sheth R et al
-\url{https://www.nature.com/articles/s41564-020-0685-1} can be used.
-This approach is useful for identifying short term growth behaviors of taxa.
-}
-\examples{
-
-# Load time series data
-data(minimalgut)
-tse <- minimalgut
-
-short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
-
-# Subset samples by Time_label and StudyIdentifier
-tse <- tse[, !(tse$Time_label \%in\% short_time_labels)]
-tse <- tse[, (tse$StudyIdentifier == "Bioreactor A")]
-
-# Get short term change
-# Case of rarefying counts
-tse <- transformAssay(tse, method = "relabundance")
-getShortTermChange(tse, assay.type = "relabundance", time.col = "Time.hr")
-
-# Case of transforming counts
-tse <- rarefyAssay(tse, assay.type = "counts")
-getShortTermChange(tse, assay.type = "subsampled", time.col = "Time.hr")
-}
diff --git a/man/addBaselineDivergence.Rd b/man/getBaselineDivergence.Rd
similarity index 100%
rename from man/addBaselineDivergence.Rd
rename to man/getBaselineDivergence.Rd
diff --git a/man/getShortTermChange.Rd b/man/getShortTermChange.Rd
new file mode 100644
index 0000000..53ba2f7
--- /dev/null
+++ b/man/getShortTermChange.Rd
@@ -0,0 +1,108 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/AllGenerics.R, R/getShortTermChange.R
+\name{addShortTermChange}
+\alias{addShortTermChange}
+\alias{getShortTermChange}
+\alias{addShortTermChange,SummarizedExperiment-method}
+\alias{getShortTermChange,SummarizedExperiment-method}
+\title{Short term changes in abundance}
+\usage{
+addShortTermChange(x, ...)
+
+getShortTermChange(x, ...)
+
+\S4method{addShortTermChange}{SummarizedExperiment}(x, name = "short_term_change", ...)
+
+\S4method{getShortTermChange}{SummarizedExperiment}(x, time.col, assay.type = "counts", group = NULL, ...)
+}
+\arguments{
+\item{x}{A
+\code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}
+object.}
+
+\item{...}{additional arguments.
+\itemize{
+\item \code{time.interval}: \code{Integer scalar}. Indicates the increment
+between time steps. By default, the function compares each sample to the
+previous one. If you need to take every second, every third, or so, time
+step, then increase this accordingly. (Default: \code{1L})
+}}
+
+\item{name}{\code{Character scalar}. Specifies a name for storing
+short term results. (Default: \code{"short_term_change"})}
+
+\item{time.col}{\code{Character scalar}. Specifies a name of the column from
+\code{colData} that identifies the sampling time points for the samples.}
+
+\item{assay.type}{\code{Character scalar}. Specifies which assay values are
+used in the dissimilarity estimation. (Default: \code{"counts"})}
+
+\item{group}{\code{Character scalar}. Specifies a name of the column from
+\code{colData} that identifies the grouping of the samples.
+(Default: \code{NULL})}
+}
+\value{
+\code{getShortTermChange} returns \code{DataFrame} object containing
+the short term change in abundance over time for a microbe.
+\code{addShortTermChange}, on the other hand, returns a
+\code{\link[SummarizedExperiment:SummarizedExperiment]{SummarizedExperiment}}
+object with these results in its \code{metadata}.
+}
+\description{
+Calculates short term changes in abundance of taxa using temporal
+abundance data.
+}
+\details{
+These functions can be utilized to calculate growth metrics for short term
+change. In specific, the functions calculate the metrics with the following
+equations:
+
+\deqn{time\_diff = time_{t} - time_{t-1}}
+
+\deqn{abundance\_diff = abundance_{t} - abundance_{t-1}}
+
+\deqn{growth\_rate = abundance\_diff - abundance_{t-1}}
+
+\deqn{rate\_of\_change = abundance\_diff - time\_diff}
+}
+\examples{
+library(miaTime)
+
+# Load time series data
+data(minimalgut)
+tse <- minimalgut
+
+# Get relative abundances
+tse <- transformAssay(tse, method = "relabundance")
+# Calculate short term changes
+df <- getShortTermChange(
+  tse, assay.type = "relabundance", time.col = "Time.hr",
+  group = "StudyIdentifier")
+
+# Calculate the logarithm of the ratio described in Ji, B.W., et al. (2020)
+tse <- transformAssay(
+    tse, assay.type = "relabundance", method = "log10", pseudocount = TRUE)
+df <- getShortTermChange(
+  tse, assay.type = "log10", time.col = "Time.hr", group = "StudyIdentifier")
+
+# Select certain bacteria
+select <- df[["FeatureID"]] \%in\% c(
+    "Ruminococcus_bromii", "Coprococcus_catus", "Akkermansia_muciniphila")
+df <- df[ select, ]
+
+# Plot results
+library(ggplot2)
+p <- ggplot(df, aes(x = Time.hr, y = growth_rate, colour = FeatureID)) +
+    geom_smooth() +
+    facet_grid(. ~ StudyIdentifier, scales = "free") +
+    scale_y_log10()
+p
+
+}
+\references{
+Ji, B.W., et al. (2020) Macroecological dynamics of gut microbiota.
+Nat Microbiol 5, 768–775 . doi: https://doi.org/10.1038/s41564-020-0685-1
+}
+\seealso{
+\code{\link[mia:getStepwiseDivergence]{mia::getStepwiseDivergence()}}
+}
diff --git a/man/addStepwiseDivergence.Rd b/man/getStepwiseDivergence.Rd
similarity index 100%
rename from man/addStepwiseDivergence.Rd
rename to man/getStepwiseDivergence.Rd

From c2e0ac09156367775d1d50af2d1be5968d669ddf Mon Sep 17 00:00:00 2001
From: TuomasBorman <tvborm@utu.fi>
Date: Sat, 25 Jan 2025 12:54:11 +0200
Subject: [PATCH 15/16] up

---
 DESCRIPTION                              |   3 +-
 NEWS                                     |   4 +-
 R/getBaselineDivergence.R                |   1 -
 R/getBimodality.R                        |   1 -
 R/getShortTermChange.R                   |  19 +-
 man/getBaselineDivergence.Rd             |   1 -
 man/getBimodality.Rd                     |   1 -
 man/getShortTermChange.Rd                |  10 +-
 tests/testthat/test-getShortTermChange.R | 236 ++++++++++++++++++++---
 9 files changed, 233 insertions(+), 43 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 3554575..9a4bcf0 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: miaTime
 Type: Package
 Title: Microbiome Time Series Analysis
-Version: 0.99.3
+Version: 0.99.4
 Authors@R:
      c(person(given = "Leo", family = "Lahti", role = c("aut"),
               email = "leo.lahti@iki.fi",
@@ -37,6 +37,7 @@ Imports:
 Suggests:
     BiocStyle,
     devtools,
+    ggplot2,
     knitr,
     lubridate,
     miaViz,
diff --git a/NEWS b/NEWS
index 4ba2d9f..238506d 100755
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,6 @@
+Changes in version 0.99.4
++ Added functions for calculating short term change metrics
+
 Changes in version 0.99.3
 + Handle repeated samples in divergence functions by calculating average
 
@@ -25,4 +28,3 @@ Changes in version 0.1.5
 + Changed the function name getTimeDivergence (deprecated) into getStepwiseDivergence
 + Added internal utility functions
 + Added SilvermanAGutData
-
diff --git a/R/getBaselineDivergence.R b/R/getBaselineDivergence.R
index 644832e..d542d4d 100644
--- a/R/getBaselineDivergence.R
+++ b/R/getBaselineDivergence.R
@@ -65,7 +65,6 @@
 #'
 #' @examples
 #' library(miaTime)
-#' library(mia)
 #'
 #' data(hitchip1006)
 #' tse <- transformAssay(hitchip1006, method = "relabundance")
diff --git a/R/getBimodality.R b/R/getBimodality.R
index 39873aa..17b4c8d 100644
--- a/R/getBimodality.R
+++ b/R/getBimodality.R
@@ -58,7 +58,6 @@
 #'
 #' @examples
 #' library(miaTime)
-#' library(mia)
 #'
 #' data(SilvermanAGutData)
 #' tse <- SilvermanAGutData
diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
index 1df5650..034876e 100644
--- a/R/getShortTermChange.R
+++ b/R/getShortTermChange.R
@@ -67,15 +67,15 @@
 #' df <- getShortTermChange(
 #'   tse, assay.type = "log10", time.col = "Time.hr", group = "StudyIdentifier")
 #'
-#' # Select certain bacteria
-#' select <- df[["FeatureID"]] %in% c(
-#'     "Ruminococcus_bromii", "Coprococcus_catus", "Akkermansia_muciniphila")
-#' df <- df[ select, ]
+#' # Select certain bacteria that have highest growth rate
+#' select <- df[["growth_rate"]] |> abs() |> order(decreasing = FALSE)
+#' select <- df[select, "FeatureID"] |> unique() |> head(3)
+#' df <- df[ which(df[["FeatureID"]] %in% select), ]
 #'
 #' # Plot results
 #' library(ggplot2)
 #' p <- ggplot(df, aes(x = Time.hr, y = growth_rate, colour = FeatureID)) +
-#'     geom_smooth() +
+#'     geom_smooth(level = 0.5) +
 #'     facet_grid(. ~ StudyIdentifier, scales = "free") +
 #'     scale_y_log10()
 #' p
@@ -88,6 +88,7 @@ NULL
 #' @export
 setMethod("addShortTermChange", signature = c(x = "SummarizedExperiment"),
     function(x, name = "short_term_change", ...){
+        temp <- .check_input(name, list("character scalar"))
         x <- .check_and_get_altExp(x, ...)
         args <- c(list(x = x), list(...)[!names(list(...)) %in% c("altexp")])
         # Calculate short term change
@@ -133,8 +134,8 @@ setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
     temp <- .check_input(time.interval, list("numeric scalar"))
     # If there are replicated samples, give warning that average is calculated
     if( anyDuplicated(df[, c("FeatureID", group, time.col)]) ){
-        warning("The dataset contains replicated samples. The average is ",
-                "calculated for each time point.", call. = FALSE)
+        message("The dataset contains replicated samples. The average ",
+                "abundance is calculated for each time point.")
     }
     # Sort data based on time
     df <- df |> arrange( !!sym(time.col) )
@@ -157,6 +158,8 @@ setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
             rate_of_change = abundance_diff / time_diff
         ) |>
         # Remove value column that includes average abundances
-        select(-value)
+        select(-value) |>
+        # Convert to DataFrame
+        DataFrame()
     return(df)
 }
diff --git a/man/getBaselineDivergence.Rd b/man/getBaselineDivergence.Rd
index d8f487f..83835f4 100644
--- a/man/getBaselineDivergence.Rd
+++ b/man/getBaselineDivergence.Rd
@@ -87,7 +87,6 @@ group (named list per group).
 }
 \examples{
 library(miaTime)
-library(mia)
 
 data(hitchip1006)
 tse <- transformAssay(hitchip1006, method = "relabundance")
diff --git a/man/getBimodality.Rd b/man/getBimodality.Rd
index 5bb53f3..d2730a2 100644
--- a/man/getBimodality.Rd
+++ b/man/getBimodality.Rd
@@ -69,7 +69,6 @@ bimodality and abundance to determine conditionally rare taxa (CRT).
 }
 \examples{
 library(miaTime)
-library(mia)
 
 data(SilvermanAGutData)
 tse <- SilvermanAGutData
diff --git a/man/getShortTermChange.Rd b/man/getShortTermChange.Rd
index 53ba2f7..dc06dc8 100644
--- a/man/getShortTermChange.Rd
+++ b/man/getShortTermChange.Rd
@@ -85,15 +85,15 @@ tse <- transformAssay(
 df <- getShortTermChange(
   tse, assay.type = "log10", time.col = "Time.hr", group = "StudyIdentifier")
 
-# Select certain bacteria
-select <- df[["FeatureID"]] \%in\% c(
-    "Ruminococcus_bromii", "Coprococcus_catus", "Akkermansia_muciniphila")
-df <- df[ select, ]
+# Select certain bacteria that have highest growth rate
+select <- df[["growth_rate"]] |> abs() |> order(decreasing = FALSE)
+select <- df[select, "FeatureID"] |> unique() |> head(3)
+df <- df[ which(df[["FeatureID"]] \%in\% select), ]
 
 # Plot results
 library(ggplot2)
 p <- ggplot(df, aes(x = Time.hr, y = growth_rate, colour = FeatureID)) +
-    geom_smooth() +
+    geom_smooth(level = 0.5) +
     facet_grid(. ~ StudyIdentifier, scales = "free") +
     scale_y_log10()
 p
diff --git a/tests/testthat/test-getShortTermChange.R b/tests/testthat/test-getShortTermChange.R
index 26ce8c8..ce9672c 100644
--- a/tests/testthat/test-getShortTermChange.R
+++ b/tests/testthat/test-getShortTermChange.R
@@ -1,25 +1,213 @@
-test_that("getShortTermChange", {
-    library(SummarizedExperiment)
-    # Load dataset
-    data(minimalgut)
-    tse <- minimalgut
-    # Check if the function handles empty input
-    empty_se <- SummarizedExperiment()
-    expect_error(getShortTermChange(empty_se), 
-                 "No data available in `x`")
-
-    # Should still return a dataframe
-    short_time_labels <- c("74.5h", "173h", "438h", "434h", "390h")
-    # Subset samples by Time_label and StudyIdentifier 
-    tse_filtered <- tse[, !(tse$Time_label %in% short_time_labels)]
-    tse_filtered <- tse_filtered[, (tse_filtered$StudyIdentifier == "Bioreactor A")]
-    
-    expect_true(all(!(tse_filtered$Time_label %in% short_time_labels)))
-    
-    result <- getShortTermChange(tse_filtered, time.col = "Time.hr")
-    # Expected output is a dataframe
-    expect_true(is.data.frame(result))  
-    expect_true("growth_diff" %in% colnames(result))
-    # Test some expected properties (e.g., that growth_diff isn't all NAs)
-    expect_false(all(is.na(result$growth_diff)))
+
+test_that("getShortTermChange errors", {
+    tse <- makeTSE(nrow = 100, ncol = 20)
+    assayNames(tse) <- "counts"
+    tse[["Time"]] <- sample(seq(1, 5), 20, replace = TRUE)
+
+    df <- getShortTermChange(tse, time.col = "Test", assay.type = "counts") |>
+        expect_error()
+    df <- getShortTermChange(tse, time.col = "group", assay.type = "counts") |>
+        expect_error()
+    df <- getShortTermChange(tse, time.col = "group", assay.type = "test") |>
+        expect_error()
+    df <- getShortTermChange(tse, time.col = "Time", group = "counts") |>
+        expect_error()
+    df <- getShortTermChange(tse, time.col = "Time", time.interval = TRUE) |>
+        expect_error()
+    tse <- addShortTermChange(tse, time.col = "Time", name = TRUE) |>
+        expect_error()
+})
+
+# Test with time interval 1 and no grouping
+test_that("getShortTermChange calculates correctly without group", {
+    tse <- makeTSE(nrow = 100, ncol = 20)
+    assayNames(tse) <- "counts"
+    tse[["Time"]] <- sample(seq(1, 5), 20, replace = TRUE)
+
+    df <- getShortTermChange(tse, time.col = "Time", assay.type = "counts") |>
+        expect_message()
+
+    # Check the structure of the returned object
+    expect_s4_class(df, "DataFrame")
+    expect_true("time_diff" %in% colnames(df))
+    expect_true("abundance_diff" %in% colnames(df))
+    expect_true("growth_rate" %in% colnames(df))
+    expect_true("rate_of_change" %in% colnames(df))
+    expect_true("Time" %in% colnames(df))
+    expect_true("FeatureID" %in% colnames(df))
+
+    # Test 10 times with random combinations
+    for( i in seq_len(10) ){
+        # Get random feature and timepoint to test
+        feat <- sample(rownames(tse), 1)
+        time_point <- sample(tse[["Time"]], 1)
+        # Get index of previous sample
+        temp_df <- colData(tse)[ order(tse[["Time"]]), ]
+        time_points <- unique(temp_df[["Time"]])
+        prev_time <- which(time_points == time_point)-1
+        ref <- rep(NA_real_, 4)
+        # Calculate reference values manually if previous sample exists
+        if( prev_time > 0 ){
+            prev_time <- time_points[ prev_time ]
+            prev_sam <- rownames(temp_df[temp_df[["Time"]] == prev_time, ])
+            sam <- rownames(temp_df[temp_df[["Time"]] == time_point, ])
+            # Take certain samples
+            tse_curr <- tse[, sam]
+            tse_prev <- tse[, prev_sam]
+            # Get values
+            time_diff <- time_point - prev_time
+            val1 <- mean(assay(tse_curr, "counts")[feat, ], na.rm = TRUE)
+            val2 <- mean(assay(tse_prev, "counts")[feat, ], na.rm = TRUE)
+            ref <- c(
+                time_diff,
+                val1 - val2,
+                (val1 - val2) / val2,
+                (val1 - val2) / time_diff
+            )
+        }
+        names(ref) <- c(
+            "time_diff", "abundance_diff", "growth_rate", "rate_of_change")
+        res <- df[df[["FeatureID"]] == feat & df[["Time"]] == time_point, ]
+        res <- unlist(res[, names(ref)])
+        #
+        expect_equal(res, ref)
+    }
+})
+
+# Test with time interval 2 and no grouping
+test_that("getShortTermChange calculates correctly with time interval 2", {
+    tse <- makeTSE(nrow = 1000, ncol = 200)
+    assayNames(tse) <- "counts"
+    tse[["Time"]] <- sample(seq(1, 5), 200, replace = TRUE)
+    tse <- transformAssay(tse, method = "relabundance")
+    df <- getShortTermChange(
+        tse, time.col = "Time", assay.type = "relabundance",
+        time.interval = 2) |>
+        expect_message()
+
+    # Check the structure of the returned object
+    expect_s4_class(df, "DataFrame")
+    expect_true("time_diff" %in% colnames(df))
+    expect_true("abundance_diff" %in% colnames(df))
+    expect_true("growth_rate" %in% colnames(df))
+    expect_true("rate_of_change" %in% colnames(df))
+    expect_true("Time" %in% colnames(df))
+    expect_true("FeatureID" %in% colnames(df))
+
+    # Test 10 times with random combinations
+    for( i in seq_len(10) ){
+        # Get random feature and timepoint to test
+        feat <- sample(rownames(tse), 1)
+        time_point <- sample(tse[["Time"]], 1)
+        # Get index of previous sample
+        temp_df <- colData(tse)[ order(tse[["Time"]]), ]
+        time_points <- unique(temp_df[["Time"]])
+        prev_time <- which(time_points == time_point)-2
+        ref <- rep(NA_real_, 4)
+        # Calculate reference values manually if previous sample exists
+        if( prev_time > 0 ){
+            prev_time <- time_points[ prev_time ]
+            prev_sam <- rownames(temp_df[temp_df[["Time"]] == prev_time, ])
+            sam <- rownames(temp_df[temp_df[["Time"]] == time_point, ])
+            # Take certain samples
+            tse_curr <- tse[, sam]
+            tse_prev <- tse[, prev_sam]
+            # Get values
+            time_diff <- time_point - prev_time
+            val1 <- mean(assay(tse_curr, "relabundance")[feat, ], na.rm = TRUE)
+            val2 <- mean(assay(tse_prev, "relabundance")[feat, ], na.rm = TRUE)
+            ref <- c(
+                time_diff,
+                val1 - val2,
+                (val1 - val2) / val2,
+                (val1 - val2) / time_diff
+            )
+        }
+        names(ref) <- c(
+            "time_diff", "abundance_diff", "growth_rate", "rate_of_change")
+        res <- df[df[["FeatureID"]] == feat & df[["Time"]] == time_point, ]
+        res <- unlist(res[, names(ref)])
+        #
+        expect_equal(res, ref)
+    }
+})
+
+# Test with grouping
+test_that("getShortTermChange calculates correctly with group", {
+    tse <- makeTSE(nrow = 1000, ncol = 200)
+    assayNames(tse) <- "counts"
+    tse[["Time"]] <- sample(seq(1, 5), 200, replace = TRUE)
+    # Remove duplicated samples
+    tse <- tse[ , !duplicated(colData(tse)[, c("group", "Time")])]
+
+    df <- getShortTermChange(
+        tse, time.col = "Time", assay.type = "counts", group = "group") |>
+        expect_no_message()
+
+    # Check the structure of the returned object
+    expect_s4_class(df, "DataFrame")
+    expect_true("time_diff" %in% colnames(df))
+    expect_true("abundance_diff" %in% colnames(df))
+    expect_true("growth_rate" %in% colnames(df))
+    expect_true("rate_of_change" %in% colnames(df))
+    expect_true("Time" %in% colnames(df))
+    expect_true("FeatureID" %in% colnames(df))
+
+    # Test 10 times with random combinations
+    for( i in seq_len(10) ){
+        # Get random feature and timepoint to test
+        feat <- sample(rownames(tse), 1)
+        time_point <- sample(tse[["Time"]], 1)
+        group <- sample(tse[["group"]], 1)
+        # Get only specific group
+        temp_df <- colData(tse)[ tse[["group"]] == group, ]
+        # Get index of previous sample
+        temp_df <- temp_df[ order(temp_df[["Time"]]), ]
+        time_points <- unique(temp_df[["Time"]])
+        prev_time <- which(time_points == time_point)-1
+        ref <- rep(NA_real_, 4)
+        # Calculate reference values manually if previous sample exists
+        if( prev_time > 0 ){
+            prev_time <- time_points[ prev_time ]
+            prev_sam <- rownames(temp_df[temp_df[["Time"]] == prev_time, ])
+            sam <- rownames(temp_df[temp_df[["Time"]] == time_point, ])
+            # Take certain samples
+            tse_curr <- tse[, sam]
+            tse_prev <- tse[, prev_sam]
+            # Get values
+            time_diff <- time_point - prev_time
+            val1 <- mean(assay(tse_curr, "counts")[feat, ], na.rm = TRUE)
+            val2 <- mean(assay(tse_prev, "counts")[feat, ], na.rm = TRUE)
+            ref <- c(
+                time_diff,
+                val1 - val2,
+                (val1 - val2) / val2,
+                (val1 - val2) / time_diff
+            )
+        }
+        names(ref) <- c(
+            "time_diff", "abundance_diff", "growth_rate", "rate_of_change")
+        select <- df[["FeatureID"]] == feat & df[["Time"]] == time_point &
+            df[["group"]] == group
+        res <- df[select, ]
+        res <- unlist(res[, names(ref)])
+        #
+        expect_equal(res, ref)
+    }
+})
+
+# Test that getShortTimeChange and addShortTimeChange are equal
+test_that("get* and add* are equal", {
+    tse <- makeTSE(nrow = 1000, ncol = 200)
+    assayNames(tse) <- "counts"
+    tse[["Time"]] <- sample(seq(1, 5), 200, replace = TRUE)
+    # Remove duplicated samples
+    tse <- tse[ , !duplicated(colData(tse)[, c("Time")])]
+
+    df <- getShortTermChange(tse, time.col = "Time") |>
+        expect_no_message()
+    tse <- addShortTermChange(tse, time.col = "Time", name = "test") |>
+        expect_no_message()
+
+    expect_equal(metadata(tse)[["test"]], df)
 })

From ed5c0987852b07b7e12b4008bfcbe09862d4020d Mon Sep 17 00:00:00 2001
From: TuomasBorman <tvborm@utu.fi>
Date: Sat, 25 Jan 2025 16:00:11 +0200
Subject: [PATCH 16/16] Fix check errors

---
 DESCRIPTION               |  4 +---
 NAMESPACE                 |  1 +
 R/getBaselineDivergence.R |  4 ++++
 R/getShortTermChange.R    | 12 +++++++-----
 R/miaTime.R               |  3 ++-
 man/getShortTermChange.Rd |  2 +-
 6 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 9a4bcf0..e3ea3a3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -23,9 +23,7 @@ biocViews: Microbiome, Software, Sequencing
 License: Artistic-2.0 | file LICENSE
 Depends:
     R (>= 4.4.0),
-    ggplot2,
-    mia,
-    SummarizedExperiment
+    mia
 Imports:
     dplyr,
     methods,
diff --git a/NAMESPACE b/NAMESPACE
index fb32bc4..b3a95c3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -31,6 +31,7 @@ importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(dplyr,summarise)
 importFrom(dplyr,summarize)
+importFrom(dplyr,sym)
 importFrom(dplyr,ungroup)
 importFrom(tidyr,pivot_wider)
 importFrom(tidyr,unnest)
diff --git a/R/getBaselineDivergence.R b/R/getBaselineDivergence.R
index d542d4d..be648b5 100644
--- a/R/getBaselineDivergence.R
+++ b/R/getBaselineDivergence.R
@@ -357,6 +357,10 @@ setMethod("addBaselineDivergence", signature = c(x = "SummarizedExperiment"),
 # This function get time difference between a sample and its reference sample
 #' @importFrom dplyr group_by mutate
 .get_time_difference <- function(x, time.col, reference, ...){
+    # This following line is to suppress "no visible binding for" messages
+    # in cmdcheck
+    temp_sample <- ":=" <- time <- .data <- NULL
+
     # Get timepoints
     time_point <- x[[time.col]]
     # Get reference time points
diff --git a/R/getShortTermChange.R b/R/getShortTermChange.R
index 034876e..3feb2b1 100644
--- a/R/getShortTermChange.R
+++ b/R/getShortTermChange.R
@@ -81,7 +81,7 @@
 #' p
 #'
 #' @seealso
-#' \code{\link[mia:getStepwiseDivergence]{mia::getStepwiseDivergence()}}
+#' \code{\link[miaTime:getStepwiseDivergence]{getStepwiseDivergence()}}
 NULL
 
 #' @rdname getShortTermChange
@@ -109,11 +109,9 @@ setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
             time.col, list("character scalar"), colnames(colData(x)))
         if( !is.numeric(x[[time.col]]) ){
             stop("'time.col' must specify numeric column from colData(x)",
-                 call. = FALSE)
+                call. = FALSE)
         }
-        #
         .check_assay_present(assay.type, x)
-        #
         temp <- .check_input(
             group, list(NULL, "character scalar"), colnames(colData(x)))
         ########################### Input check end ############################
@@ -128,9 +126,13 @@ setMethod("getShortTermChange", signature = c(x = "SummarizedExperiment"),
 ################################ HELP FUNCTIONS ################################
 
 # This function calculates the growth metrics from the data.frame.
-#' @importFrom dplyr arrange group_by summarise mutate select
+#' @importFrom dplyr arrange group_by sym summarise mutate select
 .calculate_growth_metrics <- function(
         df, assay.type, time.col, group, time.interval = 1L, ...) {
+    # This following line is to suppress "no visible binding for" messages
+    # in cmdcheck
+    FeatureID <- .data <- value <- abundance_diff <- time_diff <- NULL
+    #
     temp <- .check_input(time.interval, list("numeric scalar"))
     # If there are replicated samples, give warning that average is calculated
     if( anyDuplicated(df[, c("FeatureID", group, time.col)]) ){
diff --git a/R/miaTime.R b/R/miaTime.R
index 8a554bd..5e1b62f 100755
--- a/R/miaTime.R
+++ b/R/miaTime.R
@@ -3,7 +3,8 @@
 #' \code{miaTime} implements time series methods in \link[mia]{mia} ecosystem.
 #' @name miaTime-package
 #' @aliases miaTime
-#' @seealso \link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment} class
+#' @seealso
+#' \link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment} class
 "_PACKAGE"
 NULL
 
diff --git a/man/getShortTermChange.Rd b/man/getShortTermChange.Rd
index dc06dc8..3ea241e 100644
--- a/man/getShortTermChange.Rd
+++ b/man/getShortTermChange.Rd
@@ -104,5 +104,5 @@ Ji, B.W., et al. (2020) Macroecological dynamics of gut microbiota.
 Nat Microbiol 5, 768–775 . doi: https://doi.org/10.1038/s41564-020-0685-1
 }
 \seealso{
-\code{\link[mia:getStepwiseDivergence]{mia::getStepwiseDivergence()}}
+\code{\link[miaTime:getStepwiseDivergence]{getStepwiseDivergence()}}
 }