Expand apply, update documentation
[loomr.git] / R / loom.R
1 #' @import hdf5r
2 #' @importFrom R6 R6Class
3 NULL
4
5 #' A class for loom files
6 #'
7 #' @docType class
8 #' @name loom-class
9 #' @rdname loom-class
10 #' @aliases loom
11 #' @return Object of \code{\link{R6::R6Class}} to generate \code{loom} objects
12 #' @format An \code{\link{R6::R6Class}} object
13 #' @seealso \code{\link{hdf5r::H5File}}
14 #'
15 #' @field version Version of loomR object was created under
16 #' @field shape Shape of \code{/matrix} in columns (cells) by rows (genes)
17 #' @field chunksize Chunks set for this dataset in columns (cells) by rows (genes)
18 #' @field matrix The main data matrix, stored as columns (cells) by rows (genes)
19 #' @field layers Additional data matricies, the same shape as \code{/matrix}
20 #' @field col.attrs Extra information about cells
21 #' @field row.attrs Extra information about genes
22 #'
23 #' @section Methods:
24 #' \describe{
25 #'   \item{\code{add.layer(layer)}}{Add a data layer to this loom file, must be in column (cells) by row (genes) orientation}
26 #'   \item{\code{add.attribute(attribute, MARGIN)}}{
27 #'     Add extra information to this loom file where
28 #'     \code{attribute} is a named list where each element is a vector that is as long as one dimension of \code{/matrix} and
29 #'     \code{MARGIN} is either 1 for cells or 2 for genes
30 #'   }
31 #'   \item{\code{add.row.attribute(attribute)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 2)}}
32 #'   \item{\code{add.col.attribute(attribute)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 1)}}
33 #'   \item{\code{add.meta.data(meta.data)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 1)}}
34 #'   \item{\code{batch.scan(chunk.size, MARGIN, index.use, dataset.use, force.reset)}, \code{batch.next()}}{
35 #'     Scan a dataset in the loom file from \code{index.use[1]} to \code{index.use[2]}, iterating by \code{chunk.size}.
36 #'     \code{dataset.use} can be the name, not \code{group/name}, unless the name is present in multiple groups.
37 #'     Pass \code{MARGIN = 1} to iterate on cells or \code{MARGIN = 2} to iterate on genes for 'matrix' or any dataset in 'layers'.
38 #'     To force reevaluation of the iterator object, pass \code{force.reset = TRUE}.
39 #'     \code{MARGIN} does not need to be set for datasets in 'row_attrs' or 'col_attrs'.
40 #'     \code{chunk.size} defaults to \code{self$chunksize}, \code{MARGIN} defaults to 1,
41 #'     \code{index.use} defaults to \code{1:self$shape[MARGIN]}, \code{dataset.use} defaults to 'matrix'
42 #'   }
43 #'   \item{\code{apply(name, FUN, MARGIN, chunk.size, index.use, dataset.use, display.progress, expected, ...)}}{...}
44 #'   \item{\code{map(FUN, MARGIN, chunk.size, index.use, dataset.use, display.progress, expected, ...)}}{...}
45 #' }
46 #'
47 #' @importFrom iterators nextElem
48 #' @importFrom itertools hasNext ihasNext ichunk
49 #' @importFrom utils packageVersion txtProgressBar setTxtProgressBar
50 #'
51 #' @export
52 #'
53 loom <- R6Class(
54   # The loom class
55   # Based on the H5File class from hdf5r
56   # Not clonable (no loom$clone method), doesn't make sense since all data is on disk, not in memory
57   # Yes to portability, other packages can subclass the loom class
58   # Class is locked, other fields and methods cannot be added
59   classname = 'loom',
60   inherit = hdf5r::H5File,
61   cloneable = FALSE,
62   portable = TRUE,
63   lock_class = TRUE,
64   # Public fields and methods
65   # See above for documentation
66   public = list(
67     # Fields
68     version = NULL,
69     shape = NULL,
70     chunksize = NULL,
71     matrix = NULL,
72     layers = NULL,
73     col.attrs = NULL,
74     row.attrs = NULL,
75     # Methods
76     initialize = function(filename = NULL, mode = c('a', 'r', 'r+', 'w', 'w-'), ...) {
77       # If the file exists, run validation steps
78       do.validate <- file.exists(filename) && !(mode %in% c('w', 'w+'))
79       super$initialize(filename = filename, mode = mode, ...)
80       if (do.validate) {
81         # Run the validation steps
82         validateLoom(object = self)
83         # Store /matrix and the shape of /matrix
84         self$matrix <- self[['matrix']]
85         self$shape <- self[['matrix']]$dims
86         # Store the chunk size
87         chunks <- h5attr(x = self, which = 'chunks')
88         chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
89         chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
90         chunks <- unlist(x = strsplit(x = chunks, split = ','))
91         self$chunksize <- as.integer(x = chunks)
92         # Store version information
93         self$version <- as.character(x = tryCatch(
94           # Try getting a version
95           # If it doesn't exist, can we write to the file?
96           # If so, store the version as this version of loomR
97           # Otherwise, store the version as NA_character_
98           expr = h5attr(x = self, which = 'version'),
99           error = function(e) {
100             if (mode != 'r') {
101               version <- packageVersion(pkg = 'loomR')
102               h5attr(x = self, which = 'version') <- version
103             } else {
104               version <- NA_character_
105             }
106             return(version)
107           }
108         ))
109         # Load layers
110         private$load_layers()
111         # Load attributes
112         private$load_attributes(MARGIN = 1) # Cells (col_attrs)
113         private$load_attributes(MARGIN = 2) # Genes (row_attrs)
114       } else {
115         # Assume new HDF5 file
116         self$version <- as.character(x = packageVersion(pkg = 'loomR'))
117       }
118     },
119     add.layer = function(layers) {
120       # Value checking
121       if (!is.list(x = layers) || is.null(x = names(x = layers))) {
122         stop("'layers' must be a named list")
123       }
124       if (is.null(x = self$shape)) {
125         stop(private$err_msg)
126       }
127       # Add layers
128       for (i in 1:length(x = layers)) {
129         if (!is.matrix(x = layers[[i]])) {
130           layers[[i]] <- as.matrix(x = layers[[i]])
131         }
132         do.transpose <- FALSE
133         if (any(dim(x = layers[[i]]) != self$shape)) {
134           if (all(rev(x = dim(x = layers[[i]])) == self$shape)) {
135             do.transpose <- TRUE
136           } else {
137             stop(paste(
138               "All layers must have",
139               self$shape[1],
140               "rows for cells and",
141               self$shape[2],
142               "columns for genes"
143             ))
144           }
145         }
146         if (do.transpose) {
147           self[['layers', names(x = layers)[i]]] <- t(x = layers[[i]])
148         } else {
149           self[['layers', names(x = layers)[i]]] <- layers[[i]]
150         }
151       }
152       self$flush()
153       private$load_layers()
154       invisible(x = self)
155     },
156     add.attribute = function(attribute, MARGIN) {
157       # Value checking
158       if (!is.list(x = attribute) || is.null(x = names(x = attribute))) {
159         stop("'attribute' must be a named list")
160       }
161       if (length(x = attribute) > 1) {
162         for (i in attribute) {
163           if (!is.vector(x = attribute)) {
164             stop("All attributes must be one-dimensional vectors")
165           }
166         }
167       }
168       if (length(x = which(x = names(x = attribute) != '')) != length(x = attribute)) {
169         stop("Not all attributes had names provided")
170       }
171       if (!MARGIN %in% c(1, 2)) {
172         stop("'MARGIN' must be 1 or 2")
173       }
174       # Add the attributes as datasets for our MARGIN's group
175       if (is.null(x = self$shape)) {
176         stop(private$err_msg)
177       }
178       grp.name <- c('col_attrs', 'row_attrs')[MARGIN]
179       grp <- self[[grp.name]]
180       for (i in 1:length(x = attribute)) {
181         if (length(attribute[[i]]) != self$shape[MARGIN])
182           stop(paste(
183             "All",
184             switch(EXPR = MARGIN, '1' = 'cell', '2' = 'gene'),
185             "attributes must be of length",
186             self$shape[MARGIN]
187           ))
188         grp[[names(x = attribute)[i]]] <- attribute[[i]]
189       }
190       self$flush()
191       gc(verbose = FALSE)
192       # Load the attributes for this margin
193       private$load_attributes(MARGIN = MARGIN)
194       invisible(x = self)
195     },
196     # Add attributes for genes
197     add.row.attribute = function(attribute) {
198       self$add.attribute(attribute = attribute, MARGIN = 2)
199       invisible(x = self)
200     },
201     # Add attributes for cells
202     add.col.attribute = function(attribute) {
203       self$add.attribute(attribute = attribute, MARGIN = 1)
204       invisible(x = self)
205     },
206     # Add metadata, follows cells
207     add.meta.data = function(meta.data) {
208       self$add.col.attribute(attribute = meta.data)
209       invisible(x = self)
210     },
211     # Batch scan
212     batch.scan = function(
213       chunk.size = NULL,
214       MARGIN = 1,
215       index.use = NULL,
216       dataset.use = 'matrix',
217       force.reset = FALSE
218     ) {
219       if (is.null(x = private$it) || !grepl(pattern = dataset.use, x = private$iter.dataset) || force.reset) {
220         # Check the existence of the dataset
221         private$iter.dataset <- grep(
222           pattern = dataset.use,
223           x = list.datasets(object = self),
224           value = TRUE
225         )
226         if (length(x = private$iter.dataset) != 1) {
227           stop(paste0("Cannot find dataset '", dataset.use, "' in the loom file"))
228         }
229         if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
230           MARGIN <- 1
231         } else if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
232           MARGIN <- 2
233         }
234         # Check the margin
235         if (!(MARGIN %in% c(1, 2))) {
236           stop("MARGIN must be 1 (cells) or 2 (genes)")
237         } else {
238           private$iter.margin <- MARGIN
239         }
240         if (is.null(x = chunk.size)) {
241           chunk.size <- self$chunksize[MARGIN]
242         }
243         # Set the indices to use
244         index.use <- private$iter_range(index.use = index.use)
245         # Setup our iterator
246         private$it <- ihasNext(iterable = ichunk(
247           iterable = index.use[1]:index.use[2],
248           chunkSize = chunk.size
249         ))
250         private$iter.index <- c(index.use[1], ceiling(x = index.use[2] / chunk.size))
251       }
252       # Return the times we iterate, this isn't important, we only need the length of this vector
253       return(private$iter.index[1]:private$iter.index[2])
254     },
255     batch.next = function(return.data = TRUE) {
256       # Ensure that we have a proper iterator
257       if (!'hasNext.ihasNext' %in% suppressWarnings(expr = methods(class = class(x = private$it)))) {
258         stop("Please setup the iterator with self$batch.scan")
259       }
260       # Do the iterating
261       if (hasNext(obj = private$it)) {
262         # Get the indices for this chunk
263         chunk.indices <- unlist(x = nextElem(obj = private$it))
264         if (return.data) {
265           # If we're working with a matrix dataset, ensure chunking on proper axis
266           if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
267             to.return <- switch(
268               EXPR = private$iter.margin,
269               '1' = self[[private$iter.dataset]][chunk.indices, ],
270               '2' = self[[private$iter.dataset]][, chunk.indices]
271             )
272           } else {
273             # Otherwise, iterating over an attribute (1 dimensional)
274             to.return <- self[[private$iter.dataset]][chunk.indices]
275           }
276         } else {
277           to.return <- chunk.indices
278         }
279         # Determine if we reset the iterator
280         if (!hasNext(obj = private$it)) {
281           private$reset_batch()
282         }
283         return(to.return)
284       } else {
285         # Just in case
286         private$reset_batch()
287         return(NULL)
288       }
289     },
290     apply = function(
291       name,
292       FUN,
293       MARGIN = 1,
294       chunk.size = NULL,
295       index.use = NULL,
296       dataset.use = 'matrix',
297       display.progress = TRUE,
298       expected = NULL,
299       ...
300     ) {
301       if (!inherits(x = FUN, what = 'function')) {
302         stop("FUN must be a function")
303       }
304       # Checks datset, index, and MARGIN
305       # Sets full dataset path in private$iter.dataset
306       # Sets proper MARGIN in private$iter.margin
307       batch <- self$batch.scan(
308         chunk.size = chunk.size,
309         MARGIN = MARGIN,
310         index.use = index.use,
311         dataset.use = dataset.use,
312         force.reset = TRUE
313       )
314       # Check how we store our results
315       # And what the shape of our dataset is
316       results.matrix <- TRUE
317       dataset.matrix <- TRUE
318       if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
319         results.matrix <- FALSE
320         dataset.matrix <- FALSE
321       } else if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
322         results.matrix <- FALSE
323         dataset.matrix <- FALSE
324       }
325       if (!is.null(x = expected)) {
326         results.matrix <- switch(
327           EXPR = expected,
328           'vector' = FALSE,
329           'matrix' = TRUE,
330           stop("'expected' must be one of 'matrix', 'vector', or NULL")
331         )
332       }
333       if (display.progress) {
334         pb <- txtProgressBar(char = '=', style = 3)
335       }
336       for (i in length(x = batch)) {
337         if (display.progress) {
338           setTxtProgressBar(pb = pb, value = i / length(x = batch))
339         }
340       }
341       private$reset_batch()
342       invisible(x = NULL)
343     },
344     map = function(
345       FUN,
346       MARGIN = 1,
347       chunk.size = NULL,
348       index.use = NULL,
349       dataset.use = 'matrix',
350       display.progress = TRUE,
351       expected = NULL,
352       ...
353     ) {
354       if (!inherits(x = FUN, what = 'function')) {
355         stop("FUN must be a function")
356       }
357       # Checks datset, index, and MARGIN
358       # Sets full dataset path in private$iter.dataset
359       # Sets proper MARGIN in private$iter.margin
360       batch <- self$batch.scan(
361         chunk.size = chunk.size,
362         MARGIN = MARGIN,
363         index.use = index.use,
364         dataset.use = dataset.use,
365         force.reset = TRUE
366       )
367       # Check how we store our results
368       # And what the shape of our dataset is
369       results.matrix <- TRUE
370       dataset.matrix <- TRUE
371       if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
372         results.matrix <- FALSE
373         dataset.matrix <- FALSE
374       } else if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
375         results.matrix <- FALSE
376         dataset.matrix <- FALSE
377       }
378       if (!is.null(x = expected)) {
379         results.matrix <- switch(
380           EXPR = expected,
381           'vector' = FALSE,
382           'matrix' = TRUE,
383           stop("'expected' must be one of 'matrix', 'vector', or NULL")
384         )
385       }
386       # Determine the shape of our results
387       index.use <- private$iter_range(index.use = index.use)
388       # Create our results holding object
389       if (results.matrix) {
390         switch(
391           EXPR = private$iter.margin,
392           '1' = results <- matrix(
393             nrow = length(x = index.use[1]:index.use[2]),
394             ncol = self$shape[2]
395           ),
396           '2' = results <- matrix(
397             nrow = self$shape[1],
398             ncol = length(x = index.use[1]:index.use[2])
399           )
400         )
401       } else {
402         results <- vector(length = length(x = index.use[1]:index.use[2]))
403       }
404       if (display.progress) {
405         pb <- txtProgressBar(char = '=', style = 3)
406       }
407       for (i in 1:length(x = batch)) {
408         chunk.indices <- self$batch.next(return.data = FALSE)
409         chunk.data <- if (dataset.matrix) {
410           switch(
411             EXPR = MARGIN,
412             '1' = self[[dataset.use]][chunk.indices, ],
413             '2' = self[[dataset.use]][, chunk.indices]
414           )
415         } else {
416           self[[dataset.use]][chunk.indices]
417         }
418         chunk.data <- FUN(chunk.data, ...)
419         if (results.matrix) {
420           if (MARGIN == 1) {
421             results[chunk.indices, ] <- chunk.data
422           } else if (MARGIN == 2) {
423             results[, chunk.indices] <- chunk.data
424           }
425         } else {
426           results[chunk.indices] <- chunk.data
427         }
428         if (display.progress) {
429           setTxtProgressBar(pb = pb, value = i / length(x = batch))
430         }
431       }
432       private$reset_batch()
433       return(results)
434     }
435   ),
436   # Private fields and methods
437   # @field err_msg A simple error message if this object hasn't been created with loomR::create or loomR::connect
438   # @field it Iterator object for batch.scan and batch.next
439   # @field iter.dataset Dataset for iterating on
440   # @field iter.margin Margin to iterate over
441   # @field iter.index # Index values for iteration
442   # \describe{
443   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
444   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
445   #   \item{\code{reset_batch()}}{Reset the batch iterator fields}
446   #   \item{\code{iter_range(index.use)}}{Get the range of indices for a batch iteration}
447   # }
448   private = list(
449     # Fields
450     err_msg = "This loom object has not been created with either loomR::create or loomR::connect, please use these functions to create or connect to a loom file",
451     it = NULL,
452     iter.dataset = NULL,
453     iter.margin = NULL,
454     iter.index = NULL,
455     # Methods
456     load_attributes = function(MARGIN) {
457       attribute <- switch(
458         EXPR = MARGIN,
459         '1' = 'col_attrs',
460         '2' = 'row_attrs',
461         stop('Invalid attribute dimension')
462       )
463       group <- self[[attribute]]
464       attributes <- unlist(x = lapply(
465         X = names(x = group),
466         FUN = function(x) {
467           d <- list(group[[x]])
468           names(x = d) <- x
469           return(d)
470         }
471       ))
472       if (MARGIN == 1) {
473         self$col.attrs <- attributes
474       } else if (MARGIN == 2) {
475         self$row.attrs <- attributes
476       }
477     },
478     load_layers = function() {
479       self$layers <- unlist(x = lapply(
480         X = names(x = self[['layers']]),
481         FUN = function(n) {
482           d <- c(self[['layers', n]])
483           names(x = d) <- n
484           return(d)
485         }
486       ))
487     },
488     reset_batch = function() {
489       private$it <- NULL
490       private$iter.dataset <- NULL
491       private$iter.margin <- NULL
492       private$iter.index <- NULL
493     },
494     iter_range = function(index.use) {
495       if (is.null(private$iter.margin)) {
496         stop("Batch processing hasn't been set up")
497       }
498       if (is.null(x = index.use)) {
499         # If no index was provided, use entire range for this margin
500         index.use <- c(1, self$shape[private$iter.margin])
501       } else if (length(x = index.use) == 1) {
502         # If one index was provided, start at one and go to index
503         index.use <- c(1, index.use)
504       } else {
505         # Use index.use[1] and index.use[2]
506         index.use <- c(index.use[1], index.use[2])
507       }
508       # Ensure the indices provided fit within the range of the dataset
509       index.use[1] <- max(1, index.use[1])
510       index.use[2] <- min(index.use[2], self$shape[private$iter.margin])
511       # Ensure that index.use[1] is greater than index.use[2]
512       if (index.use[1] > index.use[2]) {
513         stop(paste0(
514           "Starting index (",
515           index.use[1],
516           ") must be lower than the ending index (",
517           index.use[2],
518           ")"
519         ))
520       }
521       return(index.use)
522     }
523   )
524 )
525
526 #' Create a loom object
527 #'
528 #' @param filename The name of the new loom file
529 #' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
530 #' @param gene.attrs A named list of vectors with extra data for genes, each vector must be as long as the number of genes in \code{data}
531 #' @param cell.attrs A named list of vectors with extra data for cells, each vector must be as long as the number of cells in \code{data}
532 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
533 #'
534 #' @return A connection to a loom file
535 #'
536 #' @importFrom utils packageVersion
537 #'
538 #' @seealso \code{\link{loom-class}}
539 #'
540 #' @export
541 #'
542 create <- function(
543   filename,
544   data,
545   gene.attrs = NULL,
546   cell.attrs = NULL,
547   layers = NULL,
548   chunk.dims = 'auto'
549 ) {
550   if (file.exists(filename)) {
551     stop(paste('File', file, 'already exists!'))
552   }
553   if (!is.matrix(x = data)) {
554     data <- as.matrix(x = data)
555   }
556   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
557     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
558   } else if (length(x = chunk.dims == 1)) {
559     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
560       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
561     }
562   } else {
563     chunk.dims <- as.integer(x = chunk.dims)
564   }
565   new.loom <- loom$new(filename = filename, mode = 'w-')
566   # Create the matrix
567   new.loom$create_dataset(
568     name = 'matrix',
569     robj = data,
570     chunk_dims = chunk.dims
571   )
572   new.loom$matrix <- new.loom[['matrix']]
573   new.loom$shape <- new.loom[['matrix']]$dims
574   # Groups
575   new.loom$create_group(name = 'layers')
576   new.loom$create_group(name = 'row_attrs')
577   new.loom$create_group(name = 'col_attrs')
578   # Check for the existance of gene or cell names
579   if (!is.null(x = colnames(x = data))) {
580     new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
581   }
582   if (!is.null(x = rownames(x = data))) {
583     new.loom$add.col.attribute(attribute = list('cell_names' = rownames(x = data)))
584   }
585   # Store some constants as HDF5 attributes
586   h5attr(x = new.loom, which = 'version') <- new.loom$version
587   h5attr(x = new.loom, which = 'chunks') <- paste0(
588     '(',
589     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
590     ')'
591   )
592   # Add layers
593   if (!is.null(x = layers)) {
594     new.loom$add.layer(layer = layers)
595   }
596   if (!is.null(x = gene.attrs)) {
597     new.loom$add.row.attribute(attribute = gene.attrs)
598   }
599   if (!is.null(x = cell.attrs)) {
600     new.loom$add.col.attribute(attribute = cell.attrs)
601   }
602   # Set last bit of information
603   chunks <- new.loom[['matrix']]$chunk_dims
604   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
605   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
606   chunks <- unlist(x = strsplit(x = chunks, split = ','))
607   new.loom$chunksize <- as.integer(x = chunks)
608   # Return the connection
609   return(new.loom)
610 }
611
612 #' Validate a loom object
613 #'
614 #' @param object A loom object
615 #'
616 #' @return None, errors out if object is an invalid loom connection
617 #'
618 #' @seealso \code{\link{loom-class}}
619 #'
620 #' @export
621 #'
622 validateLoom <- function(object) {
623   # A loom file is a specific HDF5
624   # We need a dataset in /matrix that's a two-dimensional dense matrix
625   root.datasets <- list.datasets(object = object, path = '/', recursive = FALSE)
626   if (length(x = root.datasets) != 1) {
627     stop("There can only be one dataset at the root of the loom file")
628   }
629   if (root.datasets != 'matrix') {
630     stop("The root dataset must be called '/matrix'")
631   }
632   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
633   required.groups <- c('row_attrs', 'col_attrs', 'layers')
634   dim.matrix <- object[[root.datasets]]$dims # Columns x Rows
635   names(x = dim.matrix) <- required.groups[c(2, 1)]
636   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
637   group.msg <- paste0(
638     "There can only be three groups in the loom file: '",
639     paste(required.groups, collapse = "', '"),
640     "'"
641   )
642   if (length(x = root.groups) != 3) {
643     stop(group.msg)
644   }
645   if (!all(required.groups %in% root.groups)) {
646     stop(group.msg)
647   }
648   unlist(x = sapply(
649     X = required.groups[1:2],
650     FUN = function(group) {
651       if (length(x = list.groups(object = object[[group]], recursive = FALSE)) > 0) {
652         stop(paste("Group", group, "cannot have subgroups"))
653       }
654       if (length(x = list.attributes(object = object[[group]])) > 0) {
655         stop(paste("Group", group, "cannot have subattributes"))
656       }
657       for (dataset in list.datasets(object = object[[group]])) {
658         if (object[[paste(group, dataset, sep = '/')]]$dims != dim.matrix[group]) {
659           stop(paste("All datasets in group", group, "must be of length", required.groups[group]))
660         }
661       }
662     }
663   ))
664   for (dataset in list.datasets(object = object[['/layers']])) {
665     if (any(object[[paste('layers', dataset, sep = '/')]]$dims != dim.matrix)) {
666       stop(paste("All datasets in '/layers' must be", dim.matrix[1], 'by', dim.matrix[2]))
667     }
668   }
669 }
670
671 #' Connect to a loom file
672 #'
673 #' @param filename The loom file to connect to
674 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write
675 #'
676 #' @return A loom file connection
677 #'
678 #' @seealso \code{\link{loom-class}}
679 #'
680 #' @export
681 #'
682 connect <- function(filename, mode = "r") {
683   if (!(mode %in% c('r', 'r+'))) {
684     stop("'mode' must be one of 'r' or 'r+'")
685   }
686   new.loom <- loom$new(filename = filename, mode = mode)
687   return(new.loom)
688 }
689
690 CreateLoomFromSeurat <- function(object, filename) {
691   object.data=t(object@raw.data[rownames(object@data),object@cell.names])
692   object.meta.data=object@meta.data
693   row_attrs=list(); col_attrs=list()
694   gene.names=colnames(object.data)
695   object.meta.data$ident = object@ident
696   object.meta.data$CellID = object@cell.names
697   for(i in 1:ncol(object.meta.data)) {
698     col_attrs[[colnames(object.meta.data)[i]]]=object.meta.data[,i]
699   }
700   row_attrs[["Gene"]]=gene.names
701   create(filename,object.data,gene.attrs = row_attrs, cell.attrs = col_attrs)
702 }
703
704 #' Map a function or a series of functions over a loom file
705 #'
706 #' @param X A loom object
707 #' @param MARGIN The dimmension to map over, pass 1 for cells or 2 for genes
708 #' @param FUN A function to map to the loom file
709 #' @param chunk.size Chunk size to use, defaults to \code{loomfile$chunksize[MARGIN]}
710 #' @param index.use Indices of the dataset to use, defaults to \code{1:loomfile$shape[MARGIN]}
711 #' @param dataset.use Dataset to use, defauts to 'matrix'
712 #' @param display.progress Display a progress bar
713 #' @param expected Shape of expected results. Can pass either 'matrix' or 'vector'; defaults to shape of 'dataset.use'
714 #' @param ... Extra parameters for FUN
715 #'
716 #' @return The results of the map
717 #'
718 #' @importFrom utils txtProgressBar setTxtProgressBar
719 #'
720 #' @export
721 #'
722 map <- function(
723   X,
724   MARGIN = 1,
725   FUN,
726   chunk.size = NULL,
727   index.use = NULL,
728   dataset.use = 'matrix',
729   display.progress = TRUE,
730   expected = NULL,
731   ...
732 ) {
733   if (!inherits(x = X, what = 'loom')) {
734     stop("map only works on loom objects")
735   }
736   if (!inherits(x = FUN, what = 'function')) {
737     stop("FUN must be a function")
738   }
739   # Check for existance of dataset
740   if (!any(grepl(pattern = dataset.use, x = list.datasets(object = X)))) {
741     stop(paste("Cannot find dataset", dataset.use, "in the loom file"))
742   }
743   # Figure out if we're returning a vector or matrix
744   full.dataset <- grep(
745     pattern = dataset.use,
746     x = list.datasets(object = X),
747     value = TRUE
748   )
749   results.matrix <- TRUE
750   dataset.matrix <- TRUE
751   if (grepl(pattern = 'col_attrs', x = full.dataset)) {
752     MARGIN <- 1
753     results.matrix <- FALSE
754     dataset.matrix <- FALSE
755   } else if (grepl(pattern = 'row_attrs', x = full.dataset)) {
756     MARGIN <- 2
757     results.matrix <- FALSE
758     dataset.matrix <- FALSE
759   }
760   if (!is.null(x = expected)) {
761     results.matrix <- switch(
762       EXPR = expected,
763       'vector' = FALSE,
764       'matrix' = TRUE,
765       stop("'expected' must be one of 'matrix', 'vector', or NULL")
766     )
767   }
768   # Determine the shape of our results
769   if (!(MARGIN %in% c(1, 2))) {
770     stop("MARGIN must be either 1 (cells) or 2 (genes)")
771   }
772   if (is.null(x = index.use)) {
773     index.use <- c(1, X$shape[MARGIN])
774   } else if (length(x = index.use) == 1) {
775     index.use <- c(1, index.use)
776   }
777   index.use[1] <- max(1, index.use[1])
778   index.use[2] <- min(index.use[2], X$shape[MARGIN])
779   batch <- X$batch.scan(
780     chunk.size = chunk.size,
781     MARGIN = MARGIN,
782     index.use = index.use,
783     dataset.use = dataset.use,
784     force.reset = TRUE
785   )
786   # Create our results holding object
787   if (results.matrix) {
788     switch(
789       EXPR = MARGIN,
790       '1' = results <- matrix(
791         nrow = length(x = index.use[1]:index.use[2]),
792         ncol = X$shape[2]
793       ),
794       '2' = results <- matrix(
795         nrow = X$shape[1],
796         ncol = length(x = index.use[1]:index.use[2])
797       )
798     )
799   } else {
800     results <- vector(length = length(x = index.use[1]:index.use[2]))
801   }
802   if (display.progress) {
803     pb <- txtProgressBar(char = '=', style = 3)
804   }
805   for (i in 1:length(x = batch)) {
806     chunk.indices <- X$batch.next(return.data = FALSE)
807     chunk.data <- if (dataset.matrix) {
808       switch(
809         EXPR = MARGIN,
810         '1' = X[[dataset.use]][chunk.indices, ],
811         '2' = X[[dataset.use]][, chunk.indices]
812       )
813     } else {
814       X[[dataset.use]][chunk.indices]
815     }
816     chunk.data <- FUN(chunk.data, ...)
817     if (results.matrix) {
818       if (MARGIN == 1) {
819         results[chunk.indices, ] <- chunk.data
820       } else if (MARGIN == 2) {
821         results[, chunk.indices] <- chunk.data
822       }
823     } else {
824       results[chunk.indices] <- chunk.data
825     }
826     if (display.progress) {
827       setTxtProgressBar(pb = pb, value = i / length(x = batch))
828     }
829   }
830   return(results)
831 }
832
833 #need to comment
834 #need to add progress bar
835 #but otherwise, pretty cool
836 #for paul to try :
837 # f <- connect("~/Downloads/10X43_1.loom")
838 # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
839 # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
840 map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
841   n_func = length(f_list)
842   if (n_func == 1) {
843     f_list <- list(f_list)
844   }
845   if (MARGIN == 1) {
846     results <- list()
847     for (j in 1:n_func) {
848       results[[j]] <- numeric(0)
849     }
850     rows_per_chunk <- chunksize
851     ix <- 1
852     while (ix <= self@shape[1]) {
853       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
854       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
855       for (j in 1:n_func) {
856         new_results <- apply(chunk, 1, FUN = f_list[[j]])
857         results[[j]] <- c(results[[j]], new_results)
858       }
859       ix <- ix + chunksize
860     }
861   }
862   if (MARGIN == 2) {
863     results <- list()
864     for (j in 1:n_func) {
865       results[[j]] <- numeric(0)
866     }
867     cols_per_chunk <- chunksize
868     ix <- 1
869     while (ix <= self@shape[2]) {
870       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
871       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
872       for (j in 1:n_func) {
873         new_results <- apply(chunk, 2, FUN = f_list[[j]])
874         results[[j]] <- c(results[[j]], new_results)
875       }
876       ix <- ix + chunksize
877     }
878   }
879   if (n_func == 1) {
880     results <- results[[1]]
881   }
882   return(results)
883 }