Adding difference and division

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: mia
 Type: Package
-Version: 1.17.6
+Version: 1.17.7
 Authors@R:
     c(person(given = "Tuomas", family = "Borman", role = c("aut", "cre"),
              email = "tuomas.v.borman@utu.fi",
@@ -39,7 +39,8 @@ Authors@R:
       person(given = "Pande", family = "Erawijantari", role=c("ctb")),
       person(given = "Danielle", family = "Callan", role=c("ctb")),
       person(given = "Sam", family = "Hillman", role=c("ctb")),
-      person(given = "Jesse", family = "Pasanen", role=c("ctb")))
+      person(given = "Jesse", family = "Pasanen", role=c("ctb")),
+      person(given = "Eetu", family = "Tammi", role=c("ctb")))
 Title: Microbiome analysis
 Description:
     mia implements tools for microbiome analysis based on the

NAMESPACE

@@ -353,6 +353,7 @@ importFrom(MultiAssayExperiment,ExperimentList)
 importFrom(MultiAssayExperiment,MultiAssayExperiment)
 importFrom(MultiAssayExperiment,experiments)
 importFrom(MultiAssayExperiment,sampleMap)
+importFrom(Rcpp,evalCpp)
 importFrom(Rcpp,sourceCpp)
 importFrom(S4Vectors,"metadata<-")
 importFrom(S4Vectors,DataFrame)
@@ -474,3 +475,4 @@ importFrom(vegan,rrarefy)
 importFrom(vegan,scores)
 importFrom(vegan,vegdist)
 useDynLib(mia)
+useDynLib(mia, .registration = TRUE)

R/RcppExports.R

@@ -2,6 +2,10 @@
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 .faith_cpp <- function(assay, rowTree) {
-    .Call('_mia_faith_cpp', PACKAGE = 'mia', assay, rowTree)
+    .Call(`_mia_faith_cpp`, assay, rowTree)
+}
+
+.apply_transformation_difference_or_division <- function(mat, method = "difference") {
+    .Call(`_mia_apply_transformation_difference_or_division`, mat, method)
 }
 

R/transformCounts.R

@@ -51,7 +51,7 @@
 #'   values less than or equal to the threshold are replaced with \code{value}.
 #'   (Default: \code{0})
 #'   \item \code{value}: \code{Numeric scalar}. For  \code{"cutoff"}, specifies
-#'   the replacement value for counts less than or equal to the threshold. 
+#'   the replacement value for counts less than or equal to the threshold.
 #'   (Default: \code{NA})
 #'   \item \code{tree}: \code{phylo}. Phylogeny used in PhILR transformation.
 #'   If \code{NULL}, the tree is retrieved from \code{x}.
@@ -101,6 +101,14 @@
 #' log2 = log2(x)}
 #' where \eqn{x} is a single value of data.
 #'
+#' \item 'difference': Pairwise differences between features.
+#' Calculates \eqn{x - y} for all unique feature pairs across samples,
+#' where \eqn{x} and \eqn{y} are entries of the specified assay.type.
+#'
+#' \item 'division': Pairwise ratios between features.
+#' Calculates \eqn{x / y} for all unique feature pairs across samples,
+#' where \eqn{x} and \eqn{y} are entries of the specified assay.type.
+#'
 #' \item 'pseudocount': Adds only pseudocount.
 #'
 #' \item 'cutoff': In some ecological studies, only strictly positive values
@@ -250,10 +258,10 @@ setMethod("transformAssay", signature = c(x = "SingleCellExperiment"),
 .transform_assay <- function(
         x, assay.type = "counts", assay_name = NULL,
         method = c(
-            "alr", "chi.square", "clr", "css", "cutoff", "frequency",
-            "hellinger", "log", "log10", "log2", "max", "normalize",
-            "pa", "philr", "pseudocount", "range", "rank", "rclr",
-            "relabundance", "rrank", "standardize", "total", "z"),
+            "alr", "chi.square", "clr", "css", "cutoff", "difference", "-",
+            "division", "/", "frequency", "hellinger", "log", "log10", "log2",
+            "max", "normalize", "pa", "philr", "pseudocount", "range", "rank",
+            "rclr", "relabundance", "rrank", "standardize", "total", "z"),
         MARGIN = "samples",
         name = method,
         pseudocount = FALSE,
@@ -278,6 +286,9 @@ setMethod("transformAssay", signature = c(x = "SingleCellExperiment"),
             call. = FALSE)
     }
     method <- match.arg(method, several.ok = FALSE)
