R/celda_functions.R

In celda: CEllular Latent Dirichlet Allocation

.sampleLl <- function(llProbs) {
  probsSub <- exp(llProbs - max(llProbs))
  probsNorm <- probsSub / sum(probsSub)
  probsSelect <- sample.int(
    length(probsNorm),
    size = 1L,
    replace = TRUE,
    prob = probsNorm
  )
  return(probsSelect)
}


.cosineDist <- function(x) {
  x <- t(x)
  y <- (1 - .cosine(x)) / 2
  return(stats::as.dist(y))
}


.cosine <- function(x) {
  y <- x %*% t(x) / (sqrt(rowSums(x^2) %*% t(rowSums(x^2))))
  return(y)
}


.spearmanDist <- function(x) {
  y <- (1 - stats::cor(x, method = "spearman")) / 2
  return(stats::as.dist(y))
}


.hellingerDist <- function(x) {
  y <- stats::dist(t(sqrt(x)), method = "euclidean") * 1 / sqrt(2)
  return(y)
}


.normalizeLogProbs <- function(llProbs) {
  llProbs <- exp(sweep(llProbs, 1, base::apply(llProbs, 1, max), "-"))
  probs <- sweep(llProbs, 1, rowSums(llProbs), "/")
  return(probs)
}


#' @title Normalization of count data
#' @description Performs normalization, transformation, and/or scaling of a
#'  counts matrix
#' @param counts Integer matrix. Rows represent features and columns represent
#'  cells.
#' @param normalize Character.
#'  Divides counts by the library sizes for each cell. One of 'proportion',
#'  'cpm', 'median', or 'mean'. 'proportion' uses the total counts for each
#'  cell as the library size. 'cpm' divides the library size of each cell by
#'  one million to produce counts per million. 'median' divides the library
#'  size of each cell by the median library size across all cells. 'mean'
#'  divides the library size of each cell by the mean library size across all
#'  cells.
#' @param scaleFactor Numeric. Sets the scale factor for cell-level
#'  normalization. This scale factor is multiplied to each cell after the
#'  library size of each cell had been adjusted in \code{normalize}. Default
#'  \code{NULL} which means no scale factor is applied.
#' @param transformationFun Function. Applys a transformation such as
#'  \link{sqrt}, \link{log}, \link{log2}, \link{log10}, or \link{log1p}.
#'  If NULL, no transformation will be applied. Occurs after normalization.
#'  Default NULL.
#' @param scaleFun Function. Scales the rows of the normalized and transformed
#'  count matrix. For example, 'scale' can be used to z-score normalize the
#'  rows. Default NULL.
#' @param pseudocountNormalize Numeric. Add a pseudocount to counts before
#'  normalization. Default 0.
#' @param pseudocountTransform Numeric. Add a pseudocount to normalized counts
#'  before applying the transformation function. Adding a pseudocount
#'  can be useful before applying a log transformation. Default  0.
#' @return Numeric Matrix. A normalized matrix.
#' @examples
#' data(celdaCGSim)
#' normalizedCounts <- normalizeCounts(celdaCGSim$counts, "proportion",
#'   pseudocountNormalize = 1)
#' @export
normalizeCounts <- function(counts,
                            normalize = c("proportion", "cpm",
                              "median", "mean"),
                            scaleFactor = NULL,
                            transformationFun = NULL,
                            scaleFun = NULL,
                            pseudocountNormalize = 0,
                            pseudocountTransform = 0) {

  normalize <- match.arg(normalize)
  counts <- as.matrix(counts)

  if (!is.null(transformationFun) &&
    !is.function(transformationFun)) {
    stop("'transformationFun' needs to be of class 'function'")
  }
  if (!is.null(scaleFun) && !is.function(scaleFun)) {
    stop("'scaleFun' needs to be of class 'function'")
  }
  # Perform normalization
  if (normalize == "proportion") {
    norm <- fastNormProp(counts, pseudocountNormalize)
  } else {
    counts <- counts + pseudocountNormalize
    cs <- .colSums(counts, nrow(counts), ncol(counts))
    norm <- switch(
      normalize,
      "proportion" = sweep(counts, 2, cs, "/"),
      "cpm" = sweep(counts, 2, cs / 1e6, "/"),
      "median" = sweep(counts, 2, cs / stats::median(cs), "/"),
      "mean" = sweep(counts, 2, cs / mean(cs), "/")
    )
  }

  if (!is.null(scaleFactor)) {
      norm <- norm * scaleFactor
  }

  if (!is.null(transformationFun)) {
    norm <- do.call(
      transformationFun,
      list(norm + pseudocountTransform)
    )
  }
  if (!is.null(scaleFun)) {
    norm <- t(base::apply(norm, 1, scaleFun))
  }
  colnames(norm) <- colnames(counts)
  rownames(norm) <- rownames(counts)
  return(norm)
}


