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