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