+    # Division and difference methods have aliases
+    method <- if( method %in% c("-") ) "difference" else method
+    method <- if( method %in% c("/") ) "division" else method
     # Check that MARGIN is 1 or 2
     MARGIN <- .check_MARGIN(MARGIN)
     # Check pseudocount
@@ -302,7 +313,7 @@ setMethod("transformAssay", signature = c(x = "SingleCellExperiment"),
     attr(assay, "pseudocount") <- NULL
     # Calls help function that does the transformation
     # Help function is different for mia and vegan transformations
-    if( method %in% c("log10", "log2", "css") ){
+    if( method %in% c("log10", "log2", "css", "difference", "division") ){
         transformed_table <- .apply_transformation(
             assay, method, MARGIN, ...)
     } else if( method %in% c("philr") ){
@@ -338,19 +349,21 @@ setMethod("transformAssay", signature = c(x = "SingleCellExperiment"),
         method,
         log10 = .calc_log,
         log2 = .calc_log,
-        css = .calc_css
+        css = .calc_css,
+        difference = .apply_transformation_difference_or_division,
+        division = .apply_transformation_difference_or_division
     )
     # Get transformed table
-    transformed_table <- do.call(
+    assay <- do.call(
         FUN, list(mat = assay, method = method, MARGIN = MARGIN, ...) )
     # Transpose back to normal if MARGIN is row
     if( MARGIN == 1L ){
-        transformed_table <- t(transformed_table)
+        assay <- t(assay)
     }
     # Add method and margin to attributes
-    attr(transformed_table, "mia") <- method
-    attr(transformed_table, "parameters")$margin <- MARGIN
-    return(transformed_table)
+    attr(assay, "mia") <- method
+    attr(assay, "parameters")$margin <- MARGIN
+    return(assay)
 }
 
 ########################.apply_transformation_from_vegan########################
@@ -782,11 +795,42 @@ setMethod("transformAssay", signature = c(x = "SingleCellExperiment"),
     return(mat)
 }
 
+################# .apply_transformation_difference_or_division #################
+# Computes all pairwise differences (x - y) or ration (x / y) between features
+# across samples. Returns a sparse matrix with one row per feature pair.
+# Uses C++ to improve performance on large input matrices.
+#' @useDynLib mia, .registration = TRUE
+#' @importFrom Rcpp evalCpp
+NULL
+#' @keywords internal
+.apply_transformation_difference_or_division <- function(
+        mat, method, MARGIN, ...){
+    # To harmonize transformations, the matrix contains features as columns.
+    # This orientation is used in other transformations such as in vegan.
+    if( ncol(mat) > 1000L ){
+        warning("The input matrix has over 1000 features, which may cause ",
+                "performance issues.", call. = FALSE)
+    }
+    if( is.null(colnames(mat)) ){
+        warning("No names found in the matrix. Generated labels.",
+                call. = FALSE)
+        colnames(mat) <- ncol(mat) |> seq_len() |> as.character()
+    }
+    if( method %in% "division" && any(!is.na(mat) & mat == 0 ) ){
+        warning("The assay contains zero values. Consider adding pseudocount.",
+                call. = FALSE)
+    }
+    # Calculate pairwise difference or division
+    mat <- .Call(
+        `_mia_apply_transformation_difference_or_division`, mat, method)
+    return(mat)
+}
+
 # This function is used to add transformed table back to TreeSE. With most of
-# the methods it is simple: it is added to assay. However, with philr, the
-# features do not match with original ones, so we add philr-transformed data
-# to altExp. If philr was, applied to columns, we cannot use altExp so
-# we return only the transformed data.
+# the methods it is simple: it is added to assay. However, with transformations
+# that change the dimensionality (e.g. philr, difference, division), the
+# transformed table is added to altExp instead. If transformation was, applied
+# to columns, we cannot use altExp so we return only the transformed data.
 #' @importFrom stats setNames
 .add_transformed_data <- function(x, mat, name){
     rnames_ok <- nrow(x) == nrow(mat)
@@ -825,4 +869,4 @@ setMethod("transformAssay", signature = c(x = "SingleCellExperiment"),
     }
     mat[ mat <= threshold ] <- value
     return(mat)
-}
+}

src/.gitignore

@@ -1,18 +1,2 @@
-api.o
-api_s.o
-biom.o
-biom_s.o
-propstack.o
-tree.o
-tree_s.o
-tse.o
-unifrac.o
-unifrac_internal.o
-unifrac_internal_s.o
-unifrac_s.o
-RcppExports.o
-assay.o
-faith_R.o
-propmap.o
-mia.dll
-mia.so
+*.o
+*.so

src/RcppExports.cpp

@@ -22,9 +22,22 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// apply_transformation_difference_or_division
+S4 apply_transformation_difference_or_division(NumericMatrix mat, std::string method);
+RcppExport SEXP _mia_apply_transformation_difference_or_division(SEXP matSEXP, SEXP methodSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type mat(matSEXP);
+    Rcpp::traits::input_parameter< std::string >::type method(methodSEXP);
+    rcpp_result_gen = Rcpp::wrap(apply_transformation_difference_or_division(mat, method));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_mia_faith_cpp", (DL_FUNC) &_mia_faith_cpp, 2},
+    {"_mia_apply_transformation_difference_or_division", (DL_FUNC) &_mia_apply_transformation_difference_or_division, 2},
     {NULL, NULL, 0}
 };
 

src/transformCounts.cpp

@@ -0,0 +1,84 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// Helper for difference
+void calculate_difference(
+        const NumericMatrix& mat, int i, int j, int k, int pair_idx,
+        std::vector<int>& i_vec, std::vector<int>& j_vec,
+        std::vector<double>& x_vec) {
+    double diff = mat(k, i) - mat(k, j);
+    if( diff != 0.0 ){
+        i_vec.push_back(k);
+        j_vec.push_back(pair_idx);
+        x_vec.push_back(diff);
+    }
+}
+
+// Helper for division
+void calculate_division(
+        const NumericMatrix& mat, int i, int j, int k, int pair_idx,
+        std::vector<int>& i_vec, std::vector<int>& j_vec,
+        std::vector<double>& x_vec) {
+    double denom = mat(k, j);
+    if( denom != 0.0 && std::isfinite(denom) ){
+        double ratio = mat(k, i) / denom;
+        if( ratio != 0.0 && std::isfinite(ratio) ){
+            i_vec.push_back(k);
+            j_vec.push_back(pair_idx);
+            x_vec.push_back(ratio);
+        }
+    }
+}
+
+// [[Rcpp::export(name = ".apply_transformation_difference_or_division")]]
+S4 apply_transformation_difference_or_division(
+        NumericMatrix mat, std::string method = "difference") {
+    // samples = rows, features = cols
+    int n_samples  = mat.nrow();
+    int n_features = mat.ncol();
+    int n_pairs    = n_features * (n_features - 1) / 2;
+    CharacterVector original_feature_names = colnames(mat);
+
+    // Prepare output vectors
+    std::vector<int> i_vec, j_vec;
+    std::vector<double> x_vec;
+    CharacterVector colnames(n_pairs);
+
+    int pair_idx = 0;
+    for( int i = 0; i < n_features - 1; ++i ){
+        for( int j = i + 1; j < n_features; ++j ){
+            for( int k = 0; k < n_samples; ++k ){
+                if( method == "difference" ){
+                    calculate_difference(
+                        mat, i, j, k, pair_idx, i_vec, j_vec, x_vec);
+                } else if( method == "division" ){
+                    calculate_division(
+                        mat, i, j, k, pair_idx, i_vec, j_vec, x_vec);
+                } else{
+                    stop("Unknown method: use 'difference' or 'division'");
+                }
+            }
+            // Generate colname for feature pair
+            std::string name_i = as<std::string>(original_feature_names[i]);
+            std::string name_j = as<std::string>(original_feature_names[j]);
+            colnames[pair_idx] = name_i + "-" + name_j;
+            ++pair_idx;
+        }
+    }
+
+    // Create dgTMatrix (triplet sparse matrix)
+    S4 tmat("dgTMatrix");
+    // Samples as rows
+    tmat.slot("i") = IntegerVector(i_vec.begin(), i_vec.end());
+    // Pairs as columns
+    tmat.slot("j") = IntegerVector(j_vec.begin(), j_vec.end());
+    tmat.slot("x") = NumericVector(x_vec.begin(), x_vec.end());
+    tmat.slot("Dim") = IntegerVector::create(n_samples, n_pairs);
+    tmat.slot("Dimnames") = List::create(rownames(mat), colnames);
+
+    // Convert dgTMatrix to dgCMatrix
+    Function as("as");
+    S4 cmat = as(tmat, "CsparseMatrix");
+
+    return cmat;
+}