index 76910a45d87a2e2af228f87b6a030250591a4bca..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
#
@@ -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()))
+#
+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,17 +116,51 @@ 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')
+  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)
group.msg <- paste0(
-    "There can only be three groups in the loom file: '",
+    paste("There can only be", length(x = required.groups), "groups in the loom file: '"),
paste(required.groups, collapse = "', '"),
"'"
)
-  if (length(x = root.groups) != 3) {
+  reopen.msg <- paste(
+    group.msg,
+    "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)) {
+    if (all(root.groups %in% required.groups)) {
+      if (object\$mode != 'r') {
+        missing.groups <- required.groups[!(required.groups %in% root.groups)]
+        for (group in missing.groups) {
+          object\$create_group(name = group)
+        }
+        root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
+      } else {
+        object\$close_all()
+        stop(reopen.msg)
+      }
+    } else {
+      stop(group.msg)
+    }
}
if (!all(required.groups %in% root.groups)) {
stop(group.msg)
@@ -94,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)) {
@@ -356,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