c5515c6af1bfd766ac346398f15ab46d666a7bf4
[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)}}{
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 <- 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 <- 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         # Setup our iterator
244         private$it <- ihasNext(iterable = ichunk(
245           iterable = index.use[1]:index.use[2],
246           chunkSize = chunk.size
247         ))
248       }
249       # Do the iterating
250       if (hasNext(obj = private$it)) {
251         chunk.indices <- unlist(x = nextElem(obj = private$it))
252         if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
253           return(switch(
254             EXPR = private$iter.margin,
255             '1' = self[[private$iter.dataset]][chunk.indices, ],
256             '2' = self[[private$iter.dataset]][, chunk.indices]
257           ))
258         } else {
259           return(self[[private$iter.dataset]][chunk.indices])
260         }
261       } else {
262         private$it <- NULL
263         return(NULL)
264       }
265     }
266   ),
267   # Private fields and methods
268   # @field err_msg A simple error message if this object hasn't been created with loomR::create or loomR::connect
269   # \describe{
270   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
271   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
272   # }
273   private = list(
274     # Fields
275     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",
276     it = NULL,
277     iter.dataset = NULL,
278     iter.margin = NULL,
279     # Methods
280     load_attributes = function(MARGIN) {
281       attribute <- switch(
282         EXPR = MARGIN,
283         '1' = 'col_attrs',
284         '2' = 'row_attrs',
285         stop('Invalid attribute dimension')
286       )
287       group <- self[[attribute]]
288       attributes <- unlist(x = lapply(
289         X = names(x = group),
290         FUN = function(x) {
291           d <- list(group[[x]])
292           names(x = d) <- x
293           return(d)
294         }
295       ))
296       if (MARGIN == 1) {
297         self$col.attrs <- attributes
298       } else if (MARGIN == 2) {
299         self$row.attrs <- attributes
300       }
301     },
302     load_layers = function() {
303       self$layers <- unlist(x = lapply(
304         X = names(x = self[['layers']]),
305         FUN = function(n) {
306           d <- c(self[['layers', n]])
307           names(x = d) <- n
308           return(d)
309         }
310       ))
311     }
312   )
313 )
314
315 #' Create a loom object
316 #'
317 #' @param filename The name of the new loom file
318 #' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
319 #' @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}
320 #' @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}
321 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
322 #'
323 #' @return A connection to a loom file
324 #'
325 #' @importFrom utils packageVersion
326 #'
327 #' @seealso \code{\link{loom-class}}
328 #'
329 #' @export
330 #'
331 create <- function(
332   filename,
333   data,
334   gene.attrs = NULL,
335   cell.attrs = NULL,
336   layers = NULL,
337   chunk.dims = 'auto'
338 ) {
339   if (file.exists(filename)) {
340     stop(paste('File', file, 'already exists!'))
341   }
342   if (!is.matrix(x = data)) {
343     data <- as.matrix(x = data)
344   }
345   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
346     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
347   } else if (length(x = chunk.dims == 1)) {
348     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
349       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
350     }
351   } else {
352     chunk.dims <- as.integer(x = chunk.dims)
353   }
354   new.loom <- loom$new(filename = filename, mode = 'w-')
355   # Create the matrix
356   new.loom$create_dataset(
357     name = 'matrix',
358     robj = t(x = data),
359     chunk_dims = chunk.dims
360   )
361   new.loom$shape <- new.loom[['matrix']]
362   new.loom$shape <- new.loom[['matrix']]$dims
363   if (!is.null(x = colnames(x = data))) {
364     new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
365   }
366   if (!is.null(x = rownames(x = data))) {
367     new.loom$add.col.attribute(attribute = list('cell_names' = colnames(x = data)))
368   }
369   # Store some constants as HDF5 attributes
370   h5attr(x = new.loom, which = 'version') <- new.loom$version
371   h5attr(x = new.loom, which = 'chunks') <- paste0(
372     '(',
373     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
374     ')'
375   )
376   # Groups
377   new.loom$create_group(name = 'layers')
378   new.loom$create_group(name = 'row_attrs')
379   new.loom$create_group(name = 'col_attrs')
380   # Add layers
381   for (ly in layers) {
382     new.loom$add.layer(layer = ly)
383   }
384   if (!is.null(x = gene.attrs)) {
385     new.loom$add.row.attribute(attribute = gene.attrs)
386   }
387   if (!is.null(x = cell.attrs)) {
388     new.loom$add.col.attribute(attribute = cell.attrs)
389   }
390   # Set last bit of information
391   chunks <- new.loom[['matrix']]$chunk_dims
392   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
393   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
394   chunks <- unlist(x = strsplit(x = chunks, split = ','))
395   new.loom$chunksize <- as.integer(x = chunks)
396   # Return the connection
397   return(new.loom)
398 }
399
400 #' Validate a loom object
401 #'
402 #' @param object A loom object
403 #'
404 #' @return None, errors out if object is an invalid loom connection
405 #'
406 #' @seealso \code{\link{loom-class}}
407 #'
408 #' @export
409 #'
410 validateLoom <- function(object) {
411   # A loom file is a specific HDF5
412   # We need a dataset in /matrix that's a two-dimensional dense matrix
413   root.datasets <- list.datasets(object = object, path = '/', recursive = FALSE)
414   if (length(x = root.datasets) != 1) {
415     stop("There can only be one dataset at the root of the loom file")
416   }
417   if (root.datasets != 'matrix') {
418     stop("The root dataset must be called '/matrix'")
419   }
420   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
421   required.groups <- c('row_attrs', 'col_attrs', 'layers')
422   dim.matrix <- object[[root.datasets]]$dims # Columns x Rows
423   names(x = dim.matrix) <- required.groups[c(2, 1)]
424   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
425   group.msg <- paste0(
426     "There can only be three groups in the loom file: '",
427     paste(required.groups, collapse = "', '"),
428     "'"
429   )
430   if (length(x = root.groups) != 3) {
431     stop(group.msg)
432   }
433   if (!all(required.groups %in% root.groups)) {
434     stop(group.msg)
435   }
436   unlist(x = sapply(
437     X = required.groups[1:2],
438     FUN = function(group) {
439       if (length(x = list.groups(object = object[[group]], recursive = FALSE)) > 0) {
440         stop(paste("Group", group, "cannot have subgroups"))
441       }
442       if (length(x = list.attributes(object = object[[group]])) > 0) {
443         stop(paste("Group", group, "cannot have subattributes"))
444       }
445       for (dataset in list.datasets(object = object[[group]])) {
446         if (object[[paste(group, dataset, sep = '/')]]$dims != dim.matrix[group]) {
447           stop(paste("All datasets in group", group, "must be of length", required.groups[group]))
448         }
449       }
450     }
451   ))
452   for (dataset in list.datasets(object = object[['/layers']])) {
453     if (any(object[[paste('layers', dataset, sep = '/')]]$dims != dim.matrix)) {
454       stop(paste("All datasets in '/layers' must be", dim.matrix[1], 'by', dim.matrix[2]))
455     }
456   }
457 }
458
459 #' Connect to a loom file
460 #'
461 #' @param filename The loom file to connect to
462 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write
463 #'
464 #' @return A loom file connection
465 #'
466 #' @seealso \code{\link{loom-class}}
467 #'
468 #' @export
469 #'
470 connect <- function(filename, mode = "r") {
471   if (!(mode %in% c('r', 'r+'))) {
472     stop("'mode' must be one of 'r' or 'r+'")
473   }
474   new.loom <- loom$new(filename = filename, mode = mode)
475   return(new.loom)
476 }
477
478 CreateLoomFromSeurat <- function(object, filename) {
479   object.data=t(object@raw.data[rownames(object@data),object@cell.names])
480   object.meta.data=object@meta.data
481   row_attrs=list(); col_attrs=list()
482   gene.names=colnames(object.data)
483   object.meta.data$ident = object@ident
484   object.meta.data$CellID = object@cell.names
485   for(i in 1:ncol(object.meta.data)) {
486     col_attrs[[colnames(object.meta.data)[i]]]=object.meta.data[,i]
487   }
488   row_attrs[["Gene"]]=gene.names
489   create(filename,object.data,gene.attrs = row_attrs, cell.attrs = col_attrs)
490 }
491
492 #need to comment
493 #need to add progress bar
494 #but otherwise, pretty cool
495 #for paul to try :
496 # f <- connect("~/Downloads/10X43_1.loom")
497 # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
498 # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
499 map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
500   n_func = length(f_list)
501   if (n_func == 1) {
502     f_list <- list(f_list)
503   }
504   if (MARGIN == 1) {
505     results <- list()
506     for (j in 1:n_func) {
507       results[[j]] <- numeric(0)
508     }
509     rows_per_chunk <- chunksize
510     ix <- 1
511     while (ix <= self@shape[1]) {
512       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
513       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
514       for (j in 1:n_func) {
515         new_results <- apply(chunk, 1, FUN = f_list[[j]])
516         results[[j]] <- c(results[[j]], new_results)
517       }
518       ix <- ix + chunksize
519     }
520   }
521   if (MARGIN == 2) {
522     results <- list()
523     for (j in 1:n_func) {
524       results[[j]] <- numeric(0)
525     }
526     cols_per_chunk <- chunksize
527     ix <- 1
528     while (ix <= self@shape[2]) {
529       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
530       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
531       for (j in 1:n_func) {
532         new_results <- apply(chunk, 2, FUN = f_list[[j]])
533         results[[j]] <- c(results[[j]], new_results)
534       }
535       ix <- ix + chunksize
536     }
537   }
538   if (n_func == 1) {
539     results <- results[[1]]
540   }
541   return(results)
542 }