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