#' @title Recode cell cluster labels
#' @description Recode cell subpopulaton clusters using a mapping in the
#'  \code{from} and \code{to} arguments.
#' @param sce \linkS4class{SingleCellExperiment} object returned from
#'  \link{celda_C} or \link{celda_CG}. Must contain column
#'  \code{celda_cell_cluster} in
#'  \code{\link{colData}(altExp(sce, altExpName))}.
#' @param from Numeric vector. Unique values in the range of
#'  \code{seq(celdaClusters(sce, altExpName = altExpName))} that correspond to
#'  the original cluster
#'  labels in \code{sce}.
#' @param to Numeric vector. Unique values in the range of
#'  \code{seq(celdaClusters(sce, altExpName = altExpName))} that correspond to
#'  the new cluster labels.
#' @param altExpName The name for the \link{altExp} slot
#'  to use. Default "featureSubset".
#' @return \linkS4class{SingleCellExperiment} object with recoded cell
#'  cluster labels.
#' @examples
#' data(sceCeldaCG)
#' sceReorderedZ <- recodeClusterZ(sceCeldaCG, c(1, 3), c(3, 1))
#' @importFrom plyr mapvalues
#' @export
recodeClusterZ <- function(sce, from, to, altExpName = "featureSubset") {
    if (length(setdiff(from, to)) != 0) {
        stop("All values in 'from' must have a mapping in 'to'")
    }
    if (is.null(celdaClusters(sce, altExpName = altExpName))) {
        stop("Provided 'sce' argument does not have a 'celda_cell_cluster'",
            " column in 'colData(altExp(sce, altExpName))'")
    }
    celdaClusters(sce, altExpName = altExpName) <- plyr::mapvalues(
        celdaClusters(sce, altExpName = altExpName), from, to)
    return(sce)
}


# for deprecated celda model objects
.recodeClusterZ <- function(celdaMod, from, to) {
    if (length(setdiff(from, to)) != 0) {
        stop("All values in 'from' must have a mapping in 'to'")
    }
    if (is.null(celdaClusters(celdaMod)$z)) {
        stop("Provided celdaMod argument does not have a z attribute")
    }
    celdaMod@clusters$z <- plyr::mapvalues(celdaClusters(celdaMod)$z,
        from, to)
    return(celdaMod)
}


#' @title Recode feature module labels
#' @description Recode feature module clusters using a mapping in the
#'  \code{from} and \code{to} arguments.
#' @param sce \linkS4class{SingleCellExperiment} object returned from
#'  \link{celda_G} or \link{celda_CG}. Must contain column
#'  \code{celda_feature_module} in
#'  \code{\link{rowData}(altExp(sce, altExpName))}.
#' @param from Numeric vector. Unique values in the range of
#'  \code{seq(celdaModules(sce))} that correspond to the original module labels
#'  in \code{sce}.
#' @param to Numeric vector. Unique values in the range of
#'  \code{seq(celdaModules(sce))} that correspond to the new module labels.
#' @param altExpName The name for the \link{altExp} slot
#'  to use. Default "featureSubset".
#' @return @return \linkS4class{SingleCellExperiment} object with recoded
#'  feature module labels.
#' @examples
#' data(sceCeldaCG)
#' sceReorderedY <- recodeClusterY(sceCeldaCG, c(1, 3), c(3, 1))
#' @export
recodeClusterY <- function(sce, from, to, altExpName = "featureSubset") {
    if (length(setdiff(from, to)) != 0) {
        stop("All values in 'from' must have a mapping in 'to'")
    }
    if (is.null(celdaModules(sce, altExpName = altExpName))) {
        stop("Provided 'sce' argument does not have a 'celda_feature_module'",
            " column in 'rowData(altExp(sce, altExpName))'")
    }
    celdaModules(sce, altExpName = altExpName) <- plyr::mapvalues(
        celdaModules(sce, altExpName = altExpName), from, to)
    return(sce)
}


