New dissimilarity functions

NAMESPACE

@@ -7,6 +7,7 @@ export(addAlpha)
 export(addCluster)
 export(addContaminantQC)
 export(addDPCoA)
+export(addDissimilarity)
 export(addDivergence)
 export(addDominant)
 export(addLDA)
@@ -25,11 +26,8 @@ export(calculateCCA)
 export(calculateDMN)
 export(calculateDMNgroup)
 export(calculateDPCoA)
-export(calculateJSD)
 export(calculateNMDS)
-export(calculateOverlap)
 export(calculateRDA)
-export(calculateUnifrac)
 export(cluster)
 export(convertFromBIOM)
 export(convertFromDADA2)
@@ -49,6 +47,7 @@ export(getBestDMNFit)
 export(getCrossAssociation)
 export(getDMN)
 export(getDPCoA)
+export(getDissimilarity)
 export(getDominant)
 export(getExperimentCrossAssociation)
 export(getExperimentCrossCorrelation)
@@ -113,11 +112,8 @@ export(right_join)
 export(runCCA)
 export(runDMN)
 export(runDPCoA)
-export(runJSD)
 export(runNMDS)
-export(runOverlap)
 export(runRDA)
-export(runUnifrac)
 export(setTaxonomyRanks)
 export(splitByRanks)
 export(splitOn)
@@ -147,6 +143,7 @@ exportMethods(addAlpha)
 exportMethods(addCCA)
 exportMethods(addCluster)
 exportMethods(addContaminantQC)
+exportMethods(addDissimilarity)
 exportMethods(addDivergence)
 exportMethods(addDominant)
 exportMethods(addHierarchyTree)
@@ -164,9 +161,6 @@ exportMethods(agglomerateByVariable)
 exportMethods(bestDMNFit)
 exportMethods(calculateDMN)
 exportMethods(calculateDMNgroup)
-exportMethods(calculateJSD)
-exportMethods(calculateOverlap)
-exportMethods(calculateUnifrac)
 exportMethods(checkTaxonomy)
 exportMethods(cluster)
 exportMethods(convertToBIOM)
@@ -185,6 +179,7 @@ exportMethods(getCCA)
 exportMethods(getCrossAssociation)
 exportMethods(getDMN)
 exportMethods(getDPCoA)
+exportMethods(getDissimilarity)
 exportMethods(getDominant)
 exportMethods(getExperimentCrossAssociation)
 exportMethods(getExperimentCrossCorrelation)
@@ -232,7 +227,6 @@ exportMethods(rarefyAssay)
 exportMethods(relAbundanceCounts)
 exportMethods(relabundance)
 exportMethods(right_join)
-exportMethods(runOverlap)
 exportMethods(splitOn)
 exportMethods(subsampleCounts)
 exportMethods(subsetByPrevalent)

R/calculateDistance.R

