Fix issues with create, further progress on batching
[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 #' @return Object of \code{\link{R6::R6Class}} to generate \code{loom} objects
11 #' @format An \code{\link{R6::R6Class}} object
12 #' @seealso \code{\link{hdf5r::H5File}}
13 #'
14 #' @field version Version of loomR object was created under
15 #' @field shape Shape of \code{/matrix} in columns (cells) by rows (genes)
16 #' @field chunksize Chunks set for this dataset in columns (cells) by rows (genes)
17 #' @field matrix The main data matrix, stored as columns (cells) by rows (genes)
18 #' @field layers Additional data matricies, the same shape as \code{/matrix}
19 #' @field col.attrs Extra information about cells
20 #' @field row.attrs Extra information about genes
21 #'
22 #' @section Methods:
23 #' \describe{
24 #'   \item{\code{add.layer(layer)}}{Add a data layer to this loom file, must be in column (cells) by row (genes) orientation}
25 #'   \item{\code{add.attribute(attribute, MARGIN)}}{
26 #'     Add extra information to this loom file where
27 #'     \code{attribute} is a named list where each element is a vector that is as long as one dimension of \code{/matrix} and
28 #'     \code{MARGIN} is either 1 for cells or 2 for genes
29 #'   }
30 #'   \item{\code{add.row.attribute(attribute)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 2)}}
31 #'   \item{\code{add.col.attribute(attribute)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 1)}}
32 #'   \item{\code{add.meta.data(meta.data)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 1)}}
33 #'   \item{\code{batch.scan(chunk.size, MARGIN, index.use, dataset.use)}, \code{batch.next}}{
34 #'     Scan a dataset in the loom file from \code{index.use[1]} to \code{index.use[2]}, iterating by \code{chunk.size}.
35 #'     \code{dataset.use} can be the name, not \code{group/name}, unless the name is present in multiple groups.
36 #'     If the dataset is in col_attrs, pass \code{MARGIN = 1}; if in row_attrs, pass \code{MARGIN = 2}.
37 #'     Otherwise, pass \code{MARGIN = 1} to iterate on cells or \code{MARGIN = 2} to iterate on genes.
38 #'     \code{chunk.size} defaults to \code{self$chunksize}, \code{MARGIN} defaults to 1,
39 #'     \code{index.use} defaults to \code{1:self$shape[MARGIN]}, \code{dataset.use} defaults to 'matrix'
40 #'   }
41 #' }
42 #'
43 #' @importFrom iterators nextElem
44 #' @importFrom utils packageVersion
45 #' @importFrom itertools ihasNext ichunk hasNext
46 #'
47 #' @export
48 #'
49 loom <- R6Class(
50   # The loom class
51   # Based on the H5File class from hdf5r
52   # Not clonable (no loom$clone method), doesn't make sense since all data is on disk, not in memory
53   # Yes to portability, other packages can subclass the loom class
54   # Class is locked, other fields and methods cannot be added
55   classname = 'loom',
56   inherit = hdf5r::H5File,
57   cloneable = FALSE,
58   portable = TRUE,
59   lock_class = TRUE,
60   # Public fields and methods
61   # See above for documentation
62   public = list(
63     # Fields
64     version = NULL,
65     shape = NULL,
66     chunksize = NULL,
67     matrix = NULL,
68     layers = NULL,
69     col.attrs = NULL,
70     row.attrs = NULL,
71     # Methods
72     initialize = function(filename = NULL, mode = c('a', 'r', 'r+', 'w', 'w-'), ...) {
73       # If the file exists, run validation steps
74       do.validate <- file.exists(filename) && !(mode %in% c('w', 'w+'))
75       super$initialize(filename = filename, mode = mode, ...)
76       if (do.validate) {
77         # Run the validation steps
78         validateLoom(object = self)
79         # Store /matrix and the shape of /matrix
80         self$matrix <- self[['matrix']]
81         self$shape <- self[['matrix']]$dims
82         # Store the chunk size
83         chunks <- h5attr(x = self, which = 'chunks')
84         chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
85         chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
86         chunks <- unlist(x = strsplit(x = chunks, split = ','))
87         self$chunksize <- as.integer(x = chunks)
88         # Store version information
89         self$version <- as.character(x = tryCatch(
90           # Try getting a version
91           # If it doesn't exist, can we write to the file?
92           # If so, store the version as this version of loomR
93           # Otherwise, store the version as NA_character_
94           expr = h5attr(x = self, which = 'version'),
95           error = function(e) {
96             if (mode != 'r') {
97               version <- packageVersion(pkg = 'loomR')
98               h5attr(x = self, which = 'version') <- version
99             } else {
100               version <- NA_character_
101             }
102             return(version)
103           }
104         ))
105         # Load layers
106         private$load_layers()
107         # Load attributes
108         private$load_attributes(MARGIN = 1) # Cells (col_attrs)
109         private$load_attributes(MARGIN = 2) # Genes (row_attrs)
110       } else {
111         # Assume new HDF5 file
112         self$version <- as.character(x = packageVersion(pkg = 'loomR'))
113       }
114     },
115     add.layer = function(layers) {
116       # Value checking
117       if (!is.list(x = layers) || is.null(x = names(x = layers))) {
118         stop("'layers' must be a named list")
119       }
120       if (is.null(x = self$shape)) {
121         stop(private$err_msg)
122       }
123       # Add layers
124       for (i in 1:length(x = layers)) {
125         if (!is.matrix(x = layers[[i]])) {
126           layers[[i]] <- as.matrix(x = layers[[i]])
127         }
128         do.transpose <- FALSE
129         if (any(dim(x = layers[[i]]) != self$shape)) {
130           if (all(rev(x = dim(x = layers[[i]])) == self$shape)) {
131             do.transpose <- TRUE
132           } else {
133             stop(paste(
134               "All layers must have",
135               self$shape[1],
136               "rows for cells and",
137               self$shape[2],
138               "columns for genes"
139             ))
140           }
141         }
142         if (do.transpose) {
143           self[['layers', names(x = layers)[i]]] <- t(x = layers[[i]])
144         } else {
145           self[['layers', names(x = layers)[i]]] <- layers[[i]]
146         }
147       }
148       self$flush()
149       private$load_layers()
150       invisible(x = self)
151     },
152     add.attribute = function(attribute, MARGIN) {
153       # Value checking
154       if (!is.list(x = attribute) || is.null(x = names(x = attribute))) {
155         stop("'attribute' must be a named list")
156       }
157       if (length(x = attribute) > 1) {
158         for (i in attribute) {
159           if (!is.vector(x = attribute)) {
160             stop("All attributes must be one-dimensional vectors")
161           }
162         }
163       }
164       if (length(x = which(x = names(x = attribute) != '')) != length(x = attribute)) {
165         stop("Not all attributes had names provided")
166       }
167       if (!MARGIN %in% c(1, 2)) {
168         stop("'MARGIN' must be 1 or 2")
169       }
170       # Add the attributes as datasets for our MARGIN's group
171       if (is.null(x = self$shape)) {
172         stop(private$err_msg)
173       }
174       grp.name <- c('col_attrs', 'row_attrs')[MARGIN]
175       grp <- self[[grp.name]]
176       for (i in 1:length(x = attribute)) {
177         if (length(attribute[[i]]) != self$shape[MARGIN])
178           stop(paste(
179             "All",
180             switch(EXPR = MARGIN, '1' = 'cell', '2' = 'gene'),
181             "attributes must be of length",
182             self$shape[MARGIN]
183           ))
184         grp[[names(x = attribute)[i]]] <- attribute[[i]]
185       }
186       self$flush()
187       gc(verbose = FALSE)
188       # Load the attributes for this margin
189       private$load_attributes(MARGIN = MARGIN)
190       invisible(x = self)
191     },
192     # Add attributes for genes
193     add.row.attribute = function(attribute) {
194       self$add.attribute(attribute = attribute, MARGIN = 2)
195       invisible(x = self)
196     },
197     # Add attributes for cells
198     add.col.attribute = function(attribute) {
199       self$add.attribute(attribute = attribute, MARGIN = 1)
200       invisible(x = self)
201     },
202     # Add metadata, follows cells
203     add.meta.data = function(meta.data) {
204       self$add.col.attribute(attribute = meta.data)
205       invisible(x = self)
206     },
207     # Batch scan
208     batch.scan = function(
209       chunk.size = NULL,
210       MARGIN = 1,
211       index.use = NULL,
212       dataset.use = 'matrix'
213     ) {
214       if (is.null(x = private$it) || !grepl(pattern = dataset.use, x = private$iter.dataset)) {
215         # Check the existence of the dataset
216         private$iter.dataset <- grep(
217           pattern = dataset.use,
218           x = list.datasets(object = self),
219           value = TRUE
220         )
221         if (length(x = private$iter.dataset) != 1) {
222           stop(paste0("Cannot find dataset '", dataset.use, "' in the loom file"))
223         }
224         # Check the margin
225         if (!(MARGIN %in% c(1, 2))) {
226           stop("MARGIN must be 1 (cells) or 2 (genes)")
227         } else {
228           private$iter.margin <- MARGIN
229         }
230         if (is.null(x = chunk.size)) {
231           chunk.size <- self$chunksize[MARGIN]
232         }
233         # Set the indices to use
234         if (is.null(x = index.use)) {
235           index.use <- c(1, self$shape[MARGIN])
236         } else if (length(x = index.use) == 1) {
237           index.use <- index.use <- 1:index.use
238         } else {
239           index.use <- c(index.use[1], index.use[2])
240         }
241         index.use[1] <- max(1, index.use[1])
242         index.use[2] <- min(index.use[2], self$shape[MARGIN])
243         if (index.use[1] > index.use[2]) {
244           stop(paste0(
245             "Starting index (",
246             index.use[1],
247             ") must be lower than the ending index (",
248             index.use[2],
249             ")"
250           ))
251         }
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(private$iter.index[1]:private$iter.index[2])
260       # # Do the iterating
261       # if (hasNext(obj = private$it)) {
262       #   chunk.indices <- unlist(x = nextElem(obj = private$it))
263       #   if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
264       #     return(switch(
265       #       EXPR = private$iter.margin,
266       #       '1' = self[[private$iter.dataset]][chunk.indices, ],
267       #       '2' = self[[private$iter.dataset]][, chunk.indices]
268       #     ))
269       #   } else {
270       #     return(self[[private$iter.dataset]][chunk.indices])
271       #   }
272       # } else {
273       #   private$it <- NULL
274       #   return(NULL)
275       # }
276     },
277     batch.next = function() {
278       if (!'hasNext.ihasNext' %in% suppressWarnings(expr = methods(class = class(x = private$it)))) {
279         stop("Please setup the iterator with self$batch.scan")
280       }
281       # Do the iterating
282       if (hasNext(obj = private$it)) {
283         chunk.indices <- unlist(x = nextElem(obj = private$it))
284         if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
285           to.return <- switch(
286             EXPR = private$iter.margin,
287             '1' = self[[private$iter.dataset]][chunk.indices, ],
288             '2' = self[[private$iter.dataset]][, chunk.indices]
289           )
290         } else {
291           to.return <- self[[private$iter.dataset]][chunk.indices]
292         }
293         if (!hasNext(obj = private$it)) {
294           private$reset_batch()
295         }
296         return(to.return)
297       } else {
298         private$reset_batch()
299         return(NULL)
300       }
301     }
302   ),
303   # Private fields and methods
304   # @field err_msg A simple error message if this object hasn't been created with loomR::create or loomR::connect
305   # @field it
306   # @field iter.dataset
307   # @field iter.margin
308   # @field iter.index
309   # \describe{
310   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
311   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
312   #   \item{\code{reset_batch()}}{Reset the batch iterator fields}
313   # }
314   private = list(
315     # Fields
316     err_msg = "This loom object has not been created with either loomR::create or loomR::connect, please use these function to create or connect to a loom file",
317     it = NULL,
318     iter.dataset = NULL,
319     iter.margin = NULL,
320     iter.index = NULL,
321     # Methods
322     load_attributes = function(MARGIN) {
323       attribute <- switch(
324         EXPR = MARGIN,
325         '1' = 'col_attrs',
326         '2' = 'row_attrs',
327         stop('Invalid attribute dimension')
328       )
329       group <- self[[attribute]]
330       attributes <- unlist(x = lapply(
331         X = names(x = group),
332         FUN = function(x) {
333           d <- list(group[[x]])
334           names(x = d) <- x
335           return(d)
336         }
337       ))
338       if (MARGIN == 1) {
339         self$col.attrs <- attributes
340       } else if (MARGIN == 2) {
341         self$row.attrs <- attributes
342       }
343     },
344     load_layers = function() {
345       self$layers <- unlist(x = lapply(
346         X = names(x = self[['layers']]),
347         FUN = function(n) {
348           d <- c(self[['layers', n]])
349           names(x = d) <- n
350           return(d)
351         }
352       ))
353     },
354     reset_batch = function() {
355       private$it <- NULL
356       private$iter.dataset <- NULL
357       private$iter.margin <- NULL
358       private$iter.index <- NULL
359     }
360   )
361 )
362
363 #' Create a loom object
364 #'
365 #' @param filename The name of the new loom file
366 #' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
367 #' @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}
368 #' @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}
369 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
370 #'
371 #' @return A connection to a loom file
372 #'
373 #' @importFrom utils packageVersion
374 #'
375 #' @seealso \code{\link{loom-class}}
376 #'
377 #' @export
378 #'
379 create <- function(
380   filename,
381   data,
382   gene.attrs = NULL,
383   cell.attrs = NULL,
384   layers = NULL,
385   chunk.dims = 'auto'
386 ) {
387   if (file.exists(filename)) {
388     stop(paste('File', file, 'already exists!'))
389   }
390   if (!is.matrix(x = data)) {
391     data <- as.matrix(x = data)
392   }
393   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
394     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
395   } else if (length(x = chunk.dims == 1)) {
396     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
397       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
398     }
399   } else {
400     chunk.dims <- as.integer(x = chunk.dims)
401   }
402   new.loom <- loom$new(filename = filename, mode = 'w-')
403   # Create the matrix
404   new.loom$create_dataset(
405     name = 'matrix',
406     robj = data,
407     chunk_dims = chunk.dims
408   )
409   new.loom$matrix <- new.loom[['matrix']]
410   new.loom$shape <- new.loom[['matrix']]$dims
411   # Groups
412   new.loom$create_group(name = 'layers')
413   new.loom$create_group(name = 'row_attrs')
414   new.loom$create_group(name = 'col_attrs')
415   # Check for the existance of gene or cell names
416   if (!is.null(x = colnames(x = data))) {
417     new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
418   }
419   if (!is.null(x = rownames(x = data))) {
420     new.loom$add.col.attribute(attribute = list('cell_names' = rownames(x = data)))
421   }
422   # Store some constants as HDF5 attributes
423   h5attr(x = new.loom, which = 'version') <- new.loom$version
424   h5attr(x = new.loom, which = 'chunks') <- paste0(
425     '(',
426     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
427     ')'
428   )
429   # Add layers
430   if (!is.null(x = layers)) {
431     new.loom$add.layer(layer = layers)
432   }
433   if (!is.null(x = gene.attrs)) {
434     new.loom$add.row.attribute(attribute = gene.attrs)
435   }
436   if (!is.null(x = cell.attrs)) {
437     new.loom$add.col.attribute(attribute = cell.attrs)
438   }
439   # Set last bit of information
440   chunks <- new.loom[['matrix']]$chunk_dims
441   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
442   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
443   chunks <- unlist(x = strsplit(x = chunks, split = ','))
444   new.loom$chunksize <- as.integer(x = chunks)
445   # Return the connection
446   return(new.loom)
447 }
448
449 #' Validate a loom object
450 #'
451 #' @param object A loom object
452 #'
453 #' @return None, errors out if object is an invalid loom connection
454 #'
455 #' @seealso \code{\link{loom-class}}
456 #'
457 #' @export
458 #'
459 validateLoom <- function(object) {
460   # A loom file is a specific HDF5
461   # We need a dataset in /matrix that's a two-dimensional dense matrix
462   root.datasets <- list.datasets(object = object, path = '/', recursive = FALSE)
463   if (length(x = root.datasets) != 1) {
464     stop("There can only be one dataset at the root of the loom file")
465   }
466   if (root.datasets != 'matrix') {
467     stop("The root dataset must be called '/matrix'")
468   }
469   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
470   required.groups <- c('row_attrs', 'col_attrs', 'layers')
471   dim.matrix <- object[[root.datasets]]$dims # Columns x Rows
472   names(x = dim.matrix) <- required.groups[c(2, 1)]
473   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
474   group.msg <- paste0(
475     "There can only be three groups in the loom file: '",
476     paste(required.groups, collapse = "', '"),
477     "'"
478   )
479   if (length(x = root.groups) != 3) {
480     stop(group.msg)
481   }
482   if (!all(required.groups %in% root.groups)) {
483     stop(group.msg)
484   }
485   unlist(x = sapply(
486     X = required.groups[1:2],
487     FUN = function(group) {
488       if (length(x = list.groups(object = object[[group]], recursive = FALSE)) > 0) {
489         stop(paste("Group", group, "cannot have subgroups"))
490       }
491       if (length(x = list.attributes(object = object[[group]])) > 0) {
492         stop(paste("Group", group, "cannot have subattributes"))
493       }
494       for (dataset in list.datasets(object = object[[group]])) {
495         if (object[[paste(group, dataset, sep = '/')]]$dims != dim.matrix[group]) {
496           stop(paste("All datasets in group", group, "must be of length", required.groups[group]))
497         }
498       }
499     }
500   ))
501   for (dataset in list.datasets(object = object[['/layers']])) {
502     if (any(object[[paste('layers', dataset, sep = '/')]]$dims != dim.matrix)) {
503       stop(paste("All datasets in '/layers' must be", dim.matrix[1], 'by', dim.matrix[2]))
504     }
505   }
506 }
507
508 #' Connect to a loom file
509 #'
510 #' @param filename The loom file to connect to
511 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write
512 #'
513 #' @return A loom file connection
514 #'
515 #' @seealso \code{\link{loom-class}}
516 #'
517 #' @export
518 #'
519 connect <- function(filename, mode = "r") {
520   if (!(mode %in% c('r', 'r+'))) {
521     stop("'mode' must be one of 'r' or 'r+'")
522   }
523   new.loom <- loom$new(filename = filename, mode = mode)
524   return(new.loom)
525 }
526
527 CreateLoomFromSeurat <- function(object, filename) {
528   object.data=t(object@raw.data[rownames(object@data),object@cell.names])
529   object.meta.data=object@meta.data
530   row_attrs=list(); col_attrs=list()
531   gene.names=colnames(object.data)
532   object.meta.data$ident = object@ident
533   object.meta.data$CellID = object@cell.names
534   for(i in 1:ncol(object.meta.data)) {
535     col_attrs[[colnames(object.meta.data)[i]]]=object.meta.data[,i]
536   }
537   row_attrs[["Gene"]]=gene.names
538   create(filename,object.data,gene.attrs = row_attrs, cell.attrs = col_attrs)
539 }
540
541 #need to comment
542 #need to add progress bar
543 #but otherwise, pretty cool
544 #for paul to try :
545 # f <- connect("~/Downloads/10X43_1.loom")
546 # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
547 # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
548 map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
549   n_func = length(f_list)
550   if (n_func == 1) {
551     f_list <- list(f_list)
552   }
553   if (MARGIN == 1) {
554     results <- list()
555     for (j in 1:n_func) {
556       results[[j]] <- numeric(0)
557     }
558     rows_per_chunk <- chunksize
559     ix <- 1
560     while (ix <= self@shape[1]) {
561       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
562       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
563       for (j in 1:n_func) {
564         new_results <- apply(chunk, 1, FUN = f_list[[j]])
565         results[[j]] <- c(results[[j]], new_results)
566       }
567       ix <- ix + chunksize
568     }
569   }
570   if (MARGIN == 2) {
571     results <- list()
572     for (j in 1:n_func) {
573       results[[j]] <- numeric(0)
574     }
575     cols_per_chunk <- chunksize
576     ix <- 1
577     while (ix <= self@shape[2]) {
578       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
579       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
580       for (j in 1:n_func) {
581         new_results <- apply(chunk, 2, FUN = f_list[[j]])
582         results[[j]] <- c(results[[j]], new_results)
583       }
584       ix <- ix + chunksize
585     }
586   }
587   if (n_func == 1) {
588     results <- results[[1]]
589   }
590   return(results)
591 }