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