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