# for deprecated celda model objects
.recodeClusterY <- function(celdaMod, from, to) {
    if (length(setdiff(from, to)) != 0) {
        stop("All values in 'from' must have a mapping in 'to'")
    }
    if (is.null(celdaClusters(celdaMod)$y)) {
        stop("Provided celdaMod argument does not have a y attribute")
    }
    celdaMod@clusters$y <- plyr::mapvalues(celdaClusters(celdaMod)$y,
        from, to)
    return(celdaMod)
}


#' @title Check count matrix consistency
#' @description Checks if the counts matrix is the same one used to generate
#'  the celda model object by comparing dimensions and MD5 checksum.
#' @param counts Integer matrix. Rows represent features and columns represent
#'  cells.
#' @param celdaMod A \code{celdaModel} or \code{celdaList} object.
#' @param errorOnMismatch Logical. Whether to throw an error in the event of
#'  a mismatch. Default TRUE.
#' @param ... Ignored. Placeholder to prevent check warning.
#' @return Returns TRUE if provided count matrix matches the one used in the
#'  celda object and/or \code{errorOnMismatch = FALSE}, FALSE otherwise.
#' @export
setGeneric("compareCountMatrix", function(counts, celdaMod, ...) {
    standardGeneric("compareCountMatrix")})


#' @rdname compareCountMatrix
#' @examples
#' data(celdaCGSim, celdaCGMod)
#' compareCountMatrix(celdaCGSim$counts, celdaCGMod, errorOnMismatch = FALSE)
#' @export
setMethod("compareCountMatrix",
    signature(celdaMod = "celdaModel"),
    function(counts, celdaMod, errorOnMismatch = TRUE) {

        if ("y" %in% names(celdaClusters(celdaMod))) {
            if (nrow(counts) != length(celdaClusters(celdaMod)$y)) {
                stop(
                    "The provided celda object was generated from a counts",
                    " matrix with a different number of features than the one",
                    " provided."
                )
            }
        }

        if ("z" %in% names(celdaClusters(celdaMod))) {
            if (ncol(counts) != length(celdaClusters(celdaMod)$z)) {
                stop(
                    "The provided celda object was generated from a counts",
                    " matrix with a different number of cells than the one",
                    " provided."
                )
            }
        }
        celdaChecksum <- params(celdaMod)$countChecksum

        counts <- .processCounts(counts)
        # Checksums are generated in celdaGridSearch and model after processing
        count.md5 <- .createCountChecksum(counts)
        res <- isTRUE(count.md5 == celdaChecksum)
        if (res) {
            return(TRUE)
        }
        if (!res && errorOnMismatch) {
            stop(
                "There was a mismatch between the provided count matrix and",
                " the count matrix used to generate the provided celda result."
            )
        } else if (!res && !errorOnMismatch) {
            warning("There was a mismatch between the provided count matrix",
                " and the count matrix used to generate the provided celda",
                " result.")
            return(FALSE)
        }
    }
)


