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