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