#' @rdname compareCountMatrix
#' @examples
#' data(celdaCGSim, celdaCGGridSearchRes)
#' compareCountMatrix(celdaCGSim$counts, celdaCGGridSearchRes,
#'     errorOnMismatch = FALSE)
#' @export
setMethod("compareCountMatrix",
    signature(celdaMod = "celdaList"),
    function(counts, celdaMod, errorOnMismatch = TRUE) {

        if ("y" %in% names(celdaMod@resList[[1]]@clusters)) {
            if (nrow(counts) != length(celdaMod@resList[[1]]@clusters$y)) {
                stop(
                    "The provided celda object was generated from a counts",
                    " matrix with a different number of features than the one",
                    " provided."
                )
            }
        }

        if ("z" %in% names(celdaMod@resList[[1]]@clusters)) {
            if (ncol(counts) != length(celdaMod@resList[[1]]@clusters$z)) {
                stop(
                    "The provided celda object was generated from a counts",
                    " matrix with a different number of cells than the one",
                    " provided."
                )
            }
        }
        celdaChecksum <- celdaMod@countChecksum

        counts <- .processCounts(counts)
        # Checksums are generated in celdaGridSearch and model after processing
        count.md5 <- .createCountChecksum(counts)
        res <- isTRUE(count.md5 == celdaChecksum)
        if (res) {
            return(TRUE)
        }
        if (!res && errorOnMismatch) {
            stop(
                "There was a mismatch between the provided count matrix and",
                " the count matrix used to generate the provided celda result."
            )
        } else if (!res && !errorOnMismatch) {
            warning("There was a mismatch between the provided count matrix",
                " and the count matrix used to generate the provided celda",
                " result.")
            return(FALSE)
        }
    }
)


.logMessages <- function(...,
                         sep = " ",
                         logfile = NULL,
                         append = FALSE,
                         verbose = TRUE) {
  if (isTRUE(verbose)) {
    if (!is.null(logfile)) {
      if (!is.character(logfile) || length(logfile) > 1) {
        stop(
          "The log file parameter needs to be a single character",
          " string."
        )
      }
      cat(paste(..., "\n", sep = sep),
        file = logfile,
        append = append
      )
    } else {
      message(paste(..., sep = sep))
    }
  }
}


#' @title Create a color palette
#' @description Generate a palette of `n` distinct colors.
#' @param n Integer. Number of colors to generate.
#' @param hues Character vector. Colors available from `colors()`. These will
#'  be used as the base colors for the clustering scheme in HSV. Different
#'  saturations and values will be generated for each hue. Default c("red",
#'  "cyan", "orange", "blue", "yellow", "purple", "green", "magenta").
#' @param saturationRange Numeric vector. A vector of length 2 denoting the
#'  saturation for HSV. Values must be in [0,1]. Default: c(0.25, 1).
#' @param valueRange Numeric vector. A vector of length 2 denoting the range
#'  of values for HSV. Values must be in [0,1]. Default: `c(0.5, 1)`.
#' @return A vector of distinct colors that have been converted to HEX from HSV.
#' @examples
#' colorPal <- distinctColors(6) # can be used in plotting functions
#' @importFrom grDevices colors
#' @importFrom grDevices rgb2hsv
#' @importFrom grDevices hsv
#' @export
distinctColors <- function(n,
                           hues = c(
                             "red",
                             "cyan",
                             "orange",
                             "blue",
                             "yellow",
                             "purple",
                             "green",
                             "magenta"
                           ),
                           saturationRange = c(0.7, 1),
                           valueRange = c(0.7, 1)) {
  if (!(all(hues %in% grDevices::colors()))) {
    stop(
      "Only color names listed in the 'color' function can be used in",
      " 'hues'"
    )
  }
  # Convert R colors to RGB and then to HSV color format
  huesHsv <- grDevices::rgb2hsv(grDevices::col2rgb(hues))
  # Calculate all combination of saturation/value pairs
  # Note that low saturation with low value (i.e. high darkness) is too dark
  # for all hues. Likewise, high saturation with high value (i.e. low
  # darkness) is hard to distinguish. Therefore, saturation and value are
  # set to be anticorrelated
  numVs <- ceiling(n / length(hues))
  s <- seq(
    from = saturationRange[1],
    to = saturationRange[2],
    length = numVs
  )
  v <- seq(
    from = valueRange[2],
    to = valueRange[1],
    length = numVs
  )
  # Create all combination of hues with saturation/value pairs
  list <- lapply(seq(numVs), function(x) {
    rbind(huesHsv[1, ], s[x], v[x])
  })
  newHsv <- do.call(cbind, list)
  # Convert to hex
  col <- grDevices::hsv(newHsv[1, ], newHsv[2, ], newHsv[3, ])
  return(col[seq(n)])
}


