03dcc15533d9fe6c1ba2eb95802a5437db89d40e
[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, expected, ...)}}{...}
44 #'   \item{\code{map(FUN, MARGIN, chunk.size, index.use, dataset.use, 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[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       index.use = NULL,
296       dataset.use = 'matrix',
297       expected = NULL,
298       ...
299     ) {
300       invisible(x = NULL)
301     },
302     map = function(
303       FUN,
304       MARGIN = 1,
305       chunk.size = NULL,
306       index.use = NULL,
307       dataset.use = 'matrix',
308       display.progress = TRUE,
309       expected = NULL,
310       ...
311     ) {
312       if (!inherits(x = FUN, what = 'function')) {
313         stop("FUN must be a function")
314       }
315       # Checks datset, index, and MARGIN
316       # Sets full dataset path in private$iter.dataset
317       # Sets proper MARGIN in private$iter.margin
318       batch <- self$batch.scan(
319         chunk.size = chunk.size,
320         MARGIN = MARGIN,
321         index.use = index.use,
322         dataset.use = dataset.use,
323         force.reset = TRUE
324       )
325       # Check how we store our results
326       # And what the shape of our dataset is
327       results.matrix <- TRUE
328       dataset.matrix <- TRUE
329       if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
330         results.matrix <- FALSE
331         dataset.matrix <- FALSE
332       } else if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
333         results.matrix <- FALSE
334         dataset.matrix <- FALSE
335       }
336       if (!is.null(x = expected)) {
337         results.matrix <- switch(
338           EXPR = expected,
339           'vector' = FALSE,
340           'matrix' = TRUE,
341           stop("'expected' must be one of 'matrix', 'vector', or NULL")
342         )
343       }
344       # Determine the shape of our results
345       index.use <- private$iter_range(index.use = index.use)
346       # Create our results holding object
347       if (results.matrix) {
348         switch(
349           EXPR = private$iter.margin,
350           '1' = results <- matrix(
351             nrow = length(x = index.use[1]:index.use[2]),
352             ncol = self$shape[2]
353           ),
354           '2' = results <- matrix(
355             nrow = self$shape[1],
356             ncol = length(x = index.use[1]:index.use[2])
357           )
358         )
359       } else {
360         results <- vector(length = length(x = index.use[1]:index.use[2]))
361       }
362       if (display.progress) {
363         pb <- txtProgressBar(char = '=', style = 3)
364       }
365       for (i in 1:length(x = batch)) {
366         chunk.indices <- self$batch.next(return.data = FALSE)
367         chunk.data <- if (dataset.matrix) {
368           switch(
369             EXPR = MARGIN,
370             '1' = self[[dataset.use]][chunk.indices, ],
371             '2' = self[[dataset.use]][, chunk.indices]
372           )
373         } else {
374           self[[dataset.use]][chunk.indices]
375         }
376         chunk.data <- FUN(chunk.data, ...)
377         if (results.matrix) {
378           if (MARGIN == 1) {
379             results[chunk.indices, ] <- chunk.data
380           } else if (MARGIN == 2) {
381             results[, chunk.indices] <- chunk.data
382           }
383         } else {
384           results[chunk.indices] <- chunk.data
385         }
386         if (display.progress) {
387           setTxtProgressBar(pb = pb, value = i / length(x = batch))
388         }
389       }
390       private$reset_batch()
391       return(results)
392     }
393   ),
394   # Private fields and methods
395   # @field err_msg A simple error message if this object hasn't been created with loomR::create or loomR::connect
396   # @field it Iterator object for batch.scan and batch.next
397   # @field iter.dataset Dataset for iterating on
398   # @field iter.margin Margin to iterate over
399   # @field iter.index # Index values for iteration
400   # \describe{
401   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
402   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
403   #   \item{\code{reset_batch()}}{Reset the batch iterator fields}
404   #   \item{\code{iter_range(index.use)}}{Get the range of indices for a batch iteration}
405   # }
406   private = list(
407     # Fields
408     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",
409     it = NULL,
410     iter.dataset = NULL,
411     iter.margin = NULL,
412     iter.index = NULL,
413     # Methods
414     load_attributes = function(MARGIN) {
415       attribute <- switch(
416         EXPR = MARGIN,
417         '1' = 'col_attrs',
418         '2' = 'row_attrs',
419         stop('Invalid attribute dimension')
420       )
421       group <- self[[attribute]]
422       attributes <- unlist(x = lapply(
423         X = names(x = group),
424         FUN = function(x) {
425           d <- list(group[[x]])
426           names(x = d) <- x
427           return(d)
428         }
429       ))
430       if (MARGIN == 1) {
431         self$col.attrs <- attributes
432       } else if (MARGIN == 2) {
433         self$row.attrs <- attributes
434       }
435     },
436     load_layers = function() {
437       self$layers <- unlist(x = lapply(
438         X = names(x = self[['layers']]),
439         FUN = function(n) {
440           d <- c(self[['layers', n]])
441           names(x = d) <- n
442           return(d)
443         }
444       ))
445     },
446     reset_batch = function() {
447       private$it <- NULL
448       private$iter.dataset <- NULL
449       private$iter.margin <- NULL
450       private$iter.index <- NULL
451     },
452     iter_range = function(index.use) {
453       if (is.null(private$iter.margin)) {
454         stop("Batch processing hasn't been set up")
455       }
456       if (is.null(x = index.use)) {
457         # If no index was provided, use entire range for this margin
458         index.use <- c(1, self$shape[private$iter.margin])
459       } else if (length(x = index.use) == 1) {
460         # If one index was provided, start at one and go to index
461         index.use <- c(1, index.use)
462       } else {
463         # Use index.use[1] and index.use[2]
464         index.use <- c(index.use[1], index.use[2])
465       }
466       # Ensure the indices provided fit within the range of the dataset
467       index.use[1] <- max(1, index.use[1])
468       index.use[2] <- min(index.use[2], self$shape[private$iter.margin])
469       # Ensure that index.use[1] is greater than index.use[2]
470       if (index.use[1] > index.use[2]) {
471         stop(paste0(
472           "Starting index (",
473           index.use[1],
474           ") must be lower than the ending index (",
475           index.use[2],
476           ")"
477         ))
478       }
479       return(index.use)
480     }
481   )
482 )
483
484 #' Create a loom object
485 #'
486 #' @param filename The name of the new loom file
487 #' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
488 #' @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}
489 #' @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}
490 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
491 #'
492 #' @return A connection to a loom file
493 #'
494 #' @importFrom utils packageVersion
495 #'
496 #' @seealso \code{\link{loom-class}}
497 #'
498 #' @export
499 #'
500 create <- function(
501   filename,
502   data,
503   gene.attrs = NULL,
504   cell.attrs = NULL,
505   layers = NULL,
506   chunk.dims = 'auto'
507 ) {
508   if (file.exists(filename)) {
509     stop(paste('File', file, 'already exists!'))
510   }
511   if (!is.matrix(x = data)) {
512     data <- as.matrix(x = data)
513   }
514   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
515     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
516   } else if (length(x = chunk.dims == 1)) {
517     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
518       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
519     }
520   } else {
521     chunk.dims <- as.integer(x = chunk.dims)
522   }
523   new.loom <- loom$new(filename = filename, mode = 'w-')
524   # Create the matrix
525   new.loom$create_dataset(
526     name = 'matrix',
527     robj = data,
528     chunk_dims = chunk.dims
529   )
530   new.loom$matrix <- new.loom[['matrix']]
531   new.loom$shape <- new.loom[['matrix']]$dims
532   # Groups
533   new.loom$create_group(name = 'layers')
534   new.loom$create_group(name = 'row_attrs')
535   new.loom$create_group(name = 'col_attrs')
536   # Check for the existance of gene or cell names
537   if (!is.null(x = colnames(x = data))) {
538     new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
539   }
540   if (!is.null(x = rownames(x = data))) {
541     new.loom$add.col.attribute(attribute = list('cell_names' = rownames(x = data)))
542   }
543   # Store some constants as HDF5 attributes
544   h5attr(x = new.loom, which = 'version') <- new.loom$version
545   h5attr(x = new.loom, which = 'chunks') <- paste0(
546     '(',
547     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
548     ')'
549   )
550   # Add layers
551   if (!is.null(x = layers)) {
552     new.loom$add.layer(layer = layers)
553   }
554   if (!is.null(x = gene.attrs)) {
555     new.loom$add.row.attribute(attribute = gene.attrs)
556   }
557   if (!is.null(x = cell.attrs)) {
558     new.loom$add.col.attribute(attribute = cell.attrs)
559   }
560   # Set last bit of information
561   chunks <- new.loom[['matrix']]$chunk_dims
562   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
563   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
564   chunks <- unlist(x = strsplit(x = chunks, split = ','))
565   new.loom$chunksize <- as.integer(x = chunks)
566   # Return the connection
567   return(new.loom)
568 }
569
570 #' Validate a loom object
571 #'
572 #' @param object A loom object
573 #'
574 #' @return None, errors out if object is an invalid loom connection
575 #'
576 #' @seealso \code{\link{loom-class}}
577 #'
578 #' @export
579 #'
580 validateLoom <- function(object) {
581   # A loom file is a specific HDF5
582   # We need a dataset in /matrix that's a two-dimensional dense matrix
583   root.datasets <- list.datasets(object = object, path = '/', recursive = FALSE)
584   if (length(x = root.datasets) != 1) {
585     stop("There can only be one dataset at the root of the loom file")
586   }
587   if (root.datasets != 'matrix') {
588     stop("The root dataset must be called '/matrix'")
589   }
590   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
591   required.groups <- c('row_attrs', 'col_attrs', 'layers')
592   dim.matrix <- object[[root.datasets]]$dims # Columns x Rows
593   names(x = dim.matrix) <- required.groups[c(2, 1)]
594   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
595   group.msg <- paste0(
596     "There can only be three groups in the loom file: '",
597     paste(required.groups, collapse = "', '"),
598     "'"
599   )
600   if (length(x = root.groups) != 3) {
601     stop(group.msg)
602   }
603   if (!all(required.groups %in% root.groups)) {
604     stop(group.msg)
605   }
606   unlist(x = sapply(
607     X = required.groups[1:2],
608     FUN = function(group) {
609       if (length(x = list.groups(object = object[[group]], recursive = FALSE)) > 0) {
610         stop(paste("Group", group, "cannot have subgroups"))
611       }
612       if (length(x = list.attributes(object = object[[group]])) > 0) {
613         stop(paste("Group", group, "cannot have subattributes"))
614       }
615       for (dataset in list.datasets(object = object[[group]])) {
616         if (object[[paste(group, dataset, sep = '/')]]$dims != dim.matrix[group]) {
617           stop(paste("All datasets in group", group, "must be of length", required.groups[group]))
618         }
619       }
620     }
621   ))
622   for (dataset in list.datasets(object = object[['/layers']])) {
623     if (any(object[[paste('layers', dataset, sep = '/')]]$dims != dim.matrix)) {
624       stop(paste("All datasets in '/layers' must be", dim.matrix[1], 'by', dim.matrix[2]))
625     }
626   }
627 }
628
629 #' Connect to a loom file
630 #'
631 #' @param filename The loom file to connect to
632 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write
633 #'
634 #' @return A loom file connection
635 #'
636 #' @seealso \code{\link{loom-class}}
637 #'
638 #' @export
639 #'
640 connect <- function(filename, mode = "r") {
641   if (!(mode %in% c('r', 'r+'))) {
642     stop("'mode' must be one of 'r' or 'r+'")
643   }
644   new.loom <- loom$new(filename = filename, mode = mode)
645   return(new.loom)
646 }
647
648 CreateLoomFromSeurat <- function(object, filename) {
649   object.data=t(object@raw.data[rownames(object@data),object@cell.names])
650   object.meta.data=object@meta.data
651   row_attrs=list(); col_attrs=list()
652   gene.names=colnames(object.data)
653   object.meta.data$ident = object@ident
654   object.meta.data$CellID = object@cell.names
655   for(i in 1:ncol(object.meta.data)) {
656     col_attrs[[colnames(object.meta.data)[i]]]=object.meta.data[,i]
657   }
658   row_attrs[["Gene"]]=gene.names
659   create(filename,object.data,gene.attrs = row_attrs, cell.attrs = col_attrs)
660 }
661
662 #' Map a function or a series of functions over a loom file
663 #'
664 #' @param X A loom object
665 #' @param MARGIN The dimmension to map over, pass 1 for cells or 2 for genes
666 #' @param FUN A function to map to the loom file
667 #' @param chunk.size Chunk size to use, defaults to \code{loomfile$chunksize[MARGIN]}
668 #' @param index.use Indices of the dataset to use, defaults to \code{1:loomfile$shape[MARGIN]}
669 #' @param dataset.use Dataset to use, defauts to 'matrix'
670 #' @param display.progress Display a progress bar
671 #' @param expected Shape of expected results. Can pass either 'matrix' or 'vector'; defaults to shape of 'dataset.use'
672 #' @param ... Extra parameters for FUN
673 #'
674 #' @return The results of the map
675 #'
676 #' @importFrom utils txtProgressBar setTxtProgressBar
677 #'
678 #' @export
679 #'
680 map <- function(
681   X,
682   MARGIN = 1,
683   FUN,
684   chunk.size = NULL,
685   index.use = NULL,
686   dataset.use = 'matrix',
687   display.progress = TRUE,
688   expected = NULL,
689   ...
690 ) {
691   if (!inherits(x = X, what = 'loom')) {
692     stop("map only works on loom objects")
693   }
694   if (!inherits(x = FUN, what = 'function')) {
695     stop("FUN must be a function")
696   }
697   # Check for existance of dataset
698   if (!any(grepl(pattern = dataset.use, x = list.datasets(object = X)))) {
699     stop(paste("Cannot find dataset", dataset.use, "in the loom file"))
700   }
701   # Figure out if we're returning a vector or matrix
702   full.dataset <- grep(
703     pattern = dataset.use,
704     x = list.datasets(object = X),
705     value = TRUE
706   )
707   results.matrix <- TRUE
708   dataset.matrix <- TRUE
709   if (grepl(pattern = 'col_attrs', x = full.dataset)) {
710     MARGIN <- 1
711     results.matrix <- FALSE
712     dataset.matrix <- FALSE
713   } else if (grepl(pattern = 'row_attrs', x = full.dataset)) {
714     MARGIN <- 2
715     results.matrix <- FALSE
716     dataset.matrix <- FALSE
717   }
718   if (!is.null(x = expected)) {
719     results.matrix <- switch(
720       EXPR = expected,
721       'vector' = FALSE,
722       'matrix' = TRUE,
723       stop("'expected' must be one of 'matrix', 'vector', or NULL")
724     )
725   }
726   # Determine the shape of our results
727   if (!(MARGIN %in% c(1, 2))) {
728     stop("MARGIN must be either 1 (cells) or 2 (genes)")
729   }
730   if (is.null(x = index.use)) {
731     index.use <- c(1, X$shape[MARGIN])
732   } else if (length(x = index.use) == 1) {
733     index.use <- c(1, index.use)
734   }
735   index.use[1] <- max(1, index.use[1])
736   index.use[2] <- min(index.use[2], X$shape[MARGIN])
737   batch <- X$batch.scan(
738     chunk.size = chunk.size,
739     MARGIN = MARGIN,
740     index.use = index.use,
741     dataset.use = dataset.use,
742     force.reset = TRUE
743   )
744   # Create our results holding object
745   if (results.matrix) {
746     switch(
747       EXPR = MARGIN,
748       '1' = results <- matrix(
749         nrow = length(x = index.use[1]:index.use[2]),
750         ncol = X$shape[2]
751       ),
752       '2' = results <- matrix(
753         nrow = X$shape[1],
754         ncol = length(x = index.use[1]:index.use[2])
755       )
756     )
757   } else {
758     results <- vector(length = length(x = index.use[1]:index.use[2]))
759   }
760   if (display.progress) {
761     pb <- txtProgressBar(char = '=', style = 3)
762   }
763   for (i in 1:length(x = batch)) {
764     chunk.indices <- X$batch.next(return.data = FALSE)
765     chunk.data <- if (dataset.matrix) {
766       switch(
767         EXPR = MARGIN,
768         '1' = X[[dataset.use]][chunk.indices, ],
769         '2' = X[[dataset.use]][, chunk.indices]
770       )
771     } else {
772       X[[dataset.use]][chunk.indices]
773     }
774     chunk.data <- FUN(chunk.data, ...)
775     if (results.matrix) {
776       if (MARGIN == 1) {
777         results[chunk.indices, ] <- chunk.data
778       } else if (MARGIN == 2) {
779         results[, chunk.indices] <- chunk.data
780       }
781     } else {
782       results[chunk.indices] <- chunk.data
783     }
784     if (display.progress) {
785       setTxtProgressBar(pb = pb, value = i / length(x = batch))
786     }
787   }
788   return(results)
789 }
790
791 #need to comment
792 #need to add progress bar
793 #but otherwise, pretty cool
794 #for paul to try :
795 # f <- connect("~/Downloads/10X43_1.loom")
796 # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
797 # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
798 map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
799   n_func = length(f_list)
800   if (n_func == 1) {
801     f_list <- list(f_list)
802   }
803   if (MARGIN == 1) {
804     results <- list()
805     for (j in 1:n_func) {
806       results[[j]] <- numeric(0)
807     }
808     rows_per_chunk <- chunksize
809     ix <- 1
810     while (ix <= self@shape[1]) {
811       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
812       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
813       for (j in 1:n_func) {
814         new_results <- apply(chunk, 1, FUN = f_list[[j]])
815         results[[j]] <- c(results[[j]], new_results)
816       }
817       ix <- ix + chunksize
818     }
819   }
820   if (MARGIN == 2) {
821     results <- list()
822     for (j in 1:n_func) {
823       results[[j]] <- numeric(0)
824     }
825     cols_per_chunk <- chunksize
826     ix <- 1
827     while (ix <= self@shape[2]) {
828       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
829       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
830       for (j in 1:n_func) {
831         new_results <- apply(chunk, 2, FUN = f_list[[j]])
832         results[[j]] <- c(results[[j]], new_results)
833       }
834       ix <- ix + chunksize
835     }
836   }
837   if (n_func == 1) {
838     results <- results[[1]]
839   }
840   return(results)
841 }