Finish layer adding
[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)}}{Add extra information to this loom file; \code{attribute} is a named list where each element is a vector that is as long as one dimension of \code{/matrix}, \code{MARGIN} is either 1 for cells or 2 for genes}
26 #'   \item{\code{add.row.attribute(attribute)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 2)}}
27 #'   \item{\code{add.col.attribute(attribute)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 1)}}
28 #'   \item{\code{add.meta.data(meta.data)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 1)}}
29 #' }
30 #'
31 #' @importFrom utils packageVersion
32 #'
33 #' @export
34 #'
35 loom <- R6Class(
36   # The loom class
37   # Based on the H5File class from hdf5r
38   # Not clonable (no loom$clone method), doesn't make sense since all data is on disk, not in memory
39   # Yes to portability, other packages can subclass the loom class
40   # Class is locked, other fields and methods cannot be added
41   classname = 'loom',
42   inherit = hdf5r::H5File,
43   cloneable = FALSE,
44   portable = TRUE,
45   lock_class = TRUE,
46   # Public fields and methods
47   # See above for documentation
48   public = list(
49     # Fields
50     version = NULL,
51     shape = NULL,
52     chunksize = NULL,
53     matrix = NULL,
54     layers = NULL,
55     col.attrs = NULL,
56     row.attrs = NULL,
57     # Methods
58     initialize = function(filename = NULL, mode = c('a', 'r', 'r+', 'w', 'w-'), ...) {
59       # If the file exists, run validation steps
60       do.validate <- file.exists(filename) && !(mode %in% c('w', 'w+'))
61       super$initialize(filename = filename, mode = mode, ...)
62       if (do.validate) {
63         # Run the validation steps
64         validateLoom(object = self)
65         # Store /matrix and the shape of /matrix
66         self$matrix <- self[['matrix']]
67         self$shape <- self[['matrix']]$dims
68         # Store the chunk size
69         chunks <- h5attr(x = self, which = 'chunks')
70         chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
71         chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
72         chunks <- unlist(x = strsplit(x = chunks, split = ','))
73         self$chunksize <- as.integer(x = chunks)
74         # Store version information
75         self$version <- as.character(x = tryCatch(
76           # Try getting a version
77           # If it doesn't exist, can we write to the file?
78           # If so, store the version as this version of loomR
79           # Otherwise, store the version as NA_character_
80           expr = h5attr(x = self, which = 'version'),
81           error = function(e) {
82             if (mode != 'r') {
83               version <- packageVersion(pkg = 'loomR')
84               h5attr(x = self, which = 'version') <- version
85             } else {
86               version <- NA_character_
87             }
88             return(version)
89           }
90         ))
91         # Load layers
92         private$load_layers()
93         # Load attributes
94         private$load_attributes(MARGIN = 1) # Cells (col_attrs)
95         private$load_attributes(MARGIN = 2) # Genes (row_attrs)
96       } else {
97         # Assume new HDF5 file
98         self$version <- as.character(x = packageVersion(pkg = 'loomR'))
99       }
100     },
101     add.layer = function(layers) {
102       # Value checking
103       if (!is.list(x = layers) || is.null(x = names(x = layers))) {
104         stop("'layers' must be a named list")
105       }
106       if (is.null(x = self$shape)) {
107         stop(private$err_msg)
108       }
109       # Add layers
110       for (i in 1:length(x = layers)) {
111         if (!is.matrix(x = layers[[i]])) {
112           layers[[i]] <- as.matrix(x = layers[[i]])
113         }
114         do.transpose <- FALSE
115         if (any(dim(x = layers[[i]]) != self$shape)) {
116           if (all(rev(x = dim(x = layers[[i]])) == self$shape)) {
117             do.transpose <- TRUE
118           } else {
119             stop(paste(
120               "All layers must have",
121               self$shape[1],
122               "rows for cells and",
123               self$shape[2],
124               "columns for genes"
125             ))
126           }
127         }
128         if (do.transpose) {
129           self[['layers', names(x = layers)[i]]] <- t(x = layers[[i]])
130         } else {
131           self[['layers', names(x = layers)[i]]] <- layers[[i]]
132         }
133       }
134       self$flush()
135       private$load_layers()
136       invisible(x = self)
137     },
138     add.attribute = function(attribute, MARGIN) {
139       # Value checking
140       if (!is.list(x = attribute) || is.null(x = names(x = attribute))) {
141         stop("'attribute' must be a named list")
142       }
143       if (length(x = attribute) > 1) {
144         for (i in attribute) {
145           if (!is.vector(x = attribute)) {
146             stop("All attributes must be one-dimensional vectors")
147           }
148         }
149       }
150       if (length(x = which(x = names(x = attribute) != '')) != length(x = attribute)) {
151         stop("Not all attributes had names provided")
152       }
153       if (!MARGIN %in% c(1, 2)) {
154         stop("'MARGIN' must be 1 or 2")
155       }
156       # Add the attributes as datasets for our MARGIN's group
157       if (is.null(x = self$shape)) {
158         stop(private$err_msg)
159       }
160       grp.name <- c('col_attrs', 'row_attrs')[MARGIN]
161       grp <- self[[grp.name]]
162       for (i in 1:length(x = attribute)) {
163         if (length(attribute[[i]]) != self$shape[MARGIN])
164           stop(paste(
165             "All",
166             switch(EXPR = MARGIN, '1' = 'cell', '2' = 'gene'),
167             "attributes must be of length",
168             self$shape[MARGIN]
169           ))
170         grp[[names(x = attribute)[i]]] <- attribute[[i]]
171       }
172       self$flush()
173       gc(verbose = FALSE)
174       # Load the attributes for this margin
175       private$load_attributes(MARGIN = MARGIN)
176       invisible(x = self)
177     },
178     # Add attributes for genes
179     add.row.attribute = function(attribute) {
180       self$add.attribute(attribute = attribute, MARGIN = 2)
181       invisible(x = self)
182     },
183     # Add attributes for cells
184     add.col.attribute = function(attribute) {
185       self$add.attribute(attribute = attribute, MARGIN = 1)
186       invisible(x = self)
187     },
188     # Add metadata, follows cells
189     add.meta.data = function(meta.data) {
190       self$add.col.attribute(attribute = meta.data)
191       invisible(x = self)
192     }
193   ),
194   # Private fields and methods
195   # @field err_msg A simple error message if this object hasn't been created with loomR::create or loomR::connect
196   # \describe{
197   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
198   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
199   # }
200   private = list(
201     # Fields
202     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",
203     # Methods
204     load_attributes = function(MARGIN) {
205       attribute <- switch(
206         EXPR = MARGIN,
207         '1' = 'col_attrs',
208         '2' = 'row_attrs',
209         stop('Invalid attribute dimension')
210       )
211       group <- self[[attribute]]
212       attributes <- unlist(x = lapply(
213         X = names(x = group),
214         FUN = function(x) {
215           d <- list(group[[x]])
216           names(x = d) <- x
217           return(d)
218         }
219       ))
220       if (MARGIN == 1) {
221         self$col.attrs <- attributes
222       } else if (MARGIN == 2) {
223         self$row.attrs <- attributes
224       }
225     },
226     load_layers = function() {
227       self$layers <- unlist(x = lapply(
228         X = names(x = self[['layers']]),
229         FUN = function(n) {
230           d <- c(self[['layers', n]])
231           names(x = d) <- n
232           return(d)
233         }
234       ))
235     }
236   )
237 )
238
239 #' Create a loom object
240 #'
241 #' @param filename The name of the new loom file
242 #' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
243 #' @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}
244 #' @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}
245 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
246 #'
247 #' @return A connection to a loom file
248 #'
249 #' @importFrom utils packageVersion
250 #'
251 #' @seealso \code{\link{loom-class}}
252 #'
253 #' @export
254 #'
255 create <- function(
256   filename,
257   data,
258   gene.attrs = NULL,
259   cell.attrs = NULL,
260   layers = NULL,
261   chunk.dims = 'auto'
262 ) {
263   if (file.exists(filename)) {
264     stop(paste('File', file, 'already exists!'))
265   }
266   if (!is.matrix(x = data)) {
267     data <- as.matrix(x = data)
268   }
269   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
270     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
271   } else if (length(x = chunk.dims == 1)) {
272     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
273       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
274     }
275   } else {
276     chunk.dims <- as.integer(x = chunk.dims)
277   }
278   new.loom <- loom$new(filename = filename, mode = 'w-')
279   # Create the matrix
280   new.loom$create_dataset(
281     name = 'matrix',
282     robj = t(x = data),
283     chunk_dims = chunk.dims
284   )
285   new.loom$shape <- new.loom[['matrix']]
286   new.loom$shape <- new.loom[['matrix']]$dims
287   if (!is.null(x = colnames(x = data))) {
288     new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
289   }
290   if (!is.null(x = rownames(x = data))) {
291     new.loom$add.col.attribute(attribute = list('cell_names' = colnames(x = data)))
292   }
293   # Store some constants as HDF5 attributes
294   h5attr(x = new.loom, which = 'version') <- new.loom$version
295   h5attr(x = new.loom, which = 'chunks') <- paste0(
296     '(',
297     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
298     ')'
299   )
300   # Groups
301   new.loom$create_group(name = 'layers')
302   new.loom$create_group(name = 'row_attrs')
303   new.loom$create_group(name = 'col_attrs')
304   # Add layers
305   for (ly in layers) {
306     new.loom$add.layer(layer = ly)
307   }
308   if (!is.null(x = gene.attrs)) {
309     new.loom$add.row.attribute(attribute = gene.attrs)
310   }
311   if (!is.null(x = cell.attrs)) {
312     new.loom$add.col.attribute(attribute = cell.attrs)
313   }
314   # Set last bit of information
315   chunks <- new.loom[['matrix']]$chunk_dims
316   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
317   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
318   chunks <- unlist(x = strsplit(x = chunks, split = ','))
319   new.loom$chunksize <- as.integer(x = chunks)
320   # Return the connection
321   return(new.loom)
322 }
323
324 #' Validate a loom object
325 #'
326 #' @param object A loom object
327 #'
328 #' @return None, errors out if object is an invalid loom connection
329 #'
330 #' @seealso \code{\link{loom-class}}
331 #'
332 #' @export
333 #'
334 validateLoom <- function(object) {
335   # A loom file is a specific HDF5
336   # We need a dataset in /matrix that's a two-dimensional dense matrix
337   root.datasets <- list.datasets(object = object, path = '/', recursive = FALSE)
338   if (length(x = root.datasets) != 1) {
339     stop("There can only be one dataset at the root of the loom file")
340   }
341   if (root.datasets != 'matrix') {
342     stop("The root dataset must be called '/matrix'")
343   }
344   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
345   required.groups <- c('row_attrs', 'col_attrs', 'layers')
346   dim.matrix <- object[[root.datasets]]$dims # Columns x Rows
347   names(x = dim.matrix) <- required.groups[c(2, 1)]
348   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
349   group.msg <- paste0(
350     "There can only be three groups in the loom file: '",
351     paste(required.groups, collapse = "', '"),
352     "'"
353   )
354   if (length(x = root.groups) != 3) {
355     stop(group.msg)
356   }
357   if (!all(required.groups %in% root.groups)) {
358     stop(group.msg)
359   }
360   unlist(x = sapply(
361     X = required.groups[1:2],
362     FUN = function(group) {
363       if (length(x = list.groups(object = object[[group]], recursive = FALSE)) > 0) {
364         stop(paste("Group", group, "cannot have subgroups"))
365       }
366       if (length(x = list.attributes(object = object[[group]])) > 0) {
367         stop(paste("Group", group, "cannot have subattributes"))
368       }
369       for (dataset in list.datasets(object = object[[group]])) {
370         if (object[[paste(group, dataset, sep = '/')]]$dims != dim.matrix[group]) {
371           stop(paste("All datasets in group", group, "must be of length", required.groups[group]))
372         }
373       }
374     }
375   ))
376   for (dataset in list.datasets(object = object[['/layers']])) {
377     if (any(object[[paste('layers', dataset, sep = '/')]]$dims != dim.matrix)) {
378       stop(paste("All datasets in '/layers' must be", dim.matrix[1], 'by', dim.matrix[2]))
379     }
380   }
381 }
382
383 #' Connect to a loom file
384 #'
385 #' @param filename The loom file to connect to
386 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write
387 #'
388 #' @return A loom file connection
389 #'
390 #' @seealso \code{\link{loom-class}}
391 #'
392 #' @export
393 #'
394 connect <- function(filename, mode = "r") {
395   if (!(mode %in% c('r', 'r+'))) {
396     stop("'mode' must be one of 'r' or 'r+'")
397   }
398   new.loom <- loom$new(filename = filename, mode = mode)
399   return(new.loom)
400 }
401
402 #need to comment
403 #need to add progress bar
404 #but otherwise, pretty cool
405 #for paul to try :
406 # f <- connect("~/Downloads/10X43_1.loom")
407 # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
408 # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
409 map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
410   n_func = length(f_list)
411   if (n_func == 1) {
412     f_list <- list(f_list)
413   }
414   if (MARGIN == 1) {
415     results <- list()
416     for (j in 1:n_func) {
417       results[[j]] <- numeric(0)
418     }
419     rows_per_chunk <- chunksize
420     ix <- 1
421     while (ix <= self@shape[1]) {
422       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
423       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
424       for (j in 1:n_func) {
425         new_results <- apply(chunk, 1, FUN = f_list[[j]])
426         results[[j]] <- c(results[[j]], new_results)
427       }
428       ix <- ix + chunksize
429     }
430   }
431   if (MARGIN == 2) {
432     results <- list()
433     for (j in 1:n_func) {
434       results[[j]] <- numeric(0)
435     }
436     cols_per_chunk <- chunksize
437     ix <- 1
438     while (ix <= self@shape[2]) {
439       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
440       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
441       for (j in 1:n_func) {
442         new_results <- apply(chunk, 2, FUN = f_list[[j]])
443         results[[j]] <- c(results[[j]], new_results)
444       }
445       ix <- ix + chunksize
446     }
447   }
448   if (n_func == 1) {
449     results <- results[[1]]
450   }
451   return(results)
452 }