.processCounts <- function(counts) {
  counts <- as.matrix(counts)
  if (typeof(counts) != "integer") {
    counts <- round(counts)
    storage.mode(counts) <- "integer"
  }
  return(counts)
}


# Perform some simple checks on the counts matrix, to ensure celda modeling
# expectations are met
.validateCounts <- function(counts) {
  # And each row/column of the count matrix must have at least one count
  countRowSum <- rowSums(counts)
  countColSum <- colSums(counts)
  if (sum(countRowSum == 0) > 0 | sum(countColSum == 0) > 0) {
    stop(
      "Each row and column of the count matrix must have at least",
      " one count"
    )
  }
}


# Wrapper function, creates checksum for matrix.
# Feature names, cell names are not taken into account.
#' @importFrom digest digest
.createCountChecksum <- function(counts) {
  rownames(counts) <- NULL
  colnames(counts) <- NULL
  countChecksum <- digest::digest(counts, algo = "md5")
  return(countChecksum)
}


# Generate n random deviates from the Dirichlet function with shape parameters
# alpha. Adapted from gtools v3.5
.rdirichlet <- function(n, alpha) {
  l <- length(alpha)
  x <- matrix(stats::rgamma(l * n, alpha),
    ncol = l,
    byrow = TRUE
  )
  # Check for case where all sampled entries are zero due to round off
  # One entry will be randomly chosen to be one
  isZero <- rowSums(x) == 0
  assignment <- sample(seq(l), size = sum(isZero), replace = TRUE)
  x[cbind(which(isZero), assignment)] <- 1
  # Normalize
  sm <- x %*% rep(1, l)
  y <- x / as.vector(sm)
  return(y)
}


# Make sure provided sample labels are the right type,
# or generate some if none were provided
.processSampleLabels <- function(sampleLabel, numCells) {
  if (is.null(sampleLabel)) {
    sampleLabel <- as.factor(rep("Sample_1", numCells))
  } else {
    if (length(sampleLabel) != numCells) {
      stop(
        "'sampleLabel' must be the same length as the number of",
        " columns in the 'counts' matrix."
      )
    }
  }

  if (!is.factor(sampleLabel)) {
    sampleLabel <- as.factor(sampleLabel)
  }
  return(sampleLabel)
}


#' @title Output a feature module table
#' @description Creates a table that contains the list of features in
#'  each feature module.
#' @param sce A \linkS4class{SingleCellExperiment} object returned by
#'  \link{celda_G}, or \link{celda_CG}, with the matrix
#'  located in the \code{useAssay} assay slot.
#'  Rows represent features and columns represent cells.
#' @param useAssay A string specifying which \link{assay}
#'  slot to use. Default "counts".
#' @param altExpName The name for the \link{altExp} slot
#'  to use. Default "featureSubset".
#' @param outputFile File name for feature module table. If NULL, file will
#'  not be created. Default NULL.
#' @return Matrix. Contains a list of features per each column (feature module)
#' @examples
#' data(sceCeldaCG)
#' featureModuleTable(sceCeldaCG)
#' @importFrom stringi stri_list2matrix
#' @export
featureModuleTable <- function(sce, useAssay = "counts",
    altExpName = "featureSubset", outputFile = NULL) {

  factorizeMatrix <- factorizeMatrix(sce,
      useAssay = useAssay,
      altExpName = altExpName,
      type = "proportion")
  allGenes <- topRank(factorizeMatrix$proportions$module, n = nrow(sce))
  res <- as.data.frame(stringi::stri_list2matrix(allGenes$names))
  res <- apply(res, c(1, 2), function(x) {
    if (is.na(x)) {
      return("")
    } else {
      return(x)
    }
  })
  colnames(res) <- gsub(
    pattern = "V",
    replacement = "L",
    x = colnames(res)
  )
  if (is.null(outputFile)) {
    return(res)
  } else {
    utils::write.table(
      res,
      file = outputFile,
      sep = "\t",
      row.names = FALSE,
      quote = FALSE
    )
  }
}


