Update tree agglomeration (#656)

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: mia
 Type: Package
-Version: 1.15.3
+Version: 1.15.4
 Authors@R:
     c(person(given = "Tuomas", family = "Borman", role = c("aut", "cre"),
              email = "[email protected]",

NAMESPACE

@@ -347,6 +347,7 @@ importFrom(SummarizedExperiment,colData)
 importFrom(SummarizedExperiment,rowData)
 importFrom(SummarizedExperiment,rowRanges)
 importFrom(TreeSummarizedExperiment,changeTree)
+importFrom(TreeSummarizedExperiment,subsetByLeaf)
 importFrom(ape,as.phylo)
 importFrom(ape,bind.tree)
 importFrom(ape,collapse.singles)

NEWS

@@ -156,3 +156,6 @@ computation
 + transformAssay can apply transformation to altExp
 + Added CSS transformation
 + In agglomerateByVariable, splitOn and getDominant, use 'group' to specify grouping variable.
+
+Changes in version 1.15.x
++ subsetBy*: added update.tree argument

R/agglomerate.R

@@ -23,9 +23,6 @@
 #'   regarded as empty. (Default: \code{c(NA, "", " ", "\t")}). They will be
 #'   removed if \code{na.rm = TRUE} before agglomeration.
 #'
-#' @param update.tree \code{Logical scalar}. Should
-#'   \code{rowTree()} also be agglomerated? (Default: \code{FALSE})
-#'
 #' @param agglomerateTree Deprecated. Use \code{update.tree} instead.
 #' 
 #' @param agglomerate.tree Deprecated. Use \code{update.tree} instead.
@@ -450,16 +447,18 @@ setMethod(
 
 # Agglomerate all rowTrees found in TreeSE object. Get tips that represent
 # rows and remove all others.
+#' @importFrom TreeSummarizedExperiment subsetByLeaf
 .agglomerate_trees <- function(x, by = 1, ...){
     # Get right functions based on direction
     tree_names_FUN <- switch(
         by, "1" = rowTreeNames, "2" = colTreeNames, stop("."))
     links_FUN <- switch(by, "1" = rowLinks, "2" = colLinks, stop("."))
     tree_FUN <- switch(by, "1" = rowTree, "2" = colTree, stop("."))
-    # Get right argument names for changeTree call
+    # Get right argument names for subsetByLeaf call
     args_names <- switch(
-        by, "1" = c("x", "rowTree", "rowNodeLab", "whichRowTree"),
-        "2" = c("x", "colTree", "colNodeLab", "whichColTree"),
+        by,
+        "1" = c("x", "rowLeaf", "whichRowTree", "updateTree"),
+        "2" = c("x", "colLeaf", "whichColTree", "updateTree"),
         stop("."))
     # Get names of trees and links between trees and rows
     tree_names <- tree_names_FUN(x)
@@ -475,11 +474,9 @@ setMethod(
             # Get names of nodes that are preserved
             links_temp <- links_temp[["nodeLab"]]
             # Agglomerate the tree
-            tree <- .prune_tree(tree, links_temp, ...)
-            # Change the tree with agglomerated version
-            args <- list(x, tree, links_temp, name)
+            args <- list(x, links_temp, name, TRUE)
             names(args) <- args_names
-            x <- do.call(changeTree, args)
+            x <- do.call(subsetByLeaf, args)
         }
     }
     return(x)
@@ -507,12 +504,22 @@ setMethod(
         # even after pruning; these rows have still child-nodes that represent
         # other rows.
         # Suppress warning: drop all tips of the tree: returning NULL
-        suppressWarnings(
-            tree <- drop.tip(
+        tree <- tryCatch({
+            drop.tip(
                 tree, remove_tips,
                 trim.internal = FALSE,
                 collapse.singles = FALSE)
-        )
+        }, warning = function(w) {
+            # Do nothing on warning
+        }, error = function(e) {
+            # Try to prune by also pruning internal nodes. Sometimes that is the
+            # case; we need to trim also internal nodes in order to prune
+            # leaf.
+            drop.tip(
+                tree, remove_tips,
+                trim.internal = TRUE,
+                collapse.singles = FALSE)
+        })
         # If all tips were dropped, the result is NULL --> stop loop
         if( is.null(tree) ){
             warning("Pruning resulted to empty tree.", call. = FALSE)

R/getPrevalence.R

@@ -28,6 +28,9 @@
 #'
 #' @param na.rm \code{Logical scalar}. Should NA values be omitted when calculating
 #' prevalence? (Default: \code{TRUE})
+#' 
+#' @param update.tree \code{Logical scalar}. Should
+#' \code{rowTree()} also be agglomerated? (Default: \code{FALSE})
 #'
 #' @param ... additional arguments
 #' \itemize{
@@ -444,6 +447,24 @@ setMethod("subsetByPrevalent", signature = c(x = "SummarizedExperiment"),
     }
 )
 
+#' @rdname getPrevalence
+#' @export
+setMethod("subsetByPrevalent", signature = c(x = "TreeSummarizedExperiment"),
+    function(x, update.tree = FALSE, ...){
+        # Check that update.tree is logical value
+        if( !.is_a_bool(update.tree) ){
+            stop("'update.tree' must be TRUE or FALSE.", call. = FALSE)
+        }
+        #
+        x <- callNextMethod(x, ...)
+        # Agglomerate tree if specified
+        if( update.tree ){
+            x <- .agglomerate_trees(x, ...)
+        }
+    return(x)
+    }
+)
+
 ############################# subsetByRare #################################
 
 #' @rdname getPrevalence
@@ -462,6 +483,24 @@ setMethod("subsetByRare", signature = c(x = "SummarizedExperiment"),
     }
 )
 
+#' @rdname getPrevalence
+#' @export
+setMethod("subsetByRare", signature = c(x = "TreeSummarizedExperiment"),
+    function(x, update.tree = FALSE, ...){
+        # Check that update.tree is logical value
+        if( !.is_a_bool(update.tree) ){
+            stop("'update.tree' must be TRUE or FALSE.", call. = FALSE)
+        }
+        #
+        x <- callNextMethod(x, ...)
+        # Agglomerate tree if specified
+        if( update.tree ){
+            x <- .agglomerate_trees(x, ...)
+        }
+        return(x)
+    }
+)
+
 ############################# getPrevalentAbundance ############################
 
 #' @rdname getPrevalence

tests/testthat/test-5prevalence.R

@@ -310,6 +310,13 @@ test_that("subsetByPrevalent", {
     alias <- subsetByPrevalent(gp_null, detection=5, prevalence = 0.33, rank = "Phylum")
     alias <- unname(assay(alias, "counts"))
     expect_equal(alias, pr2)
+
+    # Check that tree subsetting works
+    expect_error(subsetByPrevalent(GlobalPatterns, update.tree = 1))
+    expect_error(subsetByPrevalent(GlobalPatterns, update.tree = NULL))
+    expect_error(subsetByPrevalent(GlobalPatterns, update.tree = c(TRUE, FALSE)))
+    tse_sub <- subsetByPrevalent(GlobalPatterns, prevalence = 0.4, rank = "Genus", update.tree = TRUE)
+    expect_equal(length(rowTree(tse_sub)$tip.label), nrow(tse_sub))
 
 })
 
@@ -371,7 +378,13 @@ test_that("subsetByRare", {
     alias <- subsetByRare(gp_null, detection=5, prevalence = 0.33, rank = "Phylum")
     alias <- unname(assay(alias, "counts"))
     expect_equal(alias, pr2)
-
+
+    # Check that tree subsetting works
+    expect_error(subsetByRare(GlobalPatterns, update.tree = 1))
+    expect_error(subsetByRare(GlobalPatterns, update.tree = NULL))
+    expect_error(subsetByRare(GlobalPatterns, update.tree = c(TRUE, FALSE)))
+    tse_sub <- subsetByRare(GlobalPatterns, rank = "Class", update.tree = TRUE)
+    expect_equal(length(rowTree(tse_sub)$tip.label), nrow(tse_sub))
 })
 
 test_that("agglomerateByPrevalence", {