f4e0c20412c99766a3e0057229de03b2de3ec197
[loomr.git] / R / loom.R
1 #' @include internal.R
2 #' @import hdf5r
3 #' @importFrom R6 R6Class
4 NULL
5
6 #' A class for loom files
7 #'
8 #' @docType class
9 #' @name loom-class
10 #' @rdname loom-class
11 #' @aliases loom
12 #' @return Object of \code{\link{R6::R6Class}} to generate \code{loom} objects
13 #' @format An \code{\link{R6::R6Class}} object
14 #' @seealso \code{\link{hdf5r::H5File}}
15 #'
16 #' @field version Version of loomR object was created under
17 #' @field shape Shape of \code{/matrix} in genes (columns) by cells (rows)
18 #' @field chunksize Chunks set for this dataset in columns (cells) by rows (genes)
19 #' @field matrix The main data matrix, stored as columns (cells) by rows (genes)
20 #' @field layers Additional data matricies, the same shape as \code{/matrix}
21 #' @field col.attrs Extra information about cells
22 #' @field row.attrs Extra information about genes
23 #'
24 #' @section Methods:
25 #' \describe{
26 #'   \item{\code{add.layer(layer, chunk.size, overwrite)}}{
27 #'     Add a data layer to this loom file, must be the same dimensions as \code{/matrix}
28 #'     \describe{
29 #'       \item{\code{layer}}{A named list of matrices to be added as layers}
30 #'       \item{\code{chunk.size}}{Number of rows from each layer to stream at once, defaults to 1000}
31 #'       \item{\code{overwrite}}{If a layer already exists, overwrite with new data, defaults to FALSE}
32 #'     }
33 #'   }
34 #'   \item{\code{add.attribute(attribute, MARGIN, overwrite)}}{
35 #'     Add extra information to this loom file.
36 #'     \describe{
37 #'       \item{\code{attribute}}{A named list where the first dimmension of each element as long as one dimension of \code{/matrix}}
38 #'       \item{\code{MARGIN}}{Either 1 for genes or 2 for cells}
39 #'       \item{\code{overwrite}}{Can overwrite existing attributes?}
40 #'     }
41 #'   }
42 #'   \item{\code{add.row.attribute(attribute), add.col.attribute(attribute)}}{
43 #'     Add row or column attributes
44 #'   }
45 #'   \item{\code{batch.scan(chunk.size, MARGIN, index.use, dataset.use, force.reset)}, \code{batch.next(return.data)}}{
46 #'     Scan a dataset in the loom file from \code{index.use[1]} to \code{index.use[2]}, iterating by \code{chunk.size}.
47 #'     \code{dataset.use} can be the name, not \code{group/name}, unless the name is present in multiple groups.
48 #'     Pass \code{MARGIN = 1} to iterate on genes or \code{MARGIN = 2} to iterate on cells for 'matrix' or any dataset in 'layers'.
49 #'     To force reevaluation of the iterator object, pass \code{force.reset = TRUE}.
50 #'     \code{MARGIN} does not need to be set for datasets in 'row_attrs' or 'col_attrs'.
51 #'     \code{chunk.size} defaults to \code{self$chunksize}, \code{MARGIN} defaults to 2,
52 #'     \code{index.use} defaults to \code{1:self$shape[MARGIN]}, \code{dataset.use} defaults to 'matrix'
53 #'   }
54 #'   \item{\code{apply(name, FUN, MARGIN, chunk.size, dataset.use, overwrite, display.progress, ...)}}{
55 #'     Apply a function over a dataset within the loom file, stores the results in the loom file.
56 #'     \code{name} must be the full name of the dataset ('group/name').
57 #'     \code{apply} will always use the entire dataset specified in \code{dataset.use}
58 #'     \describe{
59 #'       \item{\code{name}}{Name of dataset to store results to}
60 #'       \item{\code{FUN}}{Function to apply}
61 #'       \item{\code{MARGIN}}{Iterate over genes (1) or cells (2)}
62 #'       \item{\code{chunk.size}}{Size to chunk \code{MARGIN} by}
63 #'       \item{\code{dataset.use}}{Name of dataset to use}
64 #'       \item{\code{overwrite}}{Overite \code{name} if already exists}
65 #'       \item{\code{display.progress}}{Display progress}
66 #'       \item{\code{...}}{Extra parameters to pass to \code{FUN}}
67 #'     }
68 #'   }
69 #'   \item{\code{map(FUN, MARGIN, chunk.size, index.use, dataset.use, display.progress, expected, ...)}}{
70 #'     Map a function onto a dataset within the loom file, returns the result into R.
71 #'     The result will default to the shape of the dataset used; to change pass either 'vector' or 'matrix' to \code{expected}.
72 #'   }
73 #'   \item{\code{add.cells(matrix.data, attributes.data = NULL, layers.data = NULL, display.progress = TRUE)}}{
74 #'     Add m2 cells to a loom file.
75 #'     \describe{
76 #'       \item{\code{matrix.data}}{a list of m2 cells where each entry is a vector of length n (num genes, \code{self$shape[1]})}
77 #'       \item{\code{attributes.data}}{a list where each entry is named for one of the datasets in \code{self[['col_attrs']]}; each entry is a vector of length m2.}
78 #'       \item{\code{layers.data}}{a list where each entry is named for one of the datasets in \code{self[['layers']]}; each entry is an n-by-m2 matrix where n is the number of genes in this loom file and m2 is the number of cells being added.}
79 #'       \item{\code{display.progress}}{display progress}
80 #'     }
81 #'   }
82 #' }
83 #'
84 #' @importFrom iterators nextElem
85 #' @importFrom itertools hasNext ihasNext ichunk
86 #' @importFrom utils packageVersion txtProgressBar setTxtProgressBar
87 #'
88 #' @export
89 #'
90 loom <- R6Class(
91   # The loom class
92   # Based on the H5File class from hdf5r
93   # Not clonable (no loom$clone method), doesn't make sense since all data is on disk, not in memory
94   # Yes to portability, other packages can subclass the loom class
95   # Class is locked, other fields and methods cannot be added
96   classname = 'loom',
97   inherit = hdf5r::H5File,
98   cloneable = FALSE,
99   portable = TRUE,
100   lock_class = TRUE,
101   # Public fields and methods
102   # See above for documentation
103   public = list(
104     # Fields
105     version = NULL,
106     shape = NULL,
107     chunksize = NULL,
108     matrix = NULL,
109     layers = NULL,
110     col.attrs = NULL,
111     row.attrs = NULL,
112     # Methods
113     initialize = function(
114       filename = NULL,
115       mode = c('a', 'r', 'r+', 'w', 'w-'),
116       skip.validate = FALSE,
117       ...
118     ) {
119       # If the file exists, run validation steps
120       do.validate <- file.exists(filename) && !(mode %in% c('w', 'w+'))
121       private$skipped.validation <- skip.validate
122       super$initialize(filename = filename, mode = mode, ...)
123       if (do.validate) {
124         # Run the validation steps
125         if (skip.validate) {
126           warning("Skipping validation step, some fields are not populated")
127         } else {
128           validateLoom(object = self)
129         }
130         # Store /matrix and the shape of /matrix
131         if (skip.validate) {
132           if (getOption(x = 'verbose')) {
133             warning("Not setting matrix field")
134           }
135         } else {
136           self$matrix <- self[['matrix']]
137         }
138         self$shape <- rev(self[['matrix']]$dims)
139         # Store the chunk size
140         chunks <- tryCatch(
141           expr = h5attr(x = self, which = 'chunks'),
142           error = function(e) {
143             hchunks <- self[['matrix']]$chunk_dims
144             pchunks <- paste0('(', paste(hchunks, collapse = ', '), ')')
145             if (mode != 'r') {
146               h5attr(x = self, which = 'chunks') <- pchunks
147             }
148             return(pchunks)
149           }
150         )
151         chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
152         chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
153         chunks <- unlist(x = strsplit(x = chunks, split = ','))
154         self$chunksize <- as.integer(x = chunks)
155         # Store version information
156         self$version <- as.character(x = tryCatch(
157           # Try getting a version
158           # If it doesn't exist, can we write to the file?
159           # If so, store the version as this version of loomR
160           # Otherwise, store the version as NA_character_
161           expr = h5attr(x = self, which = 'version'),
162           error = function(e) {
163             if (mode != 'r') {
164               version <- packageVersion(pkg = 'loomR')
165               h5attr(x = self, which = 'version') <- as.character(x = version)
166             } else {
167               version <- NA_character_
168             }
169             return(version)
170           }
171         ))
172         # Load layers
173         private$load_layers()
174         # Load attributes
175         private$load_attributes(MARGIN = 1) # Genes (row_attrs)
176         private$load_attributes(MARGIN = 2) # Cells (col_attrs
177       } else {
178         # Assume new HDF5 file
179         self$version <- as.character(x = packageVersion(pkg = 'loomR'))
180       }
181     },
182     finalizer = function() {
183       self$close_all(close_self = TRUE)
184     },
185     load.fields = function() {
186       private$load_layers()
187       private$load_attributes(MARGIN = 1)
188       private$load_attributes(MARGIN = 2)
189     },
190     # Addding attributes and layers
191     add.layer = function(layers, chunk.size = 1000, overwrite = FALSE) {
192       if (self$mode == 'r') {
193         stop(private$err_mode)
194       }
195       # Value checking
196       if (!is.list(x = layers) || is.null(x = names(x = layers))) {
197         stop("'layers' must be a named list")
198       }
199       if (is.null(x = self$shape)) {
200         stop(private$err_init)
201       }
202       # Add layers
203       for (i in 1:length(x = layers)) {
204         # if (!is.matrix(x = layers[[i]])) {
205         if (!inherits(x = layers[[i]], what = c('Matrix', 'matrix'))) {
206           layers[[i]] <- as.matrix(x = layers[[i]])
207         }
208         do.transpose <- FALSE
209         shape.use <- rev(x = self$shape)
210         if (any(dim(x = layers[[i]]) != shape.use)) {
211           if (all(rev(x = dim(x = layers[[i]])) == shape.use)) {
212             do.transpose <- TRUE
213           } else {
214             stop(paste(
215               "All layers must have",
216               shape.use[1],
217               "rows for cells and",
218               shape.use[2],
219               "columns for genes"
220             ))
221           }
222         }
223         if (do.transpose) {
224           layers[[i]] <- t(x = layers[[i]])
225         }
226         layer.name <- names(x = layers)[i]
227         if (layer.name %in% list.datasets(object = self[['layers']])) {
228           if (overwrite) {
229             self[['layers']]$link_delete(name = layer.name)
230           } else {
231             stop(paste(
232               "A layer with the name",
233               layer.name,
234               "already!"
235             ))
236           }
237         }
238         dtype <- getDtype(x = layers[[i]][1, 1])
239         self[['layers']]$create_dataset(
240           name = layer.name,
241           dtype = dtype,
242           dims = dim(x = layers[[i]])
243         )
244         chunk.points <- chunkPoints(
245           data.size = dim(x = layers[[i]])[1],
246           chunk.size = chunk.size
247         )
248         # if (display.progress) {
249         #   pb <- txtProgressBar(char = '=', style = 3)
250         # }
251         for (col in 1:ncol(x = chunk.points)) {
252           row.start <- chunk.points[1, col]
253           row.end <- chunk.points[2, col]
254           self[['layers']][[layer.name]][row.start:row.end, ] <- as.matrix(
255             x = layers[[i]][row.start:row.end, ]
256           )
257           # if (display.progress) {
258           #   setTxtProgressBar(pb = pb, value = col / ncol(x = chunk.points))
259           # }
260         }
261         # self[['layers']]$create_dataset(
262         #   name = names(x = layers)[i],
263         #   robj = layers[[i]],
264         #   chunk_dims = self$chunksize
265         # )
266       }
267       self$flush()
268       gc(verbose = FALSE)
269       private$load_layers()
270       invisible(x = self)
271     },
272     add.attribute = function(attribute, MARGIN, overwrite = FALSE) {
273       if (self$mode == 'r') {
274         stop(private$err_mode)
275       }
276       # Value checking
277       if (is.data.frame(x = attribute)) {
278         attribute <- as.list(x = attribute)
279       }
280       is.actual.list <- is.list(x = attribute)
281       if (!is.actual.list || is.null(x = names(x = attribute))) {
282         stop("Attributes must be provided as a named list")
283       }
284       # if (is.data.frame(x = attribute)) {
285       #   attribute <- as.list(x = attribute)
286       # }
287       # if (!is.list(x = attribute) || is.null(x = names(x = attribute))) {
288       #   stop("'attribute' must be a named list")
289       # }
290       if (!MARGIN %in% c(1, 2)) {
291         stop("'MARGIN' must be 1 or 2")
292       }
293       length.use <- rev(x = self[['matrix']]$dims)[MARGIN]
294       dim.msg <- paste(
295         "At least one dimmension for each",
296         switch(EXPR = MARGIN, '1' = 'gene', '2' = 'cell'),
297         "attribute must be",
298         length.use
299       )
300       for (i in 1:length(x = attribute)) {
301         if (is.matrix(x = attribute[[i]]) || is.data.frame(x = attribute[[i]])) {
302           margin.use <- which(x = dim(x = attribute[[i]]) == length.use)
303           if (!length(x = margin.use)) {
304             stop(dim.msg)
305           }
306           margin.use <- margin.use[1]
307           attribute[[i]] <- switch(
308             EXPR = margin.use,
309             '1' = t(x = as.matrix(x = attribute[[i]])),
310             '2' = as.matrix(x = attribute[[i]]),
311             stop("All attributes must be one- or two-dimmensional")
312           )
313         } else {
314           if (length(x = attribute[[i]]) != length.use) {
315             stop(dim.msg)
316           }
317           attribute[[i]] <- as.vector(x = attribute[[i]])
318         }
319       }
320       if (length(x = which(x = names(x = attribute) != '')) != length(x = attribute)) {
321         stop("Not all attributes had names provided")
322       }
323       grp.name <- c('row_attrs', 'col_attrs')[MARGIN]
324       grp <- self[[grp.name]]
325       if (!overwrite) {
326         names.fail <- which(x = names(x = attribute) %in% names(x = grp))
327         if (length(x = names.fail) != 0) {
328           stop(paste0(
329             "The following names are already used for ",
330             switch(EXPR = MARGIN, '1' = 'row', '2' = 'column'),
331             " attributes: '",
332             paste(names(x = attribute)[names.fail], collapse = ''),
333             "'"
334           ))
335         }
336       }
337       # Add the attributes as datasets for our MARGIN's group
338       for (i in 1:length(x = attribute)) {
339         try(expr = grp$link_delete(name = names(x = attribute)[i]), silent = TRUE)
340         grp[[names(x = attribute)[i]]] <- attribute[[i]]
341       }
342       self$flush()
343       gc(verbose = FALSE)
344       # Load the attributes for this margin
345       private$load_attributes(MARGIN = MARGIN)
346       invisible(x = self)
347     },
348     add.row.attribute = function(attribute, overwrite = FALSE) {
349       self$add.attribute(attribute = attribute, MARGIN = 1, overwrite = overwrite)
350       invisible(x = self)
351     },
352     add.col.attribute = function(attribute, overwrite = FALSE) {
353       self$add.attribute(attribute = attribute, MARGIN = 2, overwrite = overwrite)
354       invisible(x = self)
355     },
356     add.meta.data = function(meta.data, overwrite = FALSE) {
357       self$add.col.attribute(attribute = meta.data, overwrite = overwrite)
358       invisible(x = self)
359     },
360     # Get attribute information
361     get.attribute.df = function(
362       attribute.layer = c("row", "col"),
363       attribute.names = NULL,
364       row.names = "gene_names",
365       col.names = "cell_names"
366     ) {
367       # takes multiple row_attrs or col_attrs and combines them into a data.frame,
368       # removing rows or columns that are entirely NAs.
369       #
370       # attribute.layer either "row" for row_attrs or "col" col_attrs
371       # attribute.names name of rows or columns to combine into matrix
372       # row.names either a character vector or name of an element in row_attrs
373       # col.names either a character vector or name of an element in col_attrs
374       if (!attribute.layer %in% c("row", "col")) {
375         stop("Invalid attribute.layer. Please select either 'row' or 'col'.")
376       }
377       attribute.layer <- paste0(attribute.layer, "_attrs")
378       # check that attribute names are present
379       if (!all(attribute.names %in% self[[attribute.layer]]$names)) {
380         invalid.names <- attribute.names[which(!attribute.names %in% self[[attribute.layer]]$names)]
381         stop(paste0("Invalid attribute.names: ", paste0(invalid.names, collapse = ", ")))
382       }
383       if (attribute.layer == "row_attrs") {
384         combined.df <- data.frame(
385           self[[paste0(attribute.layer, "/", attribute.names[1])]][],
386           row.names = self[[paste0(attribute.layer, "/", row.names)]][]
387         )
388       } else {
389         combined.df <- data.frame(
390           self[[paste0(attribute.layer, "/", attribute.names[1])]][],
391           row.names = self[[paste0(attribute.layer, "/", col.names)]][]
392         )
393       }
394       if (length(x = attribute.names) > 1) {
395         for (i in 2:length(x = attribute.names)) {
396           combined.df[, attribute.names[i]] <- self[[paste0(attribute.layer, "/", attribute.names[i])]][]
397         }
398       }
399       colnames(x = combined.df) <- attribute.names
400       # check if any row is all NAs
401       rows.to.remove <- unname(obj = which(x = apply(
402         X = combined.df,
403         MARGIN = 1,
404         FUN = function(x) {
405           return(all(is.na(x = x)))
406         }
407       )))
408       if (length(x = rows.to.remove) > 1) {
409         combined.df <- combined.df[-rows.to.remove, ]
410       }
411       return(combined.df)
412     },
413     # Chunking functions
414     batch.scan = function(
415       chunk.size = NULL,
416       MARGIN = 2,
417       index.use = NULL,
418       dataset.use = 'matrix',
419       force.reset = FALSE
420     ) {
421       if (is.null(x = private$it) || !grepl(pattern = dataset.use, x = private$iter.dataset) || force.reset) {
422         # Check the existence of the dataset
423         private$iter.dataset <- grep(
424           pattern = dataset.use,
425           x = list.datasets(object = self),
426           value = TRUE
427         )
428         if (length(x = private$iter.dataset) != 1) {
429           private$reset_batch()
430           stop(paste0("Cannot find dataset '", dataset.use, "' in the loom file"))
431         }
432         if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
433           MARGIN <- 1
434         }
435         else if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
436           MARGIN <- 2
437         }
438         # Check the margin
439         if (!(MARGIN %in% c(1, 2))) {
440           private$reset_batch()
441           stop("MARGIN must be 1 (genes) or 2 (cells)")
442         } else {
443           private$iter.margin <- MARGIN
444         }
445         if (is.null(x = chunk.size)) {
446           chunk.size <- rev(x = self$chunksize)[private$iter.margin]
447         }
448         private$iter.chunksize <- chunk.size
449         # Set the indices to use
450         index.use <- private$iter_range(index.use = index.use)
451         # Setup our iterator
452         private$it <- ihasNext(iterable = ichunk(
453           iterable = index.use[1]:index.use[2],
454           chunkSize = chunk.size
455         ))
456         private$iter.index <- index.use
457       }
458       # Return the times we iterate, this isn't important, we only need the length of this vector
459       return(1:ceiling((private$iter.index[2] - private$iter.index[1]) / private$iter.chunksize))
460     },
461     batch.next = function(return.data = TRUE) {
462       # Ensure that we have a proper iterator
463       if (!'hasNext.ihasNext' %in% suppressWarnings(expr = methods(class = class(x = private$it)))) {
464         stop("Please setup the iterator with self$batch.scan")
465       }
466       # Do the iterating
467       if (hasNext(obj = private$it)) {
468         # Get the indices for this chunk
469         chunk.indices <- unlist(x = nextElem(obj = private$it))
470         if (return.data) {
471           # If we're working with a matrix dataset, ensure chunking on proper axis
472           if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
473             to.return <- switch(
474               EXPR = private$iter.margin,
475               '1' = self[[private$iter.dataset]][, chunk.indices], # Genes
476               '2' = self[[private$iter.dataset]][chunk.indices, ] # Cells
477             )
478           } else {
479             # Otherwise, iterating over an attribute (1 dimensional)
480             to.return <- self[[private$iter.dataset]][chunk.indices]
481           }
482         } else {
483           to.return <- chunk.indices
484         }
485         # Determine if we reset the iterator
486         if (!hasNext(obj = private$it)) {
487           private$reset_batch()
488         }
489         return(to.return)
490       } else {
491         # Just in case
492         private$reset_batch()
493         return(NULL)
494       }
495     },
496     apply = function(
497       name,
498       FUN,
499       MARGIN = 2,
500       index.use = NULL,
501       chunk.size = NULL,
502       dataset.use = 'matrix',
503       overwrite = FALSE,
504       display.progress = TRUE,
505       ...
506     ) {
507       if (self$mode == 'r') {
508         stop(private$err_mode)
509       }
510       if (!inherits(x = FUN, what = 'function')) {
511         stop("FUN must be a function")
512       }
513       # Check that we're storing our results properly
514       results.basename <- basename(path = name) # Basename of our results
515       results.dirname <- gsub(pattern = '/', replacement = '', x = dirname(path = name)) # Groupname of our results
516       dirnames <- c('row_attrs', 'col_attrs', 'layers') # Allowed group names
517       if (name %in% list.datasets(object = self)) {
518         if (overwrite) {
519           self$link_delete(name = name)
520         } else {
521           stop(paste("A dataset with the name", name, "already exists!"))
522         }
523       }
524       # Checks datset, index, and MARGIN
525       # Sets full dataset path in private$iter.dataset
526       # Sets proper MARGIN in private$iter.margin
527       batch <- self$batch.scan(
528         chunk.size = chunk.size,
529         MARGIN = MARGIN,
530         dataset.use = dataset.use,
531         force.reset = TRUE
532       )
533       MARGIN <- private$iter.margin
534       dataset.use <- private$iter.dataset
535       # Ensure that our group name is allowed
536       name.check <- which(x = dirnames == results.dirname)
537       if (!any(name.check)) {
538         private$reset_batch()
539         stop(paste(
540           "The dataset must go into one of:",
541           paste(dirnames, collapse = ', ')
542         ))
543       }
544       # Check that our group matches our iteration
545       # ie. To store in col_attrs, MARGIN must be 1
546       if (name.check %in% c(1, 2) && name.check != private$iter.margin) {
547         private$reset_batch()
548         stop(paste(
549           "Iteration must be over",
550           c('genes', 'cells')[name.check],
551           paste0("(MARGIN = ", name.check, ")"),
552           "to store results in",
553           paste0("'", dirnames[name.check], "'")
554         ))
555       }
556       # Check how we store our results
557       dataset.matrix <- ('matrix' %in% private$iter.dataset || grepl(pattern = 'layers', x = private$iter.dataset))
558       results.matrix <- name.check == 3
559       # Get a connection to the group we're iterating over
560       group <- self[[results.dirname]]
561       # Ensure index.use is integers within the bounds of [1, self$shape[MARGIN]]
562       if (!is.null(x = index.use)) {
563         # Filter index.use to values between 1 and self$shape[MARGIN]
564         index.use <- as.integer(x = index.use)
565         index.use[index.use >= 1 & index.use <= self$shape[MARGIN]]
566         index.use <- as.vector(x = na.omit(object = index.use))
567         # If we still have values, figure out NAs, otherwise set index.use to NULL
568         if (length(x = index.use) > 0) {
569           # Do a trial run to figure out the class of NAs
570           na.use <- NA
571           if (display.progress) {
572             catn("Running trial to determine class of NAs")
573           }
574           trial <- switch(
575             EXPR = MARGIN,
576             '1' = self[[dataset.use]][, 1],
577             '2' = self[[dataset.use]][1, ]
578           )
579           trial <- FUN(trial, ...)
580           if (is.list(x = trial)) {
581             trial <- unlist(x = trial)
582           }
583           class(x = na.use) <- class(x = trial)
584         } else {
585           warning("No values passed to 'index.use' fall within the data, using all values")
586           index.use <- NULL
587         }
588       }
589       if (display.progress) {
590         pb <- txtProgressBar(char = '=', style = 3)
591       }
592       # Have to initialize the dataset differently than
593       # appending to it
594       first <- TRUE
595       for (i in 1:length(x = batch)) {
596         # Get the indices we're iterating over
597         these.indices <- self$batch.next(return.data = FALSE)
598         if (is.null(x = index.use)) {
599           chunk.indices <- these.indices
600         } else {
601           chunk.indices <- index.use[index.use %in% these.indices]
602           chunk.na <- these.indices[!(these.indices %in% chunk.indices)]
603         }
604         # Get the data and apply FUN
605         chunk.data <- if (dataset.matrix) {
606           switch(
607             EXPR = MARGIN,
608             '1' = self[[dataset.use]][, chunk.indices], # Chunk genes
609             '2' = self[[dataset.use]][chunk.indices, ] # Chunk cells
610           )
611         } else {
612           self[[private$iter.datset]][chunk.indices]
613         }
614         chunk.data <- FUN(chunk.data, ...)
615         # If this is the first iteration
616         # Initialize the dataset within group, set first to FALSE
617         if (first) {
618           if (!is.null(x = index.use)) {
619             # If we had indices to chunk on, create a holding matrix the size
620             # of what the results should be, add the
621             holding <- switch(
622               EXPR = MARGIN,
623               '1' = matrix(nrow = nrow(x = chunk.data), ncol = length(x = these.indices)),
624               '2' = matrix(nrow = length(x = these.indices), ncol = ncol(x = chunk.data))
625             )
626             switch (
627               EXPR = MARGIN,
628               '1' = holding[, chunk.indices] <- chunk.data,
629               '2' = holding[chunk.indices, ] <- chunk.data
630             )
631             chunk.data <- holding
632           }
633           group[[results.basename]] <- chunk.data
634           first <- FALSE
635         } else {
636           # If we're writign to a matrix
637           # Figure out which way we're writing the data
638           if (results.matrix) {
639             switch(
640               EXPR = MARGIN,
641               '1' = group[[results.basename]][, chunk.indices] <- chunk.data,
642               '2' = group[[results.basename]][chunk.indices, ] <- chunk.data
643             )
644             if (!is.null(x = index.use)) {
645               switch(
646                 EXPR = MARGIN,
647                 '1' = group[results.basename][, chunk.na] <- na.use,
648                 '2' = group[results.basename][chunk.na, ] <- na.use,
649               )
650             }
651           } else {
652             # Just write to the vector
653             group[[results.basename]][chunk.indices] <- chunk.data
654             if (!is.null(x = index.use)) {
655               group[[results.basename]][chunk.na] <- na.use
656             }
657           }
658         }
659         if (display.progress) {
660           setTxtProgressBar(pb = pb, value = i / length(x = batch))
661         }
662       }
663       # Clean up and allow chaining
664       private$reset_batch()
665       # Load layers and attributes
666       private$load_layers()
667       private$load_attributes(MARGIN = 1) # Genes (row_attrs)
668       private$load_attributes(MARGIN = 2) # Cells (col_attrs)
669       invisible(x = self)
670     },
671     map = function(
672       FUN,
673       MARGIN = 2,
674       chunk.size = NULL,
675       index.use = NULL,
676       dataset.use = 'matrix',
677       display.progress = TRUE,
678       expected = NULL,
679       ...
680     ) {
681       if (!inherits(x = FUN, what = 'function')) {
682         stop("FUN must be a function")
683       }
684       # Checks datset, index, and MARGIN
685       # Sets full dataset path in private$iter.dataset
686       # Sets proper MARGIN in private$iter.margin
687       batch <- self$batch.scan(
688         chunk.size = chunk.size,
689         MARGIN = MARGIN,
690         index.use = index.use,
691         dataset.use = dataset.use,
692         force.reset = TRUE
693       )
694       MARGIN <- private$iter.margin
695       dataset.use <- private$iter.dataset
696       # Check how we store our results
697       # And what the shape of our dataset is
698       results.matrix <- TRUE
699       dataset.matrix <- TRUE
700       if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
701         results.matrix <- FALSE
702         dataset.matrix <- FALSE
703       } else if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
704         results.matrix <- FALSE
705         dataset.matrix <- FALSE
706       }
707       if (!is.null(x = expected)) {
708         results.matrix <- switch(
709           EXPR = expected,
710           'vector' = FALSE,
711           'matrix' = TRUE,
712           stop("'expected' must be one of 'matrix', 'vector', or NULL")
713         )
714       }
715       # Create our results holding object
716       results <- vector(mode = "list", length = length(x = batch))
717       if (display.progress) {
718         pb <- txtProgressBar(char = '=', style = 3)
719       }
720       for (i in 1:length(x = batch)) {
721         chunk.indices <- self$batch.next(return.data = FALSE)
722         chunk.data <- if (dataset.matrix) {
723           switch(
724             EXPR = MARGIN,
725             '1' = self[[dataset.use]][, chunk.indices], # Chunk genes
726             '2' = self[[dataset.use]][chunk.indices, ] # Chunk cells
727           )
728         } else {
729           self[[dataset.use]][chunk.indices]
730         }
731         results[[i]] <- FUN(chunk.data, ...)
732         if (display.progress) {
733           setTxtProgressBar(pb = pb, value = i / length(x = batch))
734         }
735       }
736       # Bring result list into vector or matrix format
737       if (class(results[[1]]) == 'numeric') {
738         results <- unlist(x = results, use.names = FALSE)
739       } else {
740         if (MARGIN == 1) {
741           results <- Reduce(f = cbind, x = results)
742         } else if (MARGIN == 2) {
743           results <- Reduce(f = rbind, x = results)
744         }
745       }
746
747       private$reset_batch()
748       return(results)
749     },
750     # Functions that modify `/matrix'
751     add.cells = function(
752       matrix.data,
753       attributes.data = NULL,
754       layers.data = NULL,
755       display.progress = TRUE
756     ) {
757       if (self$mode == 'r') {
758         stop(private$err_mode)
759       }
760       # Check inputs
761       n <- self[['matrix']]$dims[2]
762       if (display.progress) {
763         cate("Checking inputs...")
764       }
765       matrix.data <- check.matrix_data(x = matrix.data, n = n)
766       layers.data <- check.layers(
767         x = layers.data,
768         n = n,
769         layers.names = names(x = self[['layers']])
770       )
771       attributes.data <- check.col_attrs(
772         x = attributes.data,
773         attrs.names = names(x = self[['col_attrs']])
774       )
775       # Get the number of cells we're adding
776       num.cells <- c(
777         nCells.matrix_data(x = matrix.data),
778         nCells.layers(x = layers.data),
779         nCells.col_attrs(x = attributes.data)
780       )
781       num.cells <- max(num.cells)
782       # Flesh out the input data to add
783       if (display.progress) {
784         cate(paste("Adding", num.cells, "to this loom file"))
785       }
786       matrix.data <- addCells.matrix_data(x = matrix.data, n = n, m2 = num.cells)
787       layers.data <- addCells.layers(x = layers.data, n = n, m2 = num.cells)
788       attributes.data <- addCells.col_attrs(x = attributes.data, m2 = num.cells)
789       # Add the input to the loom file
790       dims.fill <- self[['matrix']]$dims[1]
791       dims.fill <- (dims.fill + 1L):(dims.fill + num.cells)
792       # Matrix data
793       if (display.progress) {
794         cate("Adding data to /matrix")
795       }
796       matrix.data <- t(x = as.matrix(x = data.frame(matrix.data)))
797       self[['matrix']][dims.fill, ] <- matrix.data
798       # Layer data
799       if (display.progress) {
800         cate("Adding data to /layers")
801         pb <- new.pb()
802         counter <- 0
803       }
804       layers.names <- names(x = self[['layers']])
805       for (i in layers.names) {
806         self[['layers']][[i]][dims.fill, ] <- t(x = layers.data[[i]])
807         if (display.progress) {
808           counter <- counter + 1
809           setTxtProgressBar(pb = pb, value = counter / length(x = layers.names))
810         }
811       }
812       # Column attributes
813       if (display.progress) {
814         cate("Adding data to /col_attrs")
815         pb <- new.pb()
816         counter <- 0
817       }
818       attrs.names <- names(x = self[['col_attrs']])
819       for (i in attrs.names) {
820         self[['col_attrs']][[i]][dims.fill] <- attributes.data[[i]]
821         if (display.progress) {
822           counter <- counter + 1
823           setTxtProgressBar(pb = pb, value = counter / length(x = attrs.names))
824         }
825       }
826       # # Load layers and attributes
827       # private$load_layers()
828       # private$load_attributes(MARGIN = 1)
829       # private$load_attributes(MARGIN = 2)
830       # self$shape <- self[['matrix']]$dims
831     },
832     add.loom = function(
833       other,
834       other.key = NULL,
835       self.key = NULL,
836       ...
837     ) {
838       if (self$mode == 'r') {
839         stop(private$err_mode)
840       }
841       # Connect to the other loom file
842       if (inherits(x = other, what = 'loom')) {
843         ofile <- other
844       } else if (is.character(x = other)) {
845         ofile <- connect(filename = other)
846       } else {
847         stop("'other' must be either a loom object or a path to a loom file")
848       }
849       # If we have row keys to use
850       if (!is.null(x = other.key) && !is.null(x = self.key)) {
851         other.key <- basename(path = other.key)
852         self.key <- basename(path = self.key)
853         tryCatch(
854           expr = other.key <- other[['row_attrs']][[other.key]][],
855           error = function(e) {
856             if (is.character(x = other)) {
857               ofile$close
858             }
859             stop("Failed to find the gene names dataset in the other loom file")
860           }
861         )
862         tryCatch(
863           expr = self.key <- self[['row_attrs']][[self.key]][],
864           error = function(e) {
865             if (is.character(x = other)) {
866               ofile$close
867             }
868             stop("Failed to find the gene names dataset in this loom file")
869           }
870         )
871         # Match rows
872         rows.use <- match(x = other.key, table = self.key, nomatch = 0)
873         rows.use <- rows.use[rows.use > 0]
874       } else {
875         cate("Adding the loom file as-is, assuming in the correct order")
876         Sys.sleep(time = 3)
877         rows.use <- 1:other[['matrix']]$dims[2]
878       }
879       if (max(rows.use) > self[['matrix']]$dims[2]) {
880         stop("More genes in the other loom file than present in this one")
881       } else if (max(rows.use) < self[['matrix']]$dims[2]) {
882         cate("...")
883       }
884       # Clean up
885       if (is.character(x = other)) {
886         ofile$close_all()
887       }
888     }
889   ),
890   # Private fields and methods
891   # @field err_init An error message for if this object hasn't been created with loomR::create or loomR::connect
892   # @field err_mode An error message for if this object is in read-only mode
893   # @field it Iterator object for batch.scan and batch.next
894   # @field iter.dataset Dataset for iterating on
895   # @field iter.margin Margin to iterate over
896   # @field iter.index Index values for iteration
897   # @field skipped.validation Was validation skipped?
898   # \describe{
899   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
900   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
901   #   \item{\code{reset_batch()}}{Reset the batch iterator fields}
902   #   \item{\code{iter_range(index.use)}}{Get the range of indices for a batch iteration}
903   # }
904   private = list(
905     # Fields
906     err_init = "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",
907     err_mode = "Cannot modify a loom file in read-only mode",
908     it = NULL,
909     iter.chunksize = NULL,
910     iter.dataset = NULL,
911     iter.margin = NULL,
912     iter.index = NULL,
913     skipped.validation = FALSE,
914     # Methods
915     load_attributes = function(MARGIN) {
916       attribute <- switch(
917         EXPR = MARGIN,
918         '1' = 'row_attrs',
919         '2' = 'col_attrs',
920         stop('Invalid attribute dimension')
921       )
922       group <- self[[attribute]]
923       attributes <- unlist(x = lapply(
924         X = names(x = group),
925         FUN = function(x) {
926           d <- list(group[[x]])
927           names(x = d) <- x
928           return(d)
929         }
930       ))
931       if (MARGIN == 1) {
932         self$row.attrs <- attributes
933       } else if (MARGIN == 2) {
934         self$col.attrs <- attributes
935       }
936     },
937     load_layers = function() {
938       if (private$skipped.validation) {
939         if (getOption(x = 'verbose')) {
940           warning("Not loading layers field")
941         }
942       } else {
943         self$layers <- unlist(x = lapply(
944           X = names(x = self[['layers']]),
945           FUN = function(n) {
946             d <- list(self[[paste('layers', n, sep = '/')]])
947             names(x = d) <- n
948             return(d)
949           }
950         ))
951       }
952     },
953     reset_batch = function() {
954       private$it <- NULL
955       private$iter.chunksize <- NULL
956       private$iter.dataset <- NULL
957       private$iter.margin <- NULL
958       private$iter.index <- NULL
959     },
960     iter_range = function(index.use) {
961       if (is.null(private$iter.margin)) {
962         stop("Batch processing hasn't been set up")
963       }
964       if (is.null(x = index.use)) {
965         # If no index was provided, use entire range for this margin
966         index.use <- c(1, self$shape[private$iter.margin])
967       } else if (length(x = index.use) == 1) {
968         # If one index was provided, start at one and go to index
969         index.use <- c(1, index.use)
970       } else {
971         # Use index.use[1] and index.use[2]
972         index.use <- c(index.use[1], index.use[2])
973       }
974       # Ensure the indices provided fit within the range of the dataset
975       index.use[1] <- max(1, index.use[1])
976       index.use[2] <- min(index.use[2], self$shape[private$iter.margin])
977       # Ensure that index.use[1] is greater than index.use[2]
978       if (index.use[1] > index.use[2]) {
979         stop(paste0(
980           "Starting index (",
981           index.use[1],
982           ") must be lower than the ending index (",
983           index.use[2],
984           ")"
985         ))
986       }
987       return(index.use)
988     }
989   )
990 )
991
992 #' Create a loom object
993 #'
994 #' @param filename The name of the new loom file
995 #' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
996 #' @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}
997 #' @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}
998 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
999 #' @param chunk.size How many rows of \code{data} should we stream to the loom file at any given time?
1000 #' @param overwrite Overwrite an already existing loom file?
1001 #'
1002 #' @return A connection to a loom file
1003 #'
1004 #' @importFrom utils packageVersion txtProgressBar setTxtProgressBar
1005 #'
1006 #' @seealso \code{\link{loom-class}}
1007 #'
1008 #' @export
1009 #'
1010 create <- function(
1011   filename,
1012   data,
1013   gene.attrs = NULL,
1014   cell.attrs = NULL,
1015   layers = NULL,
1016   chunk.dims = 'auto',
1017   chunk.size = 1000,
1018   overwrite = FALSE,
1019   display.progress = TRUE
1020 ) {
1021   mode <- ifelse(test = overwrite, yes = 'w', no = 'w-')
1022   if (file.exists(filename) && !overwrite) {
1023     stop(paste('File', filename, 'already exists!'))
1024   }
1025   if (!inherits(x = data, what = c('matrix', 'Matrix'))) {
1026     data <- as.matrix(x = data)
1027   }
1028   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
1029     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
1030   } else if (length(x = chunk.dims) == 1) {
1031     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
1032       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
1033     }
1034   } else {
1035     chunk.dims <- as.integer(x = chunk.dims)
1036   }
1037   new.loom <- loom$new(filename = filename, mode = mode)
1038   # Create the matrix
1039   # new.loom$create_dataset(
1040   #   name = 'matrix',
1041   #   robj = data,
1042   #   chunk_dims = chunk.dims
1043   # )
1044   dtype <- getDtype(x = data[1, 1])
1045   new.loom$create_dataset(
1046     name = 'matrix',
1047     dtype = dtype,
1048     dims = dim(x = data)
1049   )
1050   chunk.points <- chunkPoints(
1051     data.size = dim(x = data)[1],
1052     chunk.size = chunk.size
1053   )
1054   if (display.progress) {
1055     pb <- txtProgressBar(char = '=', style = 3)
1056   }
1057   for (col in 1:ncol(x = chunk.points)) {
1058     row.start <- chunk.points[1, col]
1059     row.end <- chunk.points[2, col]
1060     new.loom[['matrix']][row.start:row.end, ] <- as.matrix(x = data[row.start:row.end, ])
1061     if (display.progress) {
1062       setTxtProgressBar(pb = pb, value = col / ncol(x = chunk.points))
1063     }
1064   }
1065   new.loom$matrix <- new.loom[['matrix']]
1066   new.loom$shape <- rev(x = new.loom[['matrix']]$dims)
1067   # Groups
1068   new.loom$create_group(name = 'layers')
1069   new.loom$create_group(name = 'row_attrs')
1070   new.loom$create_group(name = 'col_attrs')
1071   # Check for the existance of gene or cell names
1072   if (!is.null(x = colnames(x = data))) {
1073     new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
1074   }
1075   if (!is.null(x = rownames(x = data))) {
1076     new.loom$add.col.attribute(attribute = list('cell_names' = rownames(x = data)))
1077   }
1078   # Store some constants as HDF5 attributes
1079   h5attr(x = new.loom, which = 'version') <- new.loom$version
1080   h5attr(x = new.loom, which = 'chunks') <- paste0(
1081     '(',
1082     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
1083     ')'
1084   )
1085   # Add layers
1086   if (!is.null(x = layers)) {
1087     new.loom$add.layer(layer = layers)
1088   }
1089   if (!is.null(x = gene.attrs)) {
1090     new.loom$add.row.attribute(attribute = gene.attrs)
1091   }
1092   if (!is.null(x = cell.attrs)) {
1093     new.loom$add.col.attribute(attribute = cell.attrs)
1094   }
1095   # Set last bit of information
1096   chunks <- new.loom[['matrix']]$chunk_dims
1097   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
1098   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
1099   chunks <- unlist(x = strsplit(x = chunks, split = ','))
1100   new.loom$chunksize <- as.integer(x = chunks)
1101   # Return the connection
1102   return(new.loom)
1103 }
1104
1105 #' Connect to a loom file
1106 #'
1107 #' @param filename The loom file to connect to
1108 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write
1109 #' @param skip.validate Skip the validation steps, use only for extremely large loom files
1110 #'
1111 #' @return A loom file connection
1112 #'
1113 #' @seealso \code{\link{loom-class}}
1114 #'
1115 #' @export
1116 #'
1117 connect <- function(filename, mode = "r", skip.validate = FALSE) {
1118   if (!(mode %in% c('r', 'r+'))) {
1119     stop("'mode' must be one of 'r' or 'r+'")
1120   }
1121   new.loom <- loom$new(filename = filename, mode = mode, skip.validate = skip.validate)
1122   return(new.loom)
1123 }
1124
1125 #' Subset a loom file
1126 #'
1127 #' @param x A loom object
1128 #' @param m Rows (genes) to subset, defaults to all rows
1129 #' @param n Columns (cells) to subset, defaults to all columns
1130 #' @param filename Filename for new loom object, defaults to ...
1131 #' @param overwrite Overwrite \code{filename} if already exists?
1132 #' @param display.progress Display progress as we're copying over data
1133 #'
1134 #' @return A loom object connected to \code{filename}
1135 #'
1136 #' @importFrom utils setTxtProgressBar
1137 #'
1138 #' @export subset.loom
1139 #' @method subset loom
1140 #'
1141 subset.loom <- function(
1142   x,
1143   m = NULL,
1144   n = NULL,
1145   filename = NULL,
1146   overwrite = FALSE,
1147   display.progress = TRUE,
1148   ...
1149 ) {
1150   # Set some defaults
1151   if (is.null(x = m)) {
1152     m <- 1:x$shape[1]
1153   }
1154   if (is.null(x = n)) {
1155     n <- 1:x$shape[2]
1156   }
1157   if (is.null(x = filename)) {
1158     filename <- paste(
1159       unlist(x = strsplit(x = x$filename, split = '.', fixed = TRUE)),
1160       collapse = '_subset.'
1161     )
1162   }
1163   if (length(x = m) == 1) {
1164     m <- 1:m
1165   }
1166   if (length(x = n) == 1) {
1167     n <- 1:n
1168   }
1169   # Ensure that m and n are within the bounds of the loom file
1170   if (max(m) > x$shape[1] || max(n) > x$shape[2]) {
1171     stop(paste(
1172       "'m' and 'n' must be less than",
1173       x$shape[1],
1174       "and",
1175       x$shape[2],
1176       "respectively"
1177     ))
1178   }
1179   extension <- rev(x = unlist(x = strsplit(x = filename, split = '.', fixed = TRUE)))[1]
1180   # Ensure that we call our new file a loom file
1181   if (extension != 'loom') {
1182     filename <- paste0(filename, '.loom')
1183   }
1184   if (display.progress) {
1185     catn("Writing new loom file to", filename)
1186   }
1187   # Make the loom file
1188   new.loom <- create(
1189     filename = filename,
1190     # data = t(x = x[['matrix']][n, m]),
1191     data = x[['matrix']][n, m],
1192     overwrite = overwrite
1193   )
1194   # Add row attributes
1195   row.attrs <- list.datasets(object = x, path = 'row_attrs', full.names = TRUE)
1196   if (length(x = row.attrs) > 0) {
1197     if (display.progress) {
1198       catn("\nAdding", length(x = row.attrs), "row attributes")
1199       pb <- new.pb()
1200       counter <- 0
1201     }
1202     for (row.attr in row.attrs) {
1203       row.list <- list(x[[row.attr]][m])
1204       names(x = row.list) <- basename(path = row.attr)
1205       new.loom$add.row.attribute(attribute = row.list)
1206       if (display.progress) {
1207         counter <- counter + 1
1208         setTxtProgressBar(pb = pb, value = counter / length(x = row.attrs))
1209       }
1210     }
1211   } else {
1212     warning("No row attributes found")
1213   }
1214   # Add col attributes
1215   col.attrs <- list.datasets(object = x, path = 'col_attrs', full.names = TRUE)
1216   if (length(x = col.attrs) > 0) {
1217     if (display.progress) {
1218       catn("\nAdding", length(x = col.attrs), "column attributes")
1219       pb <- new.pb()
1220       counter <- 0
1221     }
1222     for (col.attr in col.attrs) {
1223       col.list <- list(x[[col.attr]][n])
1224       names(x = col.list) <- basename(path = col.attr)
1225       new.loom$add.col.attribute(attribute = col.list)
1226       if (display.progress) {
1227         counter <- counter + 1
1228         setTxtProgressBar(pb = pb, value = counter / length(x = col.attrs))
1229       }
1230     }
1231   } else {
1232     warning("No column attributes found")
1233   }
1234   # Add layers
1235   layers <- list.datasets(object = x, path = 'layers', full.names = TRUE)
1236   if (length(x = layers) > 0) {
1237     if (display.progress) {
1238       catn("\nAdding", length(x = layers), "layers")
1239       pb <- new.pb()
1240       counter <- 0
1241     }
1242     for (layer in layers) {
1243       layer.list <- list(x[[layer]][n, m])
1244       names(x = layer.list) <- basename(path = layer)
1245       new.loom$add.layer(layers = layer.list)
1246       if (display.progress) {
1247         counter <- counter + 1
1248         setTxtProgressBar(pb = pb, value = counter / length(x = layers))
1249       }
1250     }
1251   } else {
1252     warning("No layers found")
1253   }
1254   new.loom$flush()
1255   return(new.loom)
1256 }
1257
1258 # #need to comment
1259 # #need to add progress bar
1260 # #but otherwise, pretty cool
1261 # #for paul to try :
1262 # # f <- connect("~/Downloads/10X43_1.loom")
1263 # # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
1264 # # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
1265 # map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
1266 #   n_func = length(f_list)
1267 #   if (n_func == 1) {
1268 #     f_list <- list(f_list)
1269 #   }
1270 #   if (MARGIN == 1) {
1271 #     results <- list()
1272 #     for (j in 1:n_func) {
1273 #       results[[j]] <- numeric(0)
1274 #     }
1275 #     rows_per_chunk <- chunksize
1276 #     ix <- 1
1277 #     while (ix <= self@shape[1]) {
1278 #       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
1279 #       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
1280 #       for (j in 1:n_func) {
1281 #         new_results <- apply(chunk, 1, FUN = f_list[[j]])
1282 #         results[[j]] <- c(results[[j]], new_results)
1283 #       }
1284 #       ix <- ix + chunksize
1285 #     }
1286 #   }
1287 #   if (MARGIN == 2) {
1288 #     results <- list()
1289 #     for (j in 1:n_func) {
1290 #       results[[j]] <- numeric(0)
1291 #     }
1292 #     cols_per_chunk <- chunksize
1293 #     ix <- 1
1294 #     while (ix <= self@shape[2]) {
1295 #       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
1296 #       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
1297 #       for (j in 1:n_func) {
1298 #         new_results <- apply(chunk, 2, FUN = f_list[[j]])
1299 #         results[[j]] <- c(results[[j]], new_results)
1300 #       }
1301 #       ix <- ix + chunksize
1302 #     }
1303 #   }
1304 #   if (n_func == 1) {
1305 #     results <- results[[1]]
1306 #   }
1307 #   return(results)
1308 # }