#' @title Retrieve row index for a set of features
#' @description This will return indices of features among the rownames
#' or rowData of a data.frame, matrix, or a \linkS4class{SummarizedExperiment}
#' object including a \linkS4class{SingleCellExperiment}.
#' Partial matching (i.e. grepping) can be used by setting
#' \code{exactMatch = FALSE}.
#' @param features Character vector of feature names to find in the rows of
#' \code{x}.
#' @param x A data.frame, matrix, or \linkS4class{SingleCellExperiment}
#' object to search.
#' @param by Character. Where to search for features in \code{x}. If set to
#' \code{"rownames"} then the features will be searched for among
#' \code{rownames(x)}. If \code{x} inherits from class
#' \linkS4class{SummarizedExperiment}, then \code{by} can be one of the
#' fields in the row annotation data.frame (i.e. one of
#' \code{colnames(rowData(x))}).
#' @param exactMatch Boolean. Whether to only identify exact matches
#' or to identify partial matches using \code{\link{grep}}.
#' @param removeNA Boolean. If set to \code{FALSE}, features not found in
#' \code{x} will be given \code{NA} and the returned vector will be the same
#' length as \code{features}. If set to \code{TRUE}, then the \code{NA}
#' values will be removed from the returned vector. Default \code{FALSE}.
#' @return A vector of row indices for the matching features in \code{x}.
#' @author Yusuke Koga, Joshua Campbell
#' @seealso '\link[scater]{retrieveFeatureInfo}' from package \code{'scater'}
#' and \code{link{regex}} for how to use regular expressions when
#' \code{exactMatch = FALSE}.
#' @examples
#' data(celdaCGSim)
#' retrieveFeatureIndex(c("Gene_1", "Gene_5"), celdaCGSim$counts)
#' retrieveFeatureIndex(c("Gene_1", "Gene_5"), celdaCGSim$counts,
#'                                             exactMatch = FALSE)
#' @export
retrieveFeatureIndex <- function(features,
                                 x,
                                 by = "rownames",
                                 exactMatch = TRUE,
                                 removeNA = FALSE) {

  # Extract vector to search through
  if (by == "rownames") {
    if (is.null(rownames(x))) {
      stop("'rownames' of 'x' are 'NULL'. Please set 'rownames' or change",
           " 'by' to search a different column in 'x'.")
    }
    search <- rownames(x)
  } else if (length(ncol(x)) > 0) {
    if (inherits(x, "SummarizedExperiment")) {
      if (!(by %in% colnames(rowData(x)))) {
        stop("'by' is not a column in 'rowData(x)'.")
      }
      search <- rowData(x)[, by]
    } else {
      if (!(by %in% colnames(x))) {
        stop("'by' is not a column in 'x'.")
      }
      search <- x[, by]
    }
  } else {
    search <- as.character(x)
  }

  # Match each element of 'pattern' in vector 'search'
  if (!isTRUE(exactMatch)) {
    featuresIndices <- rep(NA, length(features))
    for (i in seq_along(features)) {
      g <- grep(features[i], search)
      if (length(g) == 1) {
        featuresIndices[i] <- g
      } else if (length(g) > 1) {
        warning(
          "Feature '", features[i], "' matched multiple items in '",
          by, "': ", paste(search[g], collapse = ","),
          ". Only the first match will be selected."
        )
        featuresIndices[i] <- g[1]
      }
    }
  } else {
    featuresIndices <- match(features, search)
  }

  if (sum(is.na(featuresIndices)) > 0) {
    if (sum(is.na(featuresIndices)) == length(features)) {
      if (isTRUE(exactMatch)) {
        stop(
          "None of the provided features had matching",
          " items in '", by, "' within 'x'. ",
          "Check the spelling or try setting",
          " 'exactMatch = FALSE'."
        )
      } else {
        stop(
          "None of the provided features had matching",
          " items in '", by, "' within 'x'. ",
          "Check the spelling and make sure 'by' is set",
          " to the appropriate place in 'x'."
        )
      }
    }
    warning(
      "The following features were not present in 'x': ",
      paste(features[which(is.na(featuresIndices))],
        collapse = ","
      )
    )
  }

  if (isTRUE(removeNA)) {
    featuresIndices <- featuresIndices[!is.na(featuresIndices)]
  }
  return(featuresIndices)
}