add ability to extract a dataframe of multiple row_attr or col_attr elements
[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(attribute.layer = c("row", "col"),
285                            attribute.names = NULL,
286                            row.names = "gene_names",
287                            col.names = "cell_names"){
288       # takes multiple row_attrs or col_attrs and combines them into a data.frame,
289       # removing rows or columns that are entirely NAs.
290       #
291       # attribute.layer either "row" for row_attrs or "col" col_attrs
292       # attribute.names name of rows or columns to combine into matrix
293       # row.names either a character vector or name of an element in row_attrs
294       # col.names either a character vector or name of an element in col_attrs
295
296       if(!attribute.layer %in% c("row", "col")) {
297         stop("Invalid attribute.layer. Please select either 'row' or 'col'.")
298       }
299       attribute.layer <- paste0(attribute.layer, "_attrs")
300       # check that attribute names are present
301       if(!all(attribute.names %in% self[[attribute.layer]]$names)) {
302         invalid.names <- attribute.names[which(!attribute.names %in% self[[attribute.layer]]$names)]
303         stop(paste0("Invalid attribute.names: ", paste0(invalid.names, collapse = ", ")))
304       }
305       if(attribute.layer == "row_attrs") {
306         combined.df <- data.frame(self[[paste0(attribute.layer, "/", attribute.names[1])]][],
307                                   row.names = self[[paste0(attribute.layer, "/", row.names)]][])
308       } else {
309         combined.df <- data.frame(self[[paste0(attribute.layer, "/", attribute.names[1])]][],
310                                   row.names = self[[paste0(attribute.layer, "/", col.names)]][])
311       }
312       if (length(attribute.names) > 1) {
313         for (i in 2:length(attribute.names)) {
314           combined.df[, attribute.names[i]] <- self[[paste0(attribute.layer, "/", attribute.names[i])]][]
315         }
316       }
317       colnames(combined.df) <- attribute.names
318       # check if any row is all NAs
319       rows.to.remove <- unname(which(apply(X = combined.df, MARGIN = 1, FUN = function(x) all(is.na(x)))))
320       if (length(rows.to.remove) > 1) {
321         combined.df <- combined.df[-rows.to.remove, ]
322       }
323       return(combined.df)
324     },
325     # Chunking functions
326     batch.scan = function(
327       chunk.size = NULL,
328       MARGIN = 2,
329       index.use = NULL,
330       dataset.use = 'matrix',
331       force.reset = FALSE
332     ) {
333       if (is.null(x = private$it) || !grepl(pattern = dataset.use, x = private$iter.dataset) || force.reset) {
334         # Check the existence of the dataset
335         private$iter.dataset <- grep(
336           pattern = dataset.use,
337           x = list.datasets(object = self),
338           value = TRUE
339         )
340         if (length(x = private$iter.dataset) != 1) {
341           private$reset_batch()
342           stop(paste0("Cannot find dataset '", dataset.use, "' in the loom file"))
343         }
344         if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
345           MARGIN <- 1
346         }
347         else if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
348           MARGIN <- 2
349         }
350         # Check the margin
351         if (!(MARGIN %in% c(1, 2))) {
352           private$reset_batch()
353           stop("MARGIN must be 1 (genes) or 2 (cells)")
354         } else {
355           private$iter.margin <- MARGIN
356         }
357         if (is.null(x = chunk.size)) {
358           chunk.size <- rev(x = self$chunksize)[private$iter.margin]
359         }
360         private$iter.chunksize <- chunk.size
361         # Set the indices to use
362         index.use <- private$iter_range(index.use = index.use)
363         # Setup our iterator
364         private$it <- ihasNext(iterable = ichunk(
365           iterable = index.use[1]:index.use[2],
366           chunkSize = chunk.size
367         ))
368         private$iter.index <- index.use
369       }
370       # Return the times we iterate, this isn't important, we only need the length of this vector
371       return(1:ceiling((private$iter.index[2] - private$iter.index[1]) / private$iter.chunksize))
372     },
373     batch.next = function(return.data = TRUE) {
374       # Ensure that we have a proper iterator
375       if (!'hasNext.ihasNext' %in% suppressWarnings(expr = methods(class = class(x = private$it)))) {
376         stop("Please setup the iterator with self$batch.scan")
377       }
378       # Do the iterating
379       if (hasNext(obj = private$it)) {
380         # Get the indices for this chunk
381         chunk.indices <- unlist(x = nextElem(obj = private$it))
382         if (return.data) {
383           # If we're working with a matrix dataset, ensure chunking on proper axis
384           if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
385             to.return <- switch(
386               EXPR = private$iter.margin,
387               '1' = self[[private$iter.dataset]][, chunk.indices], # Genes
388               '2' = self[[private$iter.dataset]][chunk.indices, ] # Cells
389             )
390           } else {
391             # Otherwise, iterating over an attribute (1 dimensional)
392             to.return <- self[[private$iter.dataset]][chunk.indices]
393           }
394         } else {
395           to.return <- chunk.indices
396         }
397         # Determine if we reset the iterator
398         if (!hasNext(obj = private$it)) {
399           private$reset_batch()
400         }
401         return(to.return)
402       } else {
403         # Just in case
404         private$reset_batch()
405         return(NULL)
406       }
407     },
408     apply = function(
409       name,
410       FUN,
411       MARGIN = 2,
412       chunk.size = NULL,
413       dataset.use = 'matrix',
414       overwrite = FALSE,
415       display.progress = TRUE,
416       ...
417     ) {
418       if (self$mode == 'r') {
419         stop("Cannot write to disk in read-only mode")
420       }
421       if (!inherits(x = FUN, what = 'function')) {
422         stop("FUN must be a function")
423       }
424       # Check that we're storing our results properly
425       results.basename <- basename(path = name) # Basename of our results
426       results.dirname <- gsub(pattern = '/', replacement = '', x = dirname(path = name)) # Groupname of our results
427       dirnames <- c('row_attrs', 'col_attrs', 'layers') # Allowed group names
428       if (name %in% list.datasets(object = self)) {
429         if (overwrite) {
430           self$link_delete(name = name)
431         } else {
432           stop(paste("A dataset with the name", name, "already exists!"))
433         }
434       }
435       # Checks datset, index, and MARGIN
436       # Sets full dataset path in private$iter.dataset
437       # Sets proper MARGIN in private$iter.margin
438       batch <- self$batch.scan(
439         chunk.size = chunk.size,
440         MARGIN = MARGIN,
441         dataset.use = dataset.use,
442         force.reset = TRUE
443       )
444       MARGIN <- private$iter.margin
445       dataset.use <- private$iter.dataset
446       # Ensure that our group name is allowed
447       name.check <- which(x = dirnames == results.dirname)
448       if (!any(name.check)) {
449         private$reset_batch()
450         stop(paste(
451           "The dataset must go into one of:",
452           paste(dirnames, collapse = ', ')
453         ))
454       }
455       # Check that our group matches our iteration
456       # ie. To store in col_attrs, MARGIN must be 1
457       if (name.check %in% c(1, 2) && name.check != private$iter.margin) {
458         private$reset_batch()
459         stop(paste(
460           "Iteration must be over",
461           c('genes', 'cells')[name.check],
462           paste0("(MARGIN = ", name.check, ")"),
463           "to store results in",
464           paste0("'", dirnames[name.check], "'")
465         ))
466       }
467       # Check how we store our results
468       dataset.matrix <- ('matrix' %in% private$iter.dataset || grepl(pattern = 'layers', x = private$iter.dataset))
469       results.matrix <- name.check == 3
470       # Get a connection to the group we're iterating over
471       group <- self[[results.dirname]]
472       if (display.progress) {
473         pb <- txtProgressBar(char = '=', style = 3)
474       }
475       # Have to initialize the dataset differently than
476       first <- TRUE
477       for (i in 1:length(x = batch)) {
478         # Get the indices we're iterating over
479         chunk.indices <- self$batch.next(return.data = FALSE)
480         # Get the data and apply FUN
481         chunk.data <- if (dataset.matrix) {
482           switch(
483             EXPR = MARGIN,
484             '1' = self[[dataset.use]][, chunk.indices], # Chunk genes
485             '2' = self[[dataset.use]][chunk.indices, ] # Chunk cells
486           )
487         } else {
488           self[[private$iter.datset]][chunk.indices]
489         }
490         chunk.data <- FUN(chunk.data, ...)
491         # If this is the first iteration
492         # Initialize the dataset within group, set first to FALSE
493         if (first) {
494           group[[results.basename]] <- chunk.data
495           first <- FALSE
496         } else {
497           # If we're writign to a matrix
498           # Figure out which way we're writing the data
499           if (results.matrix) {
500             switch(
501               EXPR = MARGIN,
502               '1' = group[[results.basename]][, chunk.indices] <- chunk.data,
503               '2' = group[[results.basename]][chunk.indices, ] <- chunk.data
504             )
505           } else {
506             # Just write to the vector
507             group[[results.basename]][chunk.indices] <- chunk.data
508           }
509         }
510         if (display.progress) {
511           setTxtProgressBar(pb = pb, value = i / length(x = batch))
512         }
513       }
514       # Clean up and allow chaining
515       private$reset_batch()
516       # Load layers and attributes
517       private$load_layers()
518       private$load_attributes(MARGIN = 1) # Genes (row_attrs)
519       private$load_attributes(MARGIN = 2) # Cells (col_attrs)
520       invisible(x = self)
521     },
522     map = function(
523       FUN,
524       MARGIN = 2,
525       chunk.size = NULL,
526       index.use = NULL,
527       dataset.use = 'matrix',
528       display.progress = TRUE,
529       expected = NULL,
530       ...
531     ) {
532       if (!inherits(x = FUN, what = 'function')) {
533         stop("FUN must be a function")
534       }
535       # Checks datset, index, and MARGIN
536       # Sets full dataset path in private$iter.dataset
537       # Sets proper MARGIN in private$iter.margin
538       batch <- self$batch.scan(
539         chunk.size = chunk.size,
540         MARGIN = MARGIN,
541         index.use = index.use,
542         dataset.use = dataset.use,
543         force.reset = TRUE
544       )
545       MARGIN <- private$iter.margin
546       dataset.use <- private$iter.dataset
547       # Check how we store our results
548       # And what the shape of our dataset is
549       results.matrix <- TRUE
550       dataset.matrix <- TRUE
551       if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
552         results.matrix <- FALSE
553         dataset.matrix <- FALSE
554       } else if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
555         results.matrix <- FALSE
556         dataset.matrix <- FALSE
557       }
558       if (!is.null(x = expected)) {
559         results.matrix <- switch(
560           EXPR = expected,
561           'vector' = FALSE,
562           'matrix' = TRUE,
563           stop("'expected' must be one of 'matrix', 'vector', or NULL")
564         )
565       }
566       # Determine the shape of our results
567       index.use <- private$iter_range(index.use = index.use)
568       # Create our results holding object
569       if (results.matrix) {
570         switch(
571           EXPR = private$iter.margin,
572           '1' = results <- matrix( # Genes, make matrix with nCells rows and range(index.use) columns
573             nrow = self$shape[1],
574             ncol = length(x = index.use[1]:index.use[2])
575           ),
576           '2' = results <- matrix( # Cells, make matrix with range(index.use) rows and nGenes columns
577             nrow = length(x = index.use[1]:index.use[2]),
578             ncol = self$shape[2]
579           )
580         )
581       } else {
582         results <- vector(length = length(x = index.use[1]:index.use[2]))
583       }
584       if (display.progress) {
585         pb <- txtProgressBar(char = '=', style = 3)
586       }
587       for (i in 1:length(x = batch)) {
588         chunk.indices <- self$batch.next(return.data = FALSE)
589         chunk.data <- if (dataset.matrix) {
590           switch(
591             EXPR = MARGIN,
592             '1' = self[[dataset.use]][, chunk.indices], # Chunk genes
593             '2' = self[[dataset.use]][chunk.indices, ] # Chunk cells
594           )
595         } else {
596           self[[dataset.use]][chunk.indices]
597         }
598         chunk.data <- FUN(chunk.data, ...)
599         if (results.matrix) {
600           if (MARGIN == 1) {
601             results[, chunk.indices] <- chunk.data
602           } else if (MARGIN == 2) {
603             results[chunk.indices, ] <- chunk.data
604           }
605         } else {
606           results[chunk.indices] <- chunk.data
607         }
608         if (display.progress) {
609           setTxtProgressBar(pb = pb, value = i / length(x = batch))
610         }
611       }
612       private$reset_batch()
613       return(results)
614     },
615     # Functions that modify `/matrix'
616     add.cells = function(matrix.data, attributes.data = NULL, layers.data = NULL) {
617       # matrix.data is a vector of data for one cell or a list of data for several cells
618       # each entry in matrix.data must be the same length as number of genes
619       # attributes.data is an optional list or vector (with optional names) for col_attrs
620       # each entry in col_attrs must be the same length as the number of cells being added (NAs added for those that aren't)
621       # layers.data is an optional list (with optional names) for layers
622       # 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
623       if (is.vector(x = matrix.data) && !is.list(x = matrix.data)) {
624         matrix.data <- list(matrix.data)
625       }
626       list.check <- vapply(
627         X = list(matrix.data, attributes.data, layers.data),
628         FUN = is.list,
629         FUN.VALUE = logical(length = 1L)
630       )
631       if (!all(list.check)) {
632         stop("'matrix.data', 'attributes.data', and 'layers.data' must be lists")
633       }
634       lengths <- vector(
635         mode = 'integer',
636         length = 1 + length(x = attributes.data) + length(x = layers.data)
637       )
638       lengths[1] <- length(x = matrix.data)
639       attributes.end <- 1 + length(x = attributes.data)
640       if (attributes.end != 1) {
641         lengths[2:attributes.end] <- vapply(
642           X = attributes.data,
643           FUN = length,
644           FUN.VALUE = integer(length = 1L)
645         )
646       }
647       if (attributes.end != length(x = lengths)) {
648         lengths[(attributes.end + 1):length(x = lengths)] <- vapply(
649           X = layers.data,
650           FUN = length,
651           FUN.VALUE = integer(length = 1L)
652         )
653       }
654       num.cells.added <- max(lengths)
655       return(num.cells.added)
656       # private$load_layers()
657       # private$load_attributes(MARGIN = 1)
658       # private$load_attributes(MARGIN = 2)
659       return(lengths)
660     }
661     # add.loom = function() {}
662   ),
663   # Private fields and methods
664   # @field err_msg A simple error message if this object hasn't been created with loomR::create or loomR::connect
665   # @field it Iterator object for batch.scan and batch.next
666   # @field iter.dataset Dataset for iterating on
667   # @field iter.margin Margin to iterate over
668   # @field iter.index # Index values for iteration
669   # \describe{
670   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
671   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
672   #   \item{\code{reset_batch()}}{Reset the batch iterator fields}
673   #   \item{\code{iter_range(index.use)}}{Get the range of indices for a batch iteration}
674   # }
675   private = list(
676     # Fields
677     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",
678     it = NULL,
679     iter.chunksize = NULL,
680     iter.dataset = NULL,
681     iter.margin = NULL,
682     iter.index = NULL,
683     skipped.validation = FALSE,
684     # Methods
685     load_attributes = function(MARGIN) {
686       attribute <- switch(
687         EXPR = MARGIN,
688         '1' = 'row_attrs',
689         '2' = 'col_attrs',
690         stop('Invalid attribute dimension')
691       )
692       group <- self[[attribute]]
693       attributes <- unlist(x = lapply(
694         X = names(x = group),
695         FUN = function(x) {
696           d <- list(group[[x]])
697           names(x = d) <- x
698           return(d)
699         }
700       ))
701       if (MARGIN == 1) {
702         self$row.attrs <- attributes
703       } else if (MARGIN == 2) {
704         self$col.attrs <- attributes
705       }
706     },
707     load_layers = function() {
708       if (private$skipped.validation) {
709         if (getOption(x = 'verbose')) {
710           warning("Not loading layers field")
711         }
712       } else {
713         self$layers <- unlist(x = lapply(
714           X = names(x = self[['layers']]),
715           FUN = function(n) {
716             d <- list(self[[paste('layers', n, sep = '/')]])
717             names(x = d) <- n
718             return(d)
719           }
720         ))
721       }
722     },
723     reset_batch = function() {
724       private$it <- NULL
725       private$iter.chunksize <- NULL
726       private$iter.dataset <- NULL
727       private$iter.margin <- NULL
728       private$iter.index <- NULL
729     },
730     iter_range = function(index.use) {
731       if (is.null(private$iter.margin)) {
732         stop("Batch processing hasn't been set up")
733       }
734       if (is.null(x = index.use)) {
735         # If no index was provided, use entire range for this margin
736         index.use <- c(1, rev(x = self$shape)[private$iter.margin])
737       } else if (length(x = index.use) == 1) {
738         # If one index was provided, start at one and go to index
739         index.use <- c(1, index.use)
740       } else {
741         # Use index.use[1] and index.use[2]
742         index.use <- c(index.use[1], index.use[2])
743       }
744       # Ensure the indices provided fit within the range of the dataset
745       index.use[1] <- max(1, index.use[1])
746       index.use[2] <- min(index.use[2], rev(x = self$shape)[private$iter.margin])
747       # Ensure that index.use[1] is greater than index.use[2]
748       if (index.use[1] > index.use[2]) {
749         stop(paste0(
750           "Starting index (",
751           index.use[1],
752           ") must be lower than the ending index (",
753           index.use[2],
754           ")"
755         ))
756       }
757       return(index.use)
758     }
759   )
760 )
761
762 #' Create a loom object
763 #'
764 #' @param filename The name of the new loom file
765 #' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
766 #' @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}
767 #' @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}
768 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
769 #'
770 #' @return A connection to a loom file
771 #'
772 #' @importFrom utils packageVersion
773 #'
774 #' @seealso \code{\link{loom-class}}
775 #'
776 #' @export
777 #'
778 create <- function(
779   filename,
780   data,
781   gene.attrs = NULL,
782   cell.attrs = NULL,
783   layers = NULL,
784   chunk.dims = 'auto'
785 ) {
786   if (file.exists(filename)) {
787     stop(paste('File', file, 'already exists!'))
788   }
789   if (!is.matrix(x = data)) {
790     data <- as.matrix(x = data)
791   }
792   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
793     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
794   } else if (length(x = chunk.dims == 1)) {
795     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
796       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
797     }
798   } else {
799     chunk.dims <- as.integer(x = chunk.dims)
800   }
801   new.loom <- loom$new(filename = filename, mode = 'w-')
802   # Create the matrix
803   new.loom$create_dataset(
804     name = 'matrix',
805     robj = data,
806     chunk_dims = chunk.dims
807   )
808   new.loom$matrix <- new.loom[['matrix']]
809   new.loom$shape <- new.loom[['matrix']]$dims
810   # Groups
811   new.loom$create_group(name = 'layers')
812   new.loom$create_group(name = 'row_attrs')
813   new.loom$create_group(name = 'col_attrs')
814   # Check for the existance of gene or cell names
815   if (!is.null(x = colnames(x = data))) {
816     new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
817   }
818   if (!is.null(x = rownames(x = data))) {
819     new.loom$add.col.attribute(attribute = list('cell_names' = rownames(x = data)))
820   }
821   # Store some constants as HDF5 attributes
822   h5attr(x = new.loom, which = 'version') <- new.loom$version
823   h5attr(x = new.loom, which = 'chunks') <- paste0(
824     '(',
825     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
826     ')'
827   )
828   # Add layers
829   if (!is.null(x = layers)) {
830     new.loom$add.layer(layer = layers)
831   }
832   if (!is.null(x = gene.attrs)) {
833     new.loom$add.row.attribute(attribute = gene.attrs)
834   }
835   if (!is.null(x = cell.attrs)) {
836     new.loom$add.col.attribute(attribute = cell.attrs)
837   }
838   # Set last bit of information
839   chunks <- new.loom[['matrix']]$chunk_dims
840   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
841   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
842   chunks <- unlist(x = strsplit(x = chunks, split = ','))
843   new.loom$chunksize <- as.integer(x = chunks)
844   # Return the connection
845   return(new.loom)
846 }
847
848 # Validate a loom object
849 #
850 # @param object A loom object
851 #
852 # @return None, errors out if object is an invalid loom connection
853 #
854 # @seealso \code{\link{loom-class}}
855 #
856 #
857 validateLoom <- function(object) {
858   # A loom file is a specific HDF5
859   # We need a dataset in /matrix that's a two-dimensional dense matrix
860   root.datasets <- list.datasets(object = object, path = '/', recursive = FALSE)
861   if (length(x = root.datasets) != 1) {
862     stop("There can only be one dataset at the root of the loom file")
863   }
864   if (root.datasets != 'matrix') {
865     stop("The root dataset must be called '/matrix'")
866   }
867   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
868   required.groups <- c('row_attrs', 'col_attrs', 'layers')
869   dim.matrix <- object[[root.datasets]]$dims # Columns x Rows
870   names(x = dim.matrix) <- required.groups[c(2, 1)]
871   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
872   group.msg <- paste0(
873     "There can only be three groups in the loom file: '",
874     paste(required.groups, collapse = "', '"),
875     "'"
876   )
877   if (length(x = root.groups) != 3) {
878     stop(group.msg)
879   }
880   if (!all(required.groups %in% root.groups)) {
881     stop(group.msg)
882   }
883   unlist(x = sapply(
884     X = required.groups[1:2],
885     FUN = function(group) {
886       if (length(x = list.groups(object = object[[group]], recursive = FALSE)) > 0) {
887         stop(paste("Group", group, "cannot have subgroups"))
888       }
889       if (length(x = list.attributes(object = object[[group]])) > 0) {
890         stop(paste("Group", group, "cannot have subattributes"))
891       }
892       for (dataset in list.datasets(object = object[[group]])) {
893         if (object[[paste(group, dataset, sep = '/')]]$dims != dim.matrix[group]) {
894           stop(paste("All datasets in group", group, "must be of length", required.groups[group]))
895         }
896       }
897     }
898   ))
899   for (dataset in list.datasets(object = object[['/layers']])) {
900     if (any(object[[paste('layers', dataset, sep = '/')]]$dims != dim.matrix)) {
901       stop(paste("All datasets in '/layers' must be", dim.matrix[1], 'by', dim.matrix[2]))
902     }
903   }
904 }
905
906 #' Connect to a loom file
907 #'
908 #' @param filename The loom file to connect to
909 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write
910 #' @param skip.validate Skip the validation steps, use only for extremely large loom files
911 #'
912 #' @return A loom file connection
913 #'
914 #' @seealso \code{\link{loom-class}}
915 #'
916 #' @export
917 #'
918 connect <- function(filename, mode = "r", skip.validate = FALSE) {
919   if (!(mode %in% c('r', 'r+'))) {
920     stop("'mode' must be one of 'r' or 'r+'")
921   }
922   new.loom <- loom$new(filename = filename, mode = mode, skip.validate = skip.validate)
923   return(new.loom)
924 }
925
926 # Map a function or a series of functions over a loom file
927 #
928 # @param X A loom object
929 # @param MARGIN The dimmension to map over, pass 1 for cells or 2 for genes
930 # @param FUN A function to map to the loom file
931 # @param chunk.size Chunk size to use, defaults to \code{loomfile$chunksize[MARGIN]}
932 # @param index.use Indices of the dataset to use, defaults to \code{1:loomfile$shape[MARGIN]}
933 # @param dataset.use Dataset to use, defauts to 'matrix'
934 # @param display.progress Display a progress bar
935 # @param expected Shape of expected results. Can pass either 'matrix' or 'vector'; defaults to shape of 'dataset.use'
936 # @param ... Extra parameters for FUN
937 #
938 # @return The results of the map
939 #
940 # @importFrom utils txtProgressBar setTxtProgressBar
941 #
942 map <- function(
943   X,
944   MARGIN = 1,
945   FUN,
946   chunk.size = NULL,
947   index.use = NULL,
948   dataset.use = 'matrix',
949   display.progress = TRUE,
950   expected = NULL,
951   ...
952 ) {
953   if (!inherits(x = X, what = 'loom')) {
954     stop("map only works on loom objects")
955   }
956   if (!inherits(x = FUN, what = 'function')) {
957     stop("FUN must be a function")
958   }
959   # Check for existance of dataset
960   if (!any(grepl(pattern = dataset.use, x = list.datasets(object = X)))) {
961     stop(paste("Cannot find dataset", dataset.use, "in the loom file"))
962   }
963   # Figure out if we're returning a vector or matrix
964   full.dataset <- grep(
965     pattern = dataset.use,
966     x = list.datasets(object = X),
967     value = TRUE
968   )
969   results.matrix <- TRUE
970   dataset.matrix <- TRUE
971   if (grepl(pattern = 'col_attrs', x = full.dataset)) {
972     MARGIN <- 1
973     results.matrix <- FALSE
974     dataset.matrix <- FALSE
975   } else if (grepl(pattern = 'row_attrs', x = full.dataset)) {
976     MARGIN <- 2
977     results.matrix <- FALSE
978     dataset.matrix <- FALSE
979   }
980   if (!is.null(x = expected)) {
981     results.matrix <- switch(
982       EXPR = expected,
983       'vector' = FALSE,
984       'matrix' = TRUE,
985       stop("'expected' must be one of 'matrix', 'vector', or NULL")
986     )
987   }
988   # Determine the shape of our results
989   if (!(MARGIN %in% c(1, 2))) {
990     stop("MARGIN must be either 1 (cells) or 2 (genes)")
991   }
992   if (is.null(x = index.use)) {
993     index.use <- c(1, X$shape[MARGIN])
994   } else if (length(x = index.use) == 1) {
995     index.use <- c(1, index.use)
996   }
997   index.use[1] <- max(1, index.use[1])
998   index.use[2] <- min(index.use[2], X$shape[MARGIN])
999   batch <- X$batch.scan(
1000     chunk.size = chunk.size,
1001     MARGIN = MARGIN,
1002     index.use = index.use,
1003     dataset.use = dataset.use,
1004     force.reset = TRUE
1005   )
1006   # Create our results holding object
1007   if (results.matrix) {
1008     switch(
1009       EXPR = MARGIN,
1010       '1' = results <- matrix(
1011         nrow = length(x = index.use[1]:index.use[2]),
1012         ncol = X$shape[2]
1013       ),
1014       '2' = results <- matrix(
1015         nrow = X$shape[1],
1016         ncol = length(x = index.use[1]:index.use[2])
1017       )
1018     )
1019   } else {
1020     results <- vector(length = length(x = index.use[1]:index.use[2]))
1021   }
1022   if (display.progress) {
1023     pb <- txtProgressBar(char = '=', style = 3)
1024   }
1025   for (i in 1:length(x = batch)) {
1026     chunk.indices <- X$batch.next(return.data = FALSE)
1027     chunk.data <- if (dataset.matrix) {
1028       switch(
1029         EXPR = MARGIN,
1030         '1' = X[[dataset.use]][chunk.indices, ],
1031         '2' = X[[dataset.use]][, chunk.indices]
1032       )
1033     } else {
1034       X[[dataset.use]][chunk.indices]
1035     }
1036     chunk.data <- FUN(chunk.data, ...)
1037     if (results.matrix) {
1038       if (MARGIN == 1) {
1039         results[chunk.indices, ] <- chunk.data
1040       } else if (MARGIN == 2) {
1041         results[, chunk.indices] <- chunk.data
1042       }
1043     } else {
1044       results[chunk.indices] <- chunk.data
1045     }
1046     if (display.progress) {
1047       setTxtProgressBar(pb = pb, value = i / length(x = batch))
1048     }
1049   }
1050   return(results)
1051 }
1052
1053 # #need to comment
1054 # #need to add progress bar
1055 # #but otherwise, pretty cool
1056 # #for paul to try :
1057 # # f <- connect("~/Downloads/10X43_1.loom")
1058 # # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
1059 # # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
1060 # map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
1061 #   n_func = length(f_list)
1062 #   if (n_func == 1) {
1063 #     f_list <- list(f_list)
1064 #   }
1065 #   if (MARGIN == 1) {
1066 #     results <- list()
1067 #     for (j in 1:n_func) {
1068 #       results[[j]] <- numeric(0)
1069 #     }
1070 #     rows_per_chunk <- chunksize
1071 #     ix <- 1
1072 #     while (ix <= self@shape[1]) {
1073 #       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
1074 #       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
1075 #       for (j in 1:n_func) {
1076 #         new_results <- apply(chunk, 1, FUN = f_list[[j]])
1077 #         results[[j]] <- c(results[[j]], new_results)
1078 #       }
1079 #       ix <- ix + chunksize
1080 #     }
1081 #   }
1082 #   if (MARGIN == 2) {
1083 #     results <- list()
1084 #     for (j in 1:n_func) {
1085 #       results[[j]] <- numeric(0)
1086 #     }
1087 #     cols_per_chunk <- chunksize
1088 #     ix <- 1
1089 #     while (ix <= self@shape[2]) {
1090 #       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
1091 #       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
1092 #       for (j in 1:n_func) {
1093 #         new_results <- apply(chunk, 2, FUN = f_list[[j]])
1094 #         results[[j]] <- c(results[[j]], new_results)
1095 #       }
1096 #       ix <- ix + chunksize
1097 #     }
1098 #   }
1099 #   if (n_func == 1) {
1100 #     results <- results[[1]]
1101 #   }
1102 #   return(results)
1103 # }