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