Initialise agglomerateByModule

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: mia
 Type: Package
-Version: 1.17.7
+Version: 1.17.8
 Authors@R:
     c(person(given = "Tuomas", family = "Borman", role = c("aut", "cre"),
              email = "tuomas.v.borman@utu.fi",

NAMESPACE

@@ -25,6 +25,7 @@ export(addPrevalence)
 export(addPrevalentAbundance)
 export(addRDA)
 export(addTaxonomyTree)
+export(agglomerateByModule)
 export(agglomerateByPrevalence)
 export(agglomerateByRank)
 export(agglomerateByRanks)
@@ -192,6 +193,7 @@ exportMethods(addPrevalence)
 exportMethods(addPrevalentAbundance)
 exportMethods(addRDA)
 exportMethods(addTaxonomyTree)
+exportMethods(agglomerateByModule)
 exportMethods(agglomerateByPrevalence)
 exportMethods(agglomerateByRank)
 exportMethods(agglomerateByRanks)

NEWS

@@ -178,3 +178,4 @@ Changes in version 1.17.x
 + Added add.dimred parameter to meltSE() (2025-04-29)
 + Added transformation method that replaces values <= threshold with NA (2025-05-09)
 + importMetaPhlAn: Preserve features without taxonomy
++ Added agglomerateByModule()

R/AllGenerics.R

@@ -71,6 +71,11 @@ setGeneric("agglomerateByRank", signature = "x", function(x, ...)
 setGeneric("agglomerateByVariable", signature = "x", function(x, ...)
     standardGeneric("agglomerateByVariable"))
 
+#' @rdname agglomerate-methods
+#' @export
+setGeneric("agglomerateByModule", signature = "x", function(x, ...)
+    standardGeneric("agglomerateByModule"))
+
 #' @rdname calculateDMN
 #' @export
 setGeneric("calculateDMN", signature = c("x"), function(x, ...)

R/agglomerate.R

@@ -12,10 +12,10 @@
 #' with certain taxonomic ranks, as defined in \code{rowData}. Only available
 #' \code{\link{taxonomyRanks}} can be used.
 #'
-#' \code{agglomerateByVariable} merges data on rows or columns of a
-#' \code{SummarizedExperiment} as defined by a \code{factor} alongside the
-#' chosen dimension. This function allows agglomeration of data based on other
-#' variables than taxonomy ranks.
+#' \code{agglomerateByVariable} and \code{agglomerateByModule} merge data on
+#' rows or columns of a \code{SummarizedExperiment} as defined by a
+#' \code{factor} alongside the chosen dimension. This function allows
+#' agglomeration of data based on other variables than taxonomy ranks.
 #' Metadata from the \code{rowData} or \code{colData} are
 #' retained as defined by \code{archetype}.
 #' \code{\link[SummarizedExperiment:SummarizedExperiment-class]{assay}} are
@@ -32,17 +32,21 @@
 #' \code{\link[scuttle:sumCountsAcrossFeatures]{sumCountsAcrossFeatures}}.
 #' However, additional support for \code{TreeSummarizedExperiment} was added and
 #' science field agnostic names were used. In addition the \code{archetype}
-#' argument lets the user select how to preserve row or column data.
+#' argument lets the user select how to preserve row or column data. For merge
+#' data of assays the function from \code{scuttle} are used.
 #'
-#' For merge data of assays the function from \code{scuttle} are used.
+#' \code{agglomerateByModule} allows to agglomerate features or samples based
+#' on one or multiple variables of logical or numeric binary (0/1) type. It is
+#' particularly useful for agglomerating by taxonomic or functional modules,
+#' each defined by a logical or binary variable in the \code{rowData}, as
+#' features can belong to several modules.
 #'
 #' @return
 #' \code{agglomerateByRank} returns a taxonomically-agglomerated,
 #' optionally-pruned object of the same class as \code{x}.
-#' \code{agglomerateByVariable} returns an object of the same class as \code{x}
-#' with the specified entries merged into one entry in all relevant components.
-#' \code{agglomerateByRank} returns a taxonomically-agglomerated,
-#' optionally-pruned object of the same class as \code{x}.
+#' \code{agglomerateByVariable} and \code{agglomerateByModule} return an object
+#' of the same class as \code{x} with the specified entries merged into one
+#' entry in all relevant components.
 #'
 #' @inheritParams getPrevalence
 #'
@@ -120,9 +124,11 @@
 #' \code{factor vector}. A column name from \code{rowData(x)} or
 #' \code{colData(x)} or alternatively a vector specifying how the merging is
 #' performed. If vector, the value must be the same length as
-#' \code{nrow(x)/ncol(x)}. Rows/Cols corresponding to the same level will be
-#' merged. If \code{length(levels(group)) == nrow(x)/ncol(x)}, \code{x} will be
-#' returned unchanged.
+#' \code{nrow(x)/ncol(x)}. Rows or columns corresponding to the same level will
+#' be merged. If \code{length(levels(group)) == nrow(x)/ncol(x)}, \code{x} will
+#' be returned unchanged. For \code{agglomerateByModule}, \code{group} should
+#' specify one or several names of logical or numeric binary variables from the
+#' \code{rowData(x)/colData(x)} by which to agglomerate rows or columns.
 #'
 #' @param f Deprecated. Use \code{group} instead.
 #'
@@ -213,6 +219,35 @@
 #'     GlobalPatterns, by = "cols", colData(GlobalPatterns)$SampleType)
 #' merged
 #'
+#' ## Agglomerate by multiple modules
+#'
+#' # Generate 30 random modules
+#' N_module <- 30L
+#' modules <- sample(
+#'     c(TRUE, FALSE),
+#'     size = nrow(tse) * N_module,
+#'     prob = c(0.2, 0.8),
+#'     replace = TRUE
+#' )
+#'
+#' # Convert modules to matrix
+#' modules <- matrix(modules, nrow = nrow(tse))
+#' 
+#' # Add module names as colnames
+#' colnames(modules) <- paste0("module_", seq_len(ncol(modules)))
+#' 
+#' # Add modules to rowData
+#' rowData(tse) <- cbind(rowData(tse), modules)
+#' 
+#' # Extract module columns
+#' module_columns <- grep("module_", colnames(rowData(tse)), value = TRUE)
+#'
+#' # Agglomerate based on modules
+#' tse_module <- agglomerateByModule(tse, by = 1, group = module_columns)
+#' 
+#' # Optionally, store results into altExp slot
+#' altExp(tse, "modules") <- tse_module
+#'
 #' @seealso
 #' \code{\link[=splitOn]{splitOn}}
 #' \code{\link[=unsplitOn]{unsplitOn}}
@@ -389,6 +424,49 @@ setMethod("agglomerateByVariable", signature = c(x = "SummarizedExperiment"),
     }
 )
 
+#' @rdname agglomerate-methods
+#' @export
+#' @importFrom S4Vectors DataFrame SimpleList
+setMethod("agglomerateByModule", signature = c(x = "SummarizedExperiment"),
+    function(x, by, group, na.rm = FALSE){
+        # Check margin
+        by <- .check_MARGIN(by)
+        # Select side information based on margin
+        FUN <- switch(by, rowData, colData)
+        # Check group
+        if( !all(group %in% names(FUN(x))) ){
+            stop("some elements in 'group' did not match with any column of ",
+                as.character(substitute(FUN)), call. = FALSE)
+        }
+        # Extract modules table
+        modules <- as.matrix(FUN(x)[ , group, drop = FALSE])
+        # Check and process modules
+        modules <- .check_and_process_modules(modules)
+        # Merge assays by module
+        assays <- mapply(.agglomerate_module_assay, assayNames(x), assays(x),
+            MoreArgs = list(by = by, modules = modules, na.rm = na.rm),
+            SIMPLIFY = FALSE)
+        # Convert to SimpleList
+        assays <- assays |> SimpleList()
+        # Construct module-wise experiment
+        x <- do.call(class(x), list(
+            assays = assays,
+            colData = if (by == 1L)
+                colData(x) else DataFrame(row.names = colnames(modules)),
+            rowData = if (by == 1L)
+                DataFrame(row.names = colnames(modules)) else rowData(x),
+            metadata = metadata(x))
+        )
+        # Add new names to agglomerated dimension
+        if( by == 1L ){
+            rownames(x) <- colnames(modules)
+        } else {
+            colnames(x) <- colnames(modules)
+        }
+        return(x)
+    }
+)
+
 ################################ HELP FUNCTIONS ################################
 
 # This functions subset the data so that rows that do not have taxonomy

R/getPrevalence.R

@@ -577,7 +577,7 @@ setMethod("agglomerateByPrevalence", signature = c(x = "SummarizedExperiment"),
         }
         #
         # Check assays that they can be merged safely
-        temp <- mapply(.check_assays_for_merge, assayNames(x), assays(x))
+        mapply(.check_assay_for_merge, assayNames(x), assays(x))
         #
         x <- .merge_features(x, rank, check.assays = FALSE, ...)
         pr <- getPrevalent(x, rank = NULL, ...)

R/merge.R

@@ -8,8 +8,7 @@
     dim.type <- match.arg(dim.type)
     if(!is.character(f) && !is.factor(f)){
         stop("'f' must be a factor or character vector coercible to a ",
-            "meaningful factor.",
-            call. = FALSE)
+            "meaningful factor.", call. = FALSE)
     }
     if(i != length(f)){
         stop("'f' must have the same number of ",dim.type," as 'x'",
@@ -82,7 +81,6 @@
     if( !.is_a_bool(na.rm) ){
         stop("'na.rm' must be TRUE or FALSE.", call. = FALSE)
     }
-    #
     # Get correct functions based on whether we agglomerate rows or cols
     rowData_FUN <- switch(by, rowData, colData)
     nrow_FUN <- switch(by, nrow, ncol)
@@ -102,34 +100,14 @@
     # can control this behavior; it can specify the preserved rows for every
     # group or index.
     archetype <- .norm_archetype(f, archetype)
-
     # Get assays
     assays <- assays(x)
-    # We check whether the assays include values that cannot be summed. For
-    # instance, summing negative values do not make sense.
-    if( check.assays ){
-        temp <- lapply(seq_len(length(assays)), function(i)
-            .check_assays_for_merge(names(assays)[[i]], assays[[i]]))
-    }
-
-    # Transpose if we are merging columns
-    if( by == 2L ){
-        assays <- lapply(assays, function(mat) t(mat))
-    }
-    # Get the aggregation function based on whether user wants to exclude NAs
-    # and if there are any NAs. scuttle::sumCountsAcrossFeatures cannot handle
-    # NAs so if user wants to exclude them, we use own implementation.
-    FUN <- if( na.rm && anyNA(assays[[1]])) .sum_counts_accross_features_na else
-        sumCountsAcrossFeatures
-    # Agglomerate assays
-    assays <- lapply(assays, FUN, average = average, ids = f, BPPARAM = BPPARAM)
-    # Transpose back to original orientation
-    if( by == 2L ){
-        assays <- lapply(assays, function(mat) t(mat))
-    }
+    # Merge assays
+    assays <- mapply(.agglomerate_assay, assayNames(x), assays, MoreArgs = list(
+        ids = f, by = by, na.rm = na.rm, average = average, BPPARAM = BPPARAM,
+        check.assay = check.assays), SIMPLIFY = FALSE)
     # Convert to SimpleList
     assays <- assays |> SimpleList()
-
     # Now we have agglomerated assays, but TreeSE has still the original form.
     # We take specified rows/columns from the TreeSE.
     idx <- .get_element_pos(f, archetype = archetype)
@@ -138,7 +116,6 @@
     } else{
         x <- x[ , idx]
     }
-
     # Add assays back to TreeSE
     assays(x, withDimnames = FALSE) <- assays
     # Change row/colnames. Currently, they have same names as in original data
@@ -165,7 +142,7 @@
 
 # This functions checks if assay has negative or binary values. It does not
 # make sense to sum them, so we give warning to user.
-.check_assays_for_merge <- function(assay.type, assay){
+.check_assay_for_merge <- function(assay.type, assay){
     # Check if assays include binary or negative values
     if( all(assay == 0 | assay == 1) ){
         warning("'", assay.type, "'", " includes binary values.",
@@ -184,6 +161,82 @@
     return(assay)
 }
 
+#' @importFrom DelayedArray DelayedArray type rowsum
+#' @importFrom scuttle sumCountsAcrossFeatures
+.agglomerate_assay <- function(
+    assay.type, assay, by, ids, na.rm, average, BPPARAM, check.assay
+    ){
+    # Check assay
+    if( check.assay ){
+        .check_assay_for_merge(assay.type, assay)
+    }
+    # Transpose if we are merging columns
+    if( by == 2L ){
+        assay <- t(assay)
+    }
+    # Check if NAs are present
+    is_not_na <- !is.na(assay)
+    # Get the aggregation function based on whether user wants to exclude NAs
+    # and if there are any NAs. scuttle::sumCountsAcrossFeatures cannot handle
+    # NAs so if user wants to exclude them, we use own implementation.
+    FUN <- if( na.rm && any(!is_not_na) ) .sum_counts_accross_features_na else
+        sumCountsAcrossFeatures
+    assay <- FUN(assay, ids, average = average, BPPARAM = BPPARAM)
+    # Transpose back to original orientation
+    if( by == 2L ){
+        assay <- t(assay)
+    }
+    return(assay)
+}
+
+# This function checks that modules are in correct format, i.e., there should
+# be module information for each row/column. Columns should represent modules
+# and rows features. Each cell specifies the membership of feature to the
+# specific module.
+.check_and_process_modules <- function(modules){
+    # Determine class type and NAs
+    is_logical <- is.logical(modules)
+    is_num <- all(modules == 0 | modules == 1)
+    is_na <- is.na(modules)
+    # Check validity of module type
+    if( !is_logical && !is_num ){
+        stop("'groups' variables are not binary.", call. = FALSE)
+    }
+    # Convert to numeric binary to Boolean adjacency matrix
+    if( is_num ){
+        modules <- modules != 0
+    }
+    # Replace NAs
+    if( any(is_na) ){
+        warning("NAs were found in 'groups' variables and were removed",
+            "before agglomerating the experiment.", call. = FALSE)
+        # Zero out NA modules
+        modules[is.na(modules)] <- FALSE
+    }
+    return(modules)
+}
+
+# This function agglomerates an abundance table based on modules.
+.agglomerate_module_assay <- function(assay.type, assay, by, modules, na.rm){
+    # Check assay
+    .check_assay_for_merge(assay.type, assay)
+    # Replace NAs with 0
+    if( na.rm ){
+        assay[is.na(assay)] <- 0
+    }
+    # Transpose if we are merging columns
+    if( by == 2L ){
+        assay <- t(assay)
+    }
+    # Compute module-wise assay by cross-product
+    assay <- as.matrix(crossprod(modules, assay))
+    # Transpose back to original orientation
+    if( by == 2L ){
+        assay <- t(assay)
+    }
+    return(assay)
+}
+
 #' @importFrom Biostrings DNAStringSetList
 .merge_refseq_list <- function(sequences_list, f, names, ...){
     threshold <- list(...)[["threshold"]]