1 # Copyright 2017, 2018 Paul Hoffman <phoffman@nygenome.org>
2 #
3 # This file is part of loomR.
4 #
5 # loomR is free software: you can redistribute it and/or modify
6 # it under the terms of the GNU General Public License as published by
7 # the Free Software Foundation, either version 3 of the License, or
8 # (at your option) any later version.
9 #
10 # loomR is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 # GNU General Public License for more details.
14 #
15 # You should have received a copy of the GNU General Public License
16 # along with loomR.  If not, see <https://www.gnu.org/licenses/>.
18 # Generate chunk points
19 #
20 # @param data.size How big is the data being chunked
21 # @param chunk.size How big should each chunk be
22 #
23 # @return A matrix where each column is a chunk, row 1 is start points, row 2 is end points
24 #
25 chunkPoints <- function(data.size, chunk.size) {
26   return(vapply(
27     X = 1L:ceiling(data.size / chunk.size),
28     FUN = function(i) {
29       return(c(
30         start = (chunk.size * (i - 1L)) + 1L,
31         end = min(chunk.size * i, data.size)
32       ))
33     },
34     FUN.VALUE = numeric(length = 2L)
35   ))
36 }
39 # Convert sparse matrix pointers to column indicies
40 #
41 # A function to get the column (j) indices of a sparse matrix from the pointers (p)
42 #
43 # @param p A vector of sparse matrix pointers
44 #
45 # @return A vector of column (j) indices
46 #
47 # @examples
48 # dat <- c(0, 0, 1, 4, 0, 2, 0, 9, 0)
49 # smat <- Matrix::Matrix(data = dat, nrow = 3, sparse = TRUE)
50 # PointerToIndex(p = smat@p)
51 #
52 PointerToIndex <- function(p) {
53   # From StackOverflow
54   # https://stackoverflow.com/questions/20008200/r-constructing-sparse-matrix
55   dp <- diff(x = p)
56   j <- rep.int(x = seq_along(along.with = dp), times = dp)
57   return(j)
58 }
60 # Get HDF5 data types
61 #
62 # @param x An R object or string describing HDF5 datatype
63 #
64 # @return The corresponding HDF5 data type
65 #
66 # @ rdname getDtype
67 #
68 #' @import hdf5r
69 #
71 #
72 getDtype <- function(x) {
73   return(switch(
74     EXPR = class(x = x),
75     'numeric' = h5types\$double,
76     'integer' = h5types\$int,
77     'character' = H5T_STRING\$new(size = Inf),
78     'logical' = H5T_LOGICAL\$new(),
79     stop(paste("Unknown data type:", class(x = x)))
80   ))
81 }
83 # @describeIn getDtype A version of getDtype that works specifically for hdf5r types,
84 # useful for getDtype2(x = class(x = object[['dataset']]\$get_type()))
85 #
86 getDtype2 <- function(x) {
87   return(getDtype(x = switch(
88     EXPR = x,
89     'H5T_FLOAT' = numeric(),
90     'H5T_INTEGER' = integer(),
91     'H5T_STRING' = character(),
92     'H5T_LOGICAL' = logical(),
93     stop(paste("Unkown data type:", x))
94   )))
95 }
97 # Validate a loom object
98 #
99 # @param object A loom object
100 #
101 # @return None, errors out if object is an invalid loom connection
102 #
104 #
105 validateLoom <- function(object) {
106   if (!inherits(x = object, what = 'loom')) {
107     stop("No need to validate a non-loom object")
108   }
109   # A loom file is a specific HDF5
110   # We need a dataset in /matrix that's a two-dimensional dense matrix
111   root.datasets <- list.datasets(object = object, path = '/', recursive = FALSE)
112   if (length(x = root.datasets) != 1) {
113     stop("There can only be one dataset at the root of the loom file")
114   }
115   if (root.datasets != 'matrix') {
116     stop("The root dataset must be called '/matrix'")
117   }
118   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
119   required.groups <- c('row_attrs', 'col_attrs', 'layers', 'row_graphs', 'col_graphs')
120   dim.matrix <- object[['matrix']]\$dims # Columns x Rows
121   names(x = dim.matrix) <- required.groups[c(2, 1)]
122   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
123   group.msg <- paste0(
124     paste("There can only be", length(x = required.groups), "groups in the loom file: '"),
125     paste(required.groups, collapse = "', '"),
126     "'"
127   )
128   reopen.msg <- paste(
129     group.msg,
130     "Reopen in 'r+' mode to automatically add missing groups",
131     sep = '\n'
132   )
133   if (any(grepl(pattern = 'edges', x = root.groups))) {
134     if (object\$mode != 'r') {
135       cate("Moving edge groups to graph groups to conform to loom v2.0.1")
136       edges <- grep(pattern = 'edges', x = root.groups, value = TRUE)
137       for (group in edges) {
138         graph <- gsub(pattern = 'edges', replacement = 'graphs', x = group)
139         object\$link_move_to(dst_loc = object, dst_name = graph, src_name = group)
140       }
141       root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
142     } else {
143       object\$close_all()
144       stop(reopen.msg)
145     }
146   }
147   if (length(x = root.groups) > length(x = required.groups)) {
148     stop(group.msg)
149   } else if (length(x = root.groups) < length(x = required.groups)) {
150     if (all(root.groups %in% required.groups)) {
151       if (object\$mode != 'r') {
152         missing.groups <- required.groups[!(required.groups %in% root.groups)]
153         for (group in missing.groups) {
154           object\$create_group(name = group)
155         }
156         root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
157       } else {
158         object\$close_all()
159         stop(reopen.msg)
160       }
161     } else {
162       stop(group.msg)
163     }
164   }
165   if (!all(required.groups %in% root.groups)) {
166     stop(group.msg)
167   }
168   # Check row and column attributes
169   for (group in required.groups[1:2]) {
170     # No subgroups
171     if (length(x = list.groups(object = object[[group]], recursive = FALSE)) > 0) {
172       stop(paste("Group", group, "cannot have subgroups"))
173     }
174     # All datasets must have their first (last) dimmension equal to M(row) or N(column)
175     for (dataset in list.datasets(object = object[[group]])) {
176       dataset.dim <- object[[group]][[dataset]]\$dims
177       dataset.dim <- dataset.dim[length(x = dataset.dim)]
178       if (dataset.dim != dim.matrix[group]) {
179         print(dataset)
180         print(object[[group]][[dataset]])
181         print(dim.matrix)
182         stop("All datasets in group ", group, " must be of length ", dim.matrix[group])
183       }
184     }
185   }
186   # Check row and column graphs
187   graph.groups <- grep(pattern = 'graphs', x = required.groups, value = TRUE)
188   graph.msg <- "There can only be three datasets in a graph: 'a', 'b', and 'w'"
189   for (group in graph.groups) {
190     group.datasets <- list.datasets(object = object[[group]], recursive = FALSE)
191     if (length(x = group.datasets) > 0) {
192       stop(paste("All datasets in", group, "must be in a graph group"))
193     }
194     graphs <- list.groups(object = object[[group]], full.names = TRUE, recursive = FALSE)
195     for (graph in graphs) {
196       graph.datasets <- list.datasets(object = object[[graph]], full.names = TRUE)
197       if (length(x = graph.datasets) != 3) {
198         stop(graph.msg)
199       }
200       if (!all(basename(path = graph.datasets) %in% c('a', 'b', 'w'))) {
201         stop(graph.msg)
202       }
203       graph.lengths <- lapply(
204         X = graph.datasets,
205         FUN = function(dset) {
206           return(object[[dset]]\$dims)
207         }
208       )
209       if (length(x = unique(x = graph.lengths)) != 1) {
210         stop("All graph datasets must be the same length")
211       }
212     }
213   }
214   # Check layers
215   for (dataset in list.datasets(object = object[['/layers']])) {
216     if (any(object[[paste('layers', dataset, sep = '/')]]\$dims != dim.matrix)) {
217       stop(paste("All datasets in '/layers' must be", dim.matrix, 'by', dim.matrix))
218     }
219   }
220 }
222 # A function to determine if a vector is a vector and not a list
223 #
224 # @param x An object
225 #
226 # @return TRUE if 'x' is a vector or a factor, otherwise FALSE
227 #
228 is.actual_vector <- function(x) {
229   return((is.vector(x = x) || is.factor(x = x)) && !is.list(x = x))
230 }
232 # Check additions to /matrix
233 #
234 # @param x A list of vectors to add to /matrix
235 # @param n The number of genes needed in each cell
236 #
237 # @return 'x' as a list of vectors
238 #
239 check.matrix_data <- function(x, n) {
240   # Coerce x into a list, where each
241   # entry in the list is a new cell added
242   if (is.actual_vector(x = x)) {
243     x <- list(x)
244   } else if (is.matrix(x = x) || is.data.frame(x = x)) {
245     x <- as.list(x = x)
246   }
247   if (!is.list(x = x)) {
248     stop("Matrix data must be a list of vectors")
249   }
250   # Ensure that each entry in the list is a vector of length n
251   vector.check <- vapply(
252     X = x,
253     FUN = is.actual_vector,
254     FUN.VALUE = logical(length = 1L)
255   )
256   if (!all(vector.check)) {
257     stop('Each new cell added must be represented as a vector')
258   }
259   # Ensure each new cell added has data for the number of genes present
260   for (i in 1:length(x = x)) {
262     if (length(x = cell.add) > n) {
263       stop(paste(
264         "Cannot add genes to a loom file, the maximum number of genes allowed is",
265         n
266       ))
267     } else if (length(x = cell.add) < n) {
269     }
271   }
272   return(x)
273 }
275 # Get the number of cells being added to /matrix
276 #
277 # @param x A list of vectors to add to /matrix
278 #
279 # @return The number of cells in x
280 #
281 nCells.matrix_data <- function(x) {
282   return(length(x = x))
283 }
286 #
287 # @param x A list of vectors to add to /matrix
288 # @param n The number of genes each cell needs
289 # @param m2 The number of cells being added to the loom file
290 #
291 # @return 'x' with the proper number of cells
292 #
293 addCells.matrix_data <- function(x, n, m2) {
294   if (length(x = x) < m2) {
295     x[(length(x = x) + 1):m2] <- list(rep.int(x = NA, times = n))
296   }
297   return(x)
298 }
300 # Check additions to /layers
301 #
302 # @param x A list of matrices to add layers in /layers
303 # @param n The number of genes needed for each layer
304 # @param layers.names Names found in /layers
305 #
306 # @return 'x' as a list of matricies with 'n' rows for each layer present in /layers
307 #
308 check.layers <- function(x, n, layers.names) {
309   # Coerce x into a list of matricies
310   if (is.null(x = x)) {
311     x <- list()
312   } else if (is.matrix(x = x) || is.data.frame(x = x)) {
313     x <- list(as.matrix(x = x))
314   }
315   if (!is.list(x = x)) {
316     stop("Layers data must be a list of matricies")
317   }
318   # Ensure we have enough layer additions for each layer
319   # Manage named lists, taking only those with the same name as in /layers
320   # or is named with an empty string
321   if (!is.null(x = names(x = x))) {
322     x.use <- which(x = names(x = x) %in% layers.names | names(x = x) == '')
323     x <- x[x.use]
324   }
325   if (length(x = x) > length(x = layers.names)) {
327   } else if (length(x = x) < length(x = layers.names)) {
328     x[(length(x = x) + 1):length(x = layers.names)] <- data.frame(rep.int(x = NA, times = n))
329   }
330   # Ensure that we have all genes needed
331   for (i in 1:length(x = x)) {
332     layer <- x[[i]]
333     if (is.data.frame(x = layer)) {
334       layer <- as.matrix(x = layer)
335     } else if (is.actual_vector(x = layer)) {
336       layer <- matrix(data = layer, ncol = 1)
337     }
338     if (!is.matrix(x = layer)) {
339       stop("Layers data must be a list of matrices")
340     }
341     if (nrow(x = layer) > n) {
342       stop(paste(
343         "Cannot add genes to a loom file, the maximum number of genes allowed is",
344         n
345       ))
346     } else if (nrow(x = layer) < n) {
347       layer <- as.data.frame(x = layer)
348       layer[(nrow(x = layer) + 1):n, ] <- NA
349       layer <- as.matrix(x = layer)
350     }
351     x[[i]] <- layer
352   }
353   # Set names
354   x.unnamed <- which(x = !(names(x = x) %in% layers.names))
355   if (length(x = x.unnamed) == 0) {
356     x.unnamed <- 1:length(x = x)
357   }
358   names.unused <- which(x = !(layers.names %in% names(x = x)))
359   names(x = x)[x.unnamed] <- layers.names[names.unused]
360   return(x)
361 }
363 # Get the number of cells being added to /layers
364 #
365 # @param x A list of matricies to add to /layers
366 #
367 # @return The number of cells within each matrix
368 #
369 nCells.layers <- function(x) {
370   return(vapply(X = x, FUN = ncol, FUN.VALUE = integer(length = 1L)))
371 }
374 #
375 # @param x A list of matricies to add to /layers
376 # @param n The number of genes each cell needs
377 # @param m2 The number of cells being added to the loom file
378 #
379 # @return 'x' with the proper number of cells
380 #
381 addCells.layers <- function(x, n, m2) {
382   layers.extend <- vapply(X = x, FUN = ncol, FUN.VALUE = integer(length = 1L))
383   layers.extend <- which(x = layers.extend != m2)
384   for (i in layers.extend) {
385     layer <- x[[i]]
386     layer.new <- matrix(nrow = n, ncol = m2)
387     layer.new[, 1:ncol(x = layer)] <- layer
388     x[[i]] <- layer.new
389     gc(verbose = FALSE)
390   }
391   return(x)
392 }
394 # Check additions to /col_attrs
395 #
396 # @param x A list of vectors to add to /col_attrs
397 # @param attrs.names Names of attributes found in /col_attrs
398 #
399 # @return 'x' as a list of vectors for each attribute found in /col_attrs
400 #
401 check.col_attrs <- function(x, attrs.names) {
402   # Coerce x into a list of vectors
403   if (is.null(x = x)) {
404     x <- list()
405   } else if (is.actual_vector(x = x)) {
406     x <- list(x)
407   } else if (is.matrix(x = x) || is.data.frame(x = x)) {
408     x <- as.list(x = x)
409   }
410   if (!is.list(x = x)) {
411     stop("Attribute data must be a list of vectors")
412   }
413   # Ensure we have enough attribute additions for each col_attr
414   # Manage named lists, taking only those with the same name as in /col_attrs
415   # or is named with an empty string
416   if (!is.null(x = names(x = x))) {
417     x.use <- which(x = names(x = x) %in% attrs.names | names(x = x) == '')
418     x <- x[x.use]
419   }
420   if (length(x = x) > length(x = attrs.names)) {
422   } else if (length(x = x) < length(x = attrs.names)) {
423     x[(length(x = x) + 1):length(x = attrs.names)] <- NA
424   }
425   if (!all(vapply(X = x, FUN = is.actual_vector, FUN.VALUE = logical(length = 1L)))) {
426     stop("Attribute data must be a list of vectors")
427   }
428   # Set names
429   x.unnamed <- which(x = !(names(x = x) %in% attrs.names))
430   if (length(x = x.unnamed) == 0) {
431     x.unnamed <- 1:length(x = x)
432   }
433   names.unused <- which(x = !(attrs.names %in% names(x = x)))
434   names(x = x)[x.unnamed] <- attrs.names[names.unused]
435   return(x)
436 }
438 # Get the number of cells being added to /col_attrs
439 #
440 # @param x A list of vectors to add to /col_attrs
441 #
442 # @return The number of cells for each attribute
443 #
444 nCells.col_attrs <- function(x) {
445   return(vapply(X = x, FUN = length, FUN.VALUE = integer(length = 1L)))
446 }
449 #
450 # @param x A list of vectors to add to /col_attrs
451 # @param m2 The number of cells being added to the loom file
452 #
453 # @return 'x' with the proper number of cells
454 #
455 addCells.col_attrs <- function(x, m2) {
456   attrs.extend <- vapply(X = x, FUN = length, FUN.VALUE = integer(length = 1L))
457   attrs.extend <- which(x = attrs.extend != m2)
458   for (i in attrs.extend) {
459     attr <- x[[i]]
460     attr <- c(attr, rep.int(x = NA, times = m2 - length(x = attr)))
461     x[[i]] <- attr
462   }
463   return(x)
464 }
466 # Create a progress bar
467 #
468 # @return A progress bar
469 #
470 #' @importFrom utils txtProgressBar
471 #
472 new.pb <- function() {
473   return(txtProgressBar(style = 3, char = '='))
474 }
476 # Break a vector into a list of consecutive values
477 #
478 # @param x An integer vector
479 # @param max.length Maximum length of consecutive values;
480 # set to NULL if no maximum is desired
481 #
482 # @return A list where each value is a run of consecutive values from x
483 #
484 #' @importFrom stats na.omit
485 #
486 break.consecutive <- function(x, max.length = NULL, min.length = NULL) {
487   # Coerce x to integers, find unique values, sort it
488   x <- sort(x = unique(x = na.omit(object = as.integer(x = x))))
489   if (is.null(x = max.length)) {
490     max.length <- length(x = x)
491   }
492   if (is.null(x = min.length)) {
493     min.length <- 1L
494   }
495   # Functions to allocate a vector and increment a counter
496   alloc <- function() {return(vector(mode = 'integer', length = max.length))}
497   inc <- function(x) {return(x + 1L)}
498   # Results list, counters, and holding vector, all preallocated for speed
499   consecutive <- vector(mode = 'list', length = length(x = x))
500   index <- 1L
501   counter <- 0L
502   run <- alloc()
503   # For every value in x, if it's consecutive, add it to the holding vector
504   # Otherwise, add the run to the list of consecutive values
505   for (i in x) {
506     if (counter == max.length || (counter != 0L && counter >= min.length && i != run[counter] + 1L)) {
507       consecutive[[index]] <- run[1:counter]
508       index <- inc(x = index)
509       run <- alloc()
510       counter <- 0L
511     }
512     counter <- inc(x = counter)
513     run[counter] <- i
514   }
515   consecutive[[index]] <- run[1:counter]
516   consecutive <- consecutive[1:index]
517   return(consecutive)
518 }
520 # Calculate nUMI and nGene
521 #
522 # @param object An object
523 # @param chunk.size Number of cells to chunk over
524 # @param is.expr Expression threshold for 'detected' gene. For most datasets, particularly UMI datasets, will be set to 0 (default).
525 # 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).
526 # @param display.progres Display a progress bar?
527 #
528 # @return Stores nUMI in \code{col_attrs/nUMI}; stores nGene in \code{col_attrs/nGene}; will overwrite any existing datasets at these locations
529 #
530 #' @importFrom Matrix rowSums
531 #
532 calc.umi <- function(
533   object,
534   chunk.size = 1000,
535   is.expr = 0,
536   display.progress = TRUE
537 ) {
538   if (display.progress) {
539     cate("Calculating number of UMIs per cell")
540   }
541   object\$apply(
542     name = 'col_attrs/nUMI',
543     FUN = rowSums,
544     MARGIN = 2,
545     chunk.size = chunk.size,
546     dataset.use = 'matrix',
547     overwrite = TRUE,
548     display.progress = display.progress
549   )
550   if (display.progress) {
551     cate("Calculating number of genes expressed per cell")
552     cate("Using a threshold of", is.expr, "for gene expression")
553   }
554   object\$apply(
555     name = 'col_attrs/nGene',
556     FUN = function(mat, is.expr) {
557       return(rowSums(x = mat > is.expr))
558     },
559     MARGIN = 2,
560     chunk.size = chunk.size,
561     dataset.use = 'matrix',
562     overwrite = TRUE,
563     display.progress = display.progress,
564     is.expr = is.expr
565   )
566   invisible(x = object)
567 }
569 # Cat with a new line
570 #
571 # @param ... Text to be output
572 #
573 catn <- function(...) {
574   x = list(...)
575   if (length(x = x)) {
576     if (!is.null(x = names(x = x)) && length(x = x) == 1 && names(x = x) == 'file') {
577       cat(...)
578     } else {
579       cat(..., '\n')
580     }
581   } else {
582     cat()
583   }
584 }
586 # Cat to stderr
587 #
588 # @param ... Text to be output
589 #
590 cate <- function(...) {
591   catn(..., file = stderr())
592 }