@@ -1,6 +1,201 @@
-# calculateDistance function is removed. This function is left for other functions 
-# to use
-.calculate_distance <- function(mat, FUN = stats::dist, ...){
-    # Distance between all samples against all samples
-    do.call(FUN, c(list(mat),list(...)))
+#' Calculate dissimilarities
+#'
+#' These functions calculate dissimilarities on data stored in a 
+#'  \code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment}}
+#'  object.
+#'
+#' @param x a \code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment}}
+#'   object.
+#'
+#' @param method \code{Character scalar}. Specifies which distance to calculate.
+#' 
+#' @param name \code{Character scalar}. The name to be used to store the result 
+#'  in the reducedDims of the output. (Default: \code{method})
+#'
+#' @param assay.type \code{Character scalar}. Specifies which assay to use for 
+#'   calculation. (Default: \code{"counts"})
+#'
+#' @param assay_name \code{Character scalar}. Specifies which assay to use for 
+#'   calculation. (Please use \code{assay.type} instead. At 
+#'   some point \code{assay_name} will be disabled.)
+#'
+#' @param tree.name \code{Character scalar}. Specifies which tree will be used in 
+#'   calculation. (Default: \code{"phylo"})
+#'   
+#' @param transposed \code{Logical scalar}. Specifies if x is transposed with cells in
+#'   rows. (Default: \code{FALSE})
+#'
+#' @param ... other arguments passed onto \code{\link[vegan:vegdist]{vegdist}}
+#'
+#' @return 
+#' \code{getDissimilarity} returns a distance matrix.
+#' 
+#' \code{addDissimilarity} returns a
+#'   \code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment}}
+#'   with distance matrix added to reducedDim slot.
+#'
+#' @name getDissimilarity
+#'
+#' @export
+#'
+#' @examples
+#' data(GlobalPatterns)
+#' tse <- GlobalPatterns
+#' tse <- addDissimilarity(tse, method = "overlap", detection = 0.25)
+#' reducedDim(tse, "overlap")
+#' 
+NULL
+
+#' @rdname getDissimilarity
+#' @export
+setGeneric(
+    "getDissimilarity", signature = c("x"), function(x, method, ...)
+        standardGeneric("getDissimilarity"))
+
+#' @rdname getDissimilarity
+#' @export
+setMethod(
+    "getDissimilarity", signature = c(x = "TreeSummarizedExperiment"),
+    function(
+        x, method, assay_name = "counts", assay.type = assay_name,
+          tree.name = "phylo", transposed = FALSE, ...){
+    mat <- assay(x, assay.type)
+    if(!transposed){
+        mat <- t(mat)
+    }
+    tree <- rowTree(x, tree.name)
+    res <- getDissimilarity(mat, method = method, tree = tree, ...)
+    return(res)
+    }
+)
+
+#' @rdname getDissimilarity
+#' @export
+setMethod(
+    "getDissimilarity", signature = c(x = "SummarizedExperiment"),
+    function(
+        x, method, assay_name = "counts", assay.type = assay_name, 
+          transposed = FALSE, ...){
+    browser()
+    # Input checks
+    if(!.is_non_empty_string(assay.type)){
+        stop("'assay.type' must be a non-empty single character value",
+            call. = FALSE)
+    }
+    if(!.is_non_empty_string(method)){
+        stop("'method' must be a non-empty single character value",
+            call. = FALSE)
+    }
+    if(!.is_a_bool(transposed)){
+        stop("'na.rm' must be TRUE or FALSE.", call. = FALSE)
+    }
+    mat <- assay(x, assay.type)
+    if(!transposed){
+        mat <- t(mat)
+    }
+    if( method == "unifrac"){
+      node.label <- .get_rowlinks(x, assay.type)
+      mat <- .calculate_dissimilarity(mat = x, method = method, 
+                                      node.label = node.label, ...)
+    }
+    else{
+      mat <- .calculate_dissimilarity(mat = x, method = method, ...)
+    }
+    return(mat)
+    }
+)
+
+#' @rdname getDissimilarity
+#' @export
+setMethod(
+    "getDissimilarity", signature = c(x = "ANY"),
+    function(
+        x, method, assay_name = "counts", assay.type = assay_name, ...){
+    # Input check
+    if( !.is_a_string(method) ){
+        stop("'method' must be a single character value.", call. = FALSE)
+    }
+    #
+    # Calculate dissimilarity
+    mat <- .calculate_dissimilarity(mat = x, method = method, ...)
+    return(mat)
+    }
+)
+
+#' @rdname getDissimilarity
+#' @export
+setGeneric(
+  "addDissimilarity", signature = c("x"), function(x, method, ...)
+    standardGeneric("addDissimilarity"))
+
+#' @rdname getDissimilarity
+#' @export
+setMethod(
+  "addDissimilarity", signature = c(x = "SummarizedExperiment"),
+  function(
+    x, method, assay_name = "counts", assay.type = assay_name, name = method,
+      transposed = FALSE, ...){
+    res <- getDissimilarity(x, method = method, assay.type = assay.type, 
+                            transposed = transposed, ...)
+    # Input checks
+    if(!.is_non_empty_string(assay.type)){
+      stop("'assay.type' must be a non-empty single character value",
+           call. = FALSE)
+    }
+    if(!.is_non_empty_string(method)){
+      stop("'method' must be a non-empty single character value",
+           call. = FALSE)
+    }
+    if(!.is_non_empty_string(name)){
+      stop("'name' must be a non-empty single character value",
+           call. = FALSE)
+    }
+    if(!.is_a_bool(transposed)){
+      stop("'na.rm' must be TRUE or FALSE.", call. = FALSE)
+    }
+    if( !identical(rownames(as.matrix(res)), colnames(assay(x, assay.type))) ){
+      warning("Samples of the dissimilarity matrix should be the same as the
+            samples in columns of the assay specified with 'assay.type'. The 
+            result is not added to reducedDim.")
+      return(res)
+    }
+    else{
+      .add_values_to_reducedDims(x, as.matrix(res), name)
+    }
+
+  }
+)
+
+.calculate_dissimilarity <- function(
+        mat, method, node.label = NULL, diss.fun = NULL, tree = NULL, ...){
+    browser()
+    # input check
+    if( !(is.null(diss.fun) || is.function(diss.fun)) ){
+        stop("'diss.fun' must be NULL or a function.", call. = FALSE)
+    }
+    #
+    args <- c(list(mat, method = method, node.label = node.label, list(...)))
+    # If the dissimilarity functon is not specified, get default choice
+    if( is.null(diss.fun) ){
+        if( method %in% c("overlap") ){
+            diss.fun <- .get_overlap
+            message("'diss.fun' defaults to .get_overlap.")
+        } else if( method %in% c("unifrac")  ){
+            args[["tree"]] <- tree
+            diss.fun <- .get_unifrac
+            message("'diss.fun' defaults to .get_unifrac.")
+        } else if( method %in% c("jsd")  ){
+            diss.fun <- .get_jsd
+            message("'diss.fun' defaults to .get_jsd.")
+        } else if( requireNamespace("vegan") ){
+            diss.fun <- vegan::vegdist
+            message("'diss.fun' defaults to vegan::vegdist.")
+        } else{
+            diss.fun <- stats::dist
+            message("'diss.fun' defaults to stats::dist.")
+        }
+    }
+    # Calculate dissimilarity with specified function
+    res <- do.call(diss.fun, args)
+    return(res)
 }

R/calculateJSD.R

@@ -1,97 +1,3 @@
-#' Calculate the Jensen-Shannon Divergence
-#'
-#' This function calculates the Jensen-Shannon Divergence (JSD) in a
-#' \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}
-#' object.
-#' 
-#' @param x a numeric matrix with samples as rows or a
-#'   \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}}
-#'   object.
-#'
-#' @param assay.type \code{Character scalar}. Specifies the name of the
-#'   assay used in calculation. (Default: \code{"counts"})
-#'   
-#' @param exprs_values Deprecated. Use \code{assay.type} instead.
-#'   
-#' @param assay_name Deprecated. Use \code{assay.type} instead.
-#'
-#' @param chunkSize \code{Integer scalar}. Defines the size of data send
-#'   to the individual worker. Only has an effect, if \code{BPPARAM} defines
-#'   more than one worker. (Default: \code{nrow(x)})
-#' 
-#' @param BPPARAM A
-#'   \code{\link[BiocParallel:BiocParallelParam-class]{BiocParallelParam}}
-#'   object specifying whether the calculation should be parallelized.
-#'
-#' @param transposed \code{Logical scalar}. Is x transposed with samples in rows?
-#' (Default: \code{FALSE})
-#'
-#' @param ... optional arguments not used.
-#'
-#' @return a sample-by-sample distance matrix, suitable for NMDS, etc.
-#'
-#' @seealso
-#' \url{http://en.wikipedia.org/wiki/Jensen-Shannon_divergence}
-#'
-#' @references
-#' Jensen-Shannon Divergence and Hilbert space embedding.
-#' Bent Fuglede and Flemming Topsoe University of Copenhagen,
-#' Department of Mathematics
-#' \url{http://www.math.ku.dk/~topsoe/ISIT2004JSD.pdf}
-#'
-#' @name calculateJSD
-#'
-#' @author
-#' Susan Holmes \email{susan@@stat.stanford.edu}.
-#' Adapted for phyloseq by Paul J. McMurdie.
-#' Adapted for mia by Felix G.M. Ernst
-#'
-#' @export
-#'
-#' @examples
-#' data(enterotype)
-#' library(scater)
-#'
-#'
-#' jsd <- calculateJSD(enterotype)
-#' class(jsd)
-#' head(jsd)
-#'
-#' enterotype <- runMDS(enterotype, FUN = calculateJSD, name = "JSD",
-#'                      exprs_values = "counts")
-#' head(reducedDim(enterotype))
-#' head(attr(reducedDim(enterotype),"eig"))
-#' attr(reducedDim(enterotype),"GOF")
-NULL
-
-setGeneric("calculateJSD", signature = c("x"),
-           function(x, ...)
-             standardGeneric("calculateJSD"))
-
-#' @rdname calculateJSD
-#' @export
-setMethod("calculateJSD", signature = c(x = "ANY"),
-    function(x, ...){
-        .calculate_distance(x, FUN = runJSD, ...)
-    }
-)
-
-#' @rdname calculateJSD
-#'
-#' @importFrom SummarizedExperiment assay
-#'
-#' @export
-setMethod("calculateJSD", signature = c(x = "SummarizedExperiment"),
-    function(x, assay.type = assay_name, assay_name = exprs_values, 
-             exprs_values = "counts", transposed = FALSE, ...){
-        mat <- assay(x, assay.type)
-        if(!transposed){
-            mat <- t(mat)
-        }
-        calculateJSD(mat, ...)
-    }
-)
-
 # written by Susan Holmes \email{susan@@stat.stanford.edu}.
 # Adapted for phyloseq by Paul J. McMurdie.
 # Adapted for mia by Felix G.M. Ernst
@@ -113,15 +19,12 @@ setMethod("calculateJSD", signature = c(x = "SummarizedExperiment"),
     return(rowSums(d, na.rm = TRUE))
 }
 
-#' @rdname calculateJSD
-#'
 #' @importFrom utils combn
 #' @importFrom stats as.dist
 #' @importFrom BiocParallel SerialParam register bplapply bpisup bpstart bpstop
 #' @importFrom DelayedArray getAutoBPPARAM setAutoBPPARAM
-#'
-#' @export
-runJSD <- function(x, BPPARAM = SerialParam(), chunkSize = nrow(x)){
+#' 
+.get_jsd <- function(x, BPPARAM = SerialParam(), chunkSize = nrow(x), ...){
     # input check
     if(is.null(rownames(x))){
         rownames(x) <- seq_len(nrow(x))