Add license
[loomr.git] / R / internal.R
index 4f590f7f3daf245f44051b1dcef20250f25668df..b572fbb46e28362f3e39268fdcbc56c2fd68c18f 100644 (file)
@@ -1,3 +1,20 @@
+# Copyright 2017, 2018 Paul Hoffman <phoffman@nygenome.org>
+#
+# This file is part of loomR.
+#
+# loomR is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# loomR is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with loomR.  If not, see <https://www.gnu.org/licenses/>.
+
 # Generate chunk points
 #
 # @param data.size How big is the data being chunked
@@ -18,12 +35,36 @@ chunkPoints <- function(data.size, chunk.size) {
   ))
 }
 
+
+# Convert sparse matrix pointers to column indicies
+#
+# A function to get the column (j) indices of a sparse matrix from the pointers (p)
+#
+# @param p A vector of sparse matrix pointers
+#
+# @return A vector of column (j) indices
+#
+# @examples
+# dat <- c(0, 0, 1, 4, 0, 2, 0, 9, 0)
+# smat <- Matrix::Matrix(data = dat, nrow = 3, sparse = TRUE)
+# PointerToIndex(p = smat@p)
+#
+PointerToIndex <- function(p) {
+  # From StackOverflow
+  # https://stackoverflow.com/questions/20008200/r-constructing-sparse-matrix
+  dp <- diff(x = p)
+  j <- rep.int(x = seq_along(along.with = dp), times = dp)
+  return(j)
+}
+
 # Get HDF5 data types
 #
-# @param An R object
+# @param x An R object or string describing HDF5 datatype
 #
 # @return The corresponding HDF5 data type
 #
+# @ rdname getDtype
+#
 #' @import hdf5r
 #
 # @seealso \code\link{hdf5r::h5types}
@@ -34,11 +75,25 @@ getDtype <- function(x) {
     'numeric' = h5types$double,
     'integer' = h5types$int,
     'character' = H5T_STRING$new(size = Inf),
-    'logical' = h5types$hbool_t,
+    'logical' = H5T_LOGICAL$new(),
     stop(paste("Unknown data type:", class(x = x)))
   ))
 }
 
+# @describeIn getDtype A version of getDtype that works specifically for hdf5r types,
+# useful for getDtype2(x = class(x = object[['dataset']]$get_type())[1])
+#
+getDtype2 <- function(x) {
+  return(getDtype(x = switch(
+    EXPR = x,
+    'H5T_FLOAT' = numeric(),
+    'H5T_INTEGER' = integer(),
+    'H5T_STRING' = character(),
+    'H5T_LOGICAL' = logical(),
+    stop(paste("Unkown data type:", x))
+  )))
+}
+
 # Validate a loom object
 #
 # @param object A loom object
@@ -61,7 +116,7 @@ validateLoom <- function(object) {
     stop("The root dataset must be called '/matrix'")
   }
   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
-  required.groups <- c('row_attrs', 'col_attrs', 'layers', 'row_edges', 'col_edges')
+  required.groups <- c('row_attrs', 'col_attrs', 'layers', 'row_graphs', 'col_graphs')
   dim.matrix <- object[['matrix']]$dims # Columns x Rows
   names(x = dim.matrix) <- required.groups[c(2, 1)]
   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
@@ -75,6 +130,20 @@ validateLoom <- function(object) {
     "Reopen in 'r+' mode to automatically add missing groups",
     sep = '\n'
   )
+  if (any(grepl(pattern = 'edges', x = root.groups))) {
+    if (object$mode != 'r') {
+      cate("Moving edge groups to graph groups to conform to loom v2.0.1")
+      edges <- grep(pattern = 'edges', x = root.groups, value = TRUE)
+      for (group in edges) {
+        graph <- gsub(pattern = 'edges', replacement = 'graphs', x = group)
+        object$link_move_to(dst_loc = object, dst_name = graph, src_name = group)
+      }
+      root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
+    } else {
+      object$close_all()
+      stop(reopen.msg)
+    }
+  }
   if (length(x = root.groups) > length(x = required.groups)) {
     stop(group.msg)
   } else if (length(x = root.groups) < length(x = required.groups)) {
@@ -86,6 +155,7 @@ validateLoom <- function(object) {
         }
         root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
       } else {
+        object$close_all()
         stop(reopen.msg)
       }
     } else {
@@ -113,6 +183,34 @@ validateLoom <- function(object) {
       }
     }
   }
+  # Check row and column graphs
+  graph.groups <- grep(pattern = 'graphs', x = required.groups, value = TRUE)
+  graph.msg <- "There can only be three datasets in a graph: 'a', 'b', and 'w'"
+  for (group in graph.groups) {
+    group.datasets <- list.datasets(object = object[[group]], recursive = FALSE)
+    if (length(x = group.datasets) > 0) {
+      stop(paste("All datasets in", group, "must be in a graph group"))
+    }
+    graphs <- list.groups(object = object[[group]], full.names = TRUE, recursive = FALSE)
+    for (graph in graphs) {
+      graph.datasets <- list.datasets(object = object[[graph]], full.names = TRUE)
+      if (length(x = graph.datasets) != 3) {
+        stop(graph.msg)
+      }
+      if (!all(basename(path = graph.datasets) %in% c('a', 'b', 'w'))) {
+        stop(graph.msg)
+      }
+      graph.lengths <- lapply(
+        X = graph.datasets,
+        FUN = function(dset) {
+          return(object[[dset]]$dims)
+        }
+      )
+      if (length(x = unique(x = graph.lengths)) != 1) {
+        stop("All graph datasets must be the same length")
+      }
+    }
+  }
   # Check layers
   for (dataset in list.datasets(object = object[['/layers']])) {
     if (any(object[[paste('layers', dataset, sep = '/')]]$dims != dim.matrix)) {
@@ -375,12 +473,114 @@ new.pb <- function() {
   return(txtProgressBar(style = 3, char = '='))
 }
 
+# Break a vector into a list of consecutive values
+#
+# @param x An integer vector
+# @param max.length Maximum length of consecutive values;
+# set to NULL if no maximum is desired
+#
+# @return A list where each value is a run of consecutive values from x
+#
+#' @importFrom stats na.omit
+#
+break.consecutive <- function(x, max.length = NULL, min.length = NULL) {
+  # Coerce x to integers, find unique values, sort it
+  x <- sort(x = unique(x = na.omit(object = as.integer(x = x))))
+  if (is.null(x = max.length)) {
+    max.length <- length(x = x)
+  }
+  if (is.null(x = min.length)) {
+    min.length <- 1L
+  }
+  # Functions to allocate a vector and increment a counter
+  alloc <- function() {return(vector(mode = 'integer', length = max.length))}
+  inc <- function(x) {return(x + 1L)}
+  # Results list, counters, and holding vector, all preallocated for speed
+  consecutive <- vector(mode = 'list', length = length(x = x))
+  index <- 1L
+  counter <- 0L
+  run <- alloc()
+  # For every value in x, if it's consecutive, add it to the holding vector
+  # Otherwise, add the run to the list of consecutive values
+  for (i in x) {
+    if (counter == max.length || (counter != 0L && counter >= min.length && i != run[counter] + 1L)) {
+      consecutive[[index]] <- run[1:counter]
+      index <- inc(x = index)
+      run <- alloc()
+      counter <- 0L
+    }
+    counter <- inc(x = counter)
+    run[counter] <- i
+  }
+  consecutive[[index]] <- run[1:counter]
+  consecutive <- consecutive[1:index]
+  return(consecutive)
+}
+
+# Calculate nUMI and nGene
+#
+# @param object An object
+# @param chunk.size Number of cells to chunk over
+# @param is.expr Expression threshold for 'detected' gene. For most datasets, particularly UMI datasets, will be set to 0 (default).
+# If not, when initializing, this should be set to a level based on pre-normalized counts (i.e. require at least 5 counts to be treated as expresesd).
+# @param display.progres Display a progress bar?
+#
+# @return Stores nUMI in \code{col_attrs/nUMI}; stores nGene in \code{col_attrs/nGene}; will overwrite any existing datasets at these locations
+#
+#' @importFrom Matrix rowSums
+#
+calc.umi <- function(
+  object,
+  chunk.size = 1000,
+  is.expr = 0,
+  display.progress = TRUE
+) {
+  if (display.progress) {
+    cate("Calculating number of UMIs per cell")
+  }
+  object$apply(
+    name = 'col_attrs/nUMI',
+    FUN = rowSums,
+    MARGIN = 2,
+    chunk.size = chunk.size,
+    dataset.use = 'matrix',
+    overwrite = TRUE,
+    display.progress = display.progress
+  )
+  if (display.progress) {
+    cate("Calculating number of genes expressed per cell")
+    cate("Using a threshold of", is.expr, "for gene expression")
+  }
+  object$apply(
+    name = 'col_attrs/nGene',
+    FUN = function(mat, is.expr) {
+      return(rowSums(x = mat > is.expr))
+    },
+    MARGIN = 2,
+    chunk.size = chunk.size,
+    dataset.use = 'matrix',
+    overwrite = TRUE,
+    display.progress = display.progress,
+    is.expr = is.expr
+  )
+  invisible(x = object)
+}
+
 # Cat with a new line
 #
 # @param ... Text to be output
 #
 catn <- function(...) {
-  cat(..., '\n')
+  x = list(...)
+  if (length(x = x)) {
+    if (!is.null(x = names(x = x)) && length(x = x) == 1 && names(x = x) == 'file') {
+      cat(...)
+    } else {
+      cat(..., '\n')
+    }
+  } else {
+    cat()
+  }
 }
 
 # Cat to stderr