Merge branch 'develop' of github.com:mojaveazure/loomR into develop
[loomr.git] / R / loom.R
1 #' @include internal.R
2 #' @import hdf5r
3 #' @importFrom R6 R6Class
4 NULL
5
6 #' A class for loom files
7 #'
8 #' @docType class
9 #' @name loom-class
10 #' @rdname loom-class
11 #' @aliases loom
12 #' @return Object of \code{\link{R6::R6Class}} to generate \code{loom} objects
13 #' @format An \code{\link{R6::R6Class}} object
14 #' @seealso \code{\link{hdf5r::H5File}}
15 #'
16 #' @field version Version of loomR object was created under
17 #' @field shape Shape of \code{/matrix} in genes (columns) by cells (rows)
18 #' @field chunksize Chunks set for this dataset in columns (cells) by rows (genes)
19 #' @field matrix The main data matrix, stored as columns (cells) by rows (genes)
20 #' @field layers Additional data matricies, the same shape as \code{/matrix}
21 #' @field col.attrs Extra information about cells
22 #' @field row.attrs Extra information about genes
23 #'
24 #' @section Methods:
25 #' \describe{
26 #'   \item{\code{add.layer(layer, chunk.size, overwrite)}}{
27 #'     Add a data layer to this loom file, must be the same dimensions as \code{/matrix}
28 #'     \describe{
29 #'       \item{\code{layer}}{A named list of matrices to be added as layers}
30 #'       \item{\code{chunk.size}}{Number of rows from each layer to stream at once, defaults to 1000}
31 #'       \item{\code{overwrite}}{If a layer already exists, overwrite with new data, defaults to FALSE}
32 #'     }
33 #'   }
34 #'   \item{\code{add.attribute(attribute, MARGIN, overwrite)}}{
35 #'     Add extra information to this loom file.
36 #'     \describe{
37 #'       \item{\code{attribute}}{A named list where the first dimmension of each element as long as one dimension of \code{/matrix}}
38 #'       \item{\code{MARGIN}}{Either 1 for genes or 2 for cells}
39 #'       \item{\code{overwrite}}{Can overwrite existing attributes?}
40 #'     }
41 #'   }
42 #'   \item{\code{add.row.attribute(attribute), add.col.attribute(attribute)}}{
43 #'     Add row or column attributes
44 #'   }
45 #'   \item{\code{batch.scan(chunk.size, MARGIN, index.use, dataset.use, force.reset)}, \code{batch.next(return.data)}}{
46 #'     Scan a dataset in the loom file from \code{index.use[1]} to \code{index.use[2]}, iterating by \code{chunk.size}.
47 #'     \code{dataset.use} can be the name, not \code{group/name}, unless the name is present in multiple groups.
48 #'     Pass \code{MARGIN = 1} to iterate on genes or \code{MARGIN = 2} to iterate on cells for 'matrix' or any dataset in 'layers'.
49 #'     To force reevaluation of the iterator object, pass \code{force.reset = TRUE}.
50 #'     \code{MARGIN} does not need to be set for datasets in 'row_attrs' or 'col_attrs'.
51 #'     \code{chunk.size} defaults to \code{self$chunksize}, \code{MARGIN} defaults to 2,
52 #'     \code{index.use} defaults to \code{1:self$shape[MARGIN]}, \code{dataset.use} defaults to 'matrix'
53 #'   }
54 #'   \item{\code{apply(name, FUN, MARGIN, chunk.size, dataset.use, overwrite, display.progress, ...)}}{
55 #'     Apply a function over a dataset within the loom file, stores the results in the loom file.
56 #'     \code{name} must be the full name of the dataset ('group/name').
57 #'     \code{apply} will always use the entire dataset specified in \code{dataset.use}
58 #'     \describe{
59 #'       \item{\code{name}}{Name of dataset to store results to}
60 #'       \item{\code{FUN}}{Function to apply}
61 #'       \item{\code{MARGIN}}{Iterate over genes (1) or cells (2)}
62 #'       \item{\code{chunk.size}}{Size to chunk \code{MARGIN} by}
63 #'       \item{\code{dataset.use}}{Name of dataset to use}
64 #'       \item{\code{overwrite}}{Overite \code{name} if already exists}
65 #'       \item{\code{display.progress}}{Display progress}
66 #'       \item{\code{...}}{Extra parameters to pass to \code{FUN}}
67 #'     }
68 #'   }
69 #'   \item{\code{map(FUN, MARGIN, chunk.size, index.use, dataset.use, display.progress, expected, ...)}}{
70 #'     Map a function onto a dataset within the loom file, returns the result into R.
71 #'     The result will default to the shape of the dataset used; to change pass either 'vector' or 'matrix' to \code{expected}.
72 #'   }
73 #'   \item{\code{add.cells(matrix.data, attributes.data = NULL, layers.data = NULL, display.progress = TRUE)}}{
74 #'     Add m2 cells to a loom file.
75 #'     \describe{
76 #'       \item{\code{matrix.data}}{a list of m2 cells where each entry is a vector of length n (num genes, \code{self$shape[1]})}
77 #'       \item{\code{attributes.data}}{a list where each entry is named for one of the datasets in \code{self[['col_attrs']]}; each entry is a vector of length m2.}
78 #'       \item{\code{layers.data}}{a list where each entry is named for one of the datasets in \code{self[['layers']]}; each entry is an n-by-m2 matrix where n is the number of genes in this loom file and m2 is the number of cells being added.}
79 #'       \item{\code{display.progress}}{display progress}
80 #'     }
81 #'   }
82 #' }
83 #'
84 #' @importFrom iterators nextElem
85 #' @importFrom itertools hasNext ihasNext ichunk
86 #' @importFrom utils packageVersion txtProgressBar setTxtProgressBar
87 #'
88 #' @export
89 #'
90 loom <- R6Class(
91   # The loom class
92   # Based on the H5File class from hdf5r
93   # Not clonable (no loom$clone method), doesn't make sense since all data is on disk, not in memory
94   # Yes to portability, other packages can subclass the loom class
95   # Class is locked, other fields and methods cannot be added
96   classname = 'loom',
97   inherit = hdf5r::H5File,
98   cloneable = FALSE,
99   portable = TRUE,
100   lock_class = TRUE,
101   # Public fields and methods
102   # See above for documentation
103   public = list(
104     # Fields
105     version = NULL,
106     shape = NULL,
107     chunksize = NULL,
108     matrix = NULL,
109     layers = NULL,
110     col.attrs = NULL,
111     row.attrs = NULL,
112     # Methods
113     initialize = function(
114       filename = NULL,
115       mode = c('a', 'r', 'r+', 'w', 'w-'),
116       skip.validate = FALSE,
117       ...
118     ) {
119       # If the file exists, run validation steps
120       do.validate <- file.exists(filename) && !(mode %in% c('w', 'w+'))
121       private$skipped.validation <- skip.validate
122       super$initialize(filename = filename, mode = mode, ...)
123       if (do.validate) {
124         # Run the validation steps
125         if (skip.validate) {
126           warning("Skipping validation step, some fields are not populated")
127         } else {
128           validateLoom(object = self)
129         }
130         # Store /matrix and the shape of /matrix
131         if (skip.validate) {
132           if (getOption(x = 'verbose')) {
133             warning("Not setting matrix field")
134           }
135         } else {
136           self$matrix <- self[['matrix']]
137         }
138         self$shape <- rev(self[['matrix']]$dims)
139         # Store the chunk size
140         chunks <- tryCatch(
141           expr = h5attr(x = self, which = 'chunks'),
142           error = function(e) {
143             hchunks <- self[['matrix']]$chunk_dims
144             pchunks <- paste0('(', paste(hchunks, collapse = ', '), ')')
145             if (mode != 'r') {
146               h5attr(x = self, which = 'chunks') <- pchunks
147             }
148             return(pchunks)
149           }
150         )
151         chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
152         chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
153         chunks <- unlist(x = strsplit(x = chunks, split = ','))
154         self$chunksize <- as.integer(x = chunks)
155         # Store version information
156         self$version <- as.character(x = tryCatch(
157           # Try getting a version
158           # If it doesn't exist, can we write to the file?
159           # If so, store the version as this version of loomR
160           # Otherwise, store the version as NA_character_
161           expr = h5attr(x = self, which = 'version'),
162           error = function(e) {
163             if (mode != 'r') {
164               version <- packageVersion(pkg = 'loomR')
165               h5attr(x = self, which = 'version') <- as.character(x = version)
166             } else {
167               version <- NA_character_
168             }
169             return(version)
170           }
171         ))
172         # Load layers
173         private$load_layers()
174         # Load attributes
175         private$load_attributes(MARGIN = 1) # Genes (row_attrs)
176         private$load_attributes(MARGIN = 2) # Cells (col_attrs
177       } else {
178         # Assume new HDF5 file
179         self$version <- as.character(x = packageVersion(pkg = 'loomR'))
180       }
181     },
182     finalizer = function() {
183       self$close_all(close_self = TRUE)
184     },
185     load.fields = function() {
186       private$load_layers()
187       private$load_attributes(MARGIN = 1)
188       private$load_attributes(MARGIN = 2)
189     },
190     # Addding attributes and layers
191     add.layer = function(layers, chunk.size = 1000, overwrite = FALSE) {
192       if (self$mode == 'r') {
193         stop(private$err_mode)
194       }
195       # Value checking
196       if (!is.list(x = layers) || is.null(x = names(x = layers))) {
197         stop("'layers' must be a named list")
198       }
199       if (is.null(x = self$shape)) {
200         stop(private$err_init)
201       }
202       # Add layers
203       for (i in 1:length(x = layers)) {
204         # if (!is.matrix(x = layers[[i]])) {
205         if (!inherits(x = layers[[i]], what = c('Matrix', 'matrix'))) {
206           layers[[i]] <- as.matrix(x = layers[[i]])
207         }
208         do.transpose <- FALSE
209         shape.use <- rev(x = self$shape)
210         if (any(dim(x = layers[[i]]) != shape.use)) {
211           if (all(rev(x = dim(x = layers[[i]])) == shape.use)) {
212             do.transpose <- TRUE
213           } else {
214             stop(paste(
215               "All layers must have",
216               shape.use[1],
217               "rows for cells and",
218               shape.use[2],
219               "columns for genes"
220             ))
221           }
222         }
223         if (do.transpose) {
224           layers[[i]] <- t(x = layers[[i]])
225         }
226         layer.name <- names(x = layers)[i]
227         if (layer.name %in% list.datasets(object = self[['layers']])) {
228           if (overwrite) {
229             self[['layers']]$link_delete(name = layer.name)
230           } else {
231             stop(paste(
232               "A layer with the name",
233               layer.name,
234               "already!"
235             ))
236           }
237         }
238         dtype <- getDtype(x = layers[[i]][1, 1])
239         self[['layers']]$create_dataset(
240           name = layer.name,
241           dtype = dtype,
242           dims = dim(x = layers[[i]])
243         )
244         chunk.points <- chunkPoints(
245           data.size = dim(x = layers[[i]])[1],
246           chunk.size = chunk.size
247         )
248         # if (display.progress) {
249         #   pb <- txtProgressBar(char = '=', style = 3)
250         # }
251         for (col in 1:ncol(x = chunk.points)) {
252           row.start <- chunk.points[1, col]
253           row.end <- chunk.points[2, col]
254           self[['layers']][[layer.name]][row.start:row.end, ] <- as.matrix(
255             x = layers[[i]][row.start:row.end, ]
256           )
257           # if (display.progress) {
258           #   setTxtProgressBar(pb = pb, value = col / ncol(x = chunk.points))
259           # }
260         }
261         # self[['layers']]$create_dataset(
262         #   name = names(x = layers)[i],
263         #   robj = layers[[i]],
264         #   chunk_dims = self$chunksize
265         # )
266       }
267       self$flush()
268       gc(verbose = FALSE)
269       private$load_layers()
270       invisible(x = self)
271     },
272     add.attribute = function(attribute, MARGIN, overwrite = FALSE) {
273       if (self$mode == 'r') {
274         stop(private$err_mode)
275       }
276       # Value checking
277       if (is.data.frame(x = attribute)) {
278         attribute <- as.list(x = attribute)
279       }
280       is.actual.list <- is.list(x = attribute)
281       if (!is.actual.list || is.null(x = names(x = attribute))) {
282         stop("Attributes must be provided as a named list")
283       }
284       # if (is.data.frame(x = attribute)) {
285       #   attribute <- as.list(x = attribute)
286       # }
287       # if (!is.list(x = attribute) || is.null(x = names(x = attribute))) {
288       #   stop("'attribute' must be a named list")
289       # }
290       if (!MARGIN %in% c(1, 2)) {
291         stop("'MARGIN' must be 1 or 2")
292       }
293       length.use <- rev(x = self[['matrix']]$dims)[MARGIN]
294       dim.msg <- paste(
295         "At least one dimmension for each",
296         switch(EXPR = MARGIN, '1' = 'gene', '2' = 'cell'),
297         "attribute must be",
298         length.use
299       )
300       for (i in 1:length(x = attribute)) {
301         if (is.matrix(x = attribute[[i]]) || is.data.frame(x = attribute[[i]])) {
302           margin.use <- which(x = dim(x = attribute[[i]]) == length.use)
303           if (!length(x = margin.use)) {
304             stop(dim.msg)
305           }
306           margin.use <- margin.use[1]
307           attribute[[i]] <- switch(
308             EXPR = margin.use,
309             '1' = t(x = as.matrix(x = attribute[[i]])),
310             '2' = as.matrix(x = attribute[[i]]),
311             stop("All attributes must be one- or two-dimmensional")
312           )
313         } else {
314           if (length(x = attribute[[i]]) != length.use) {
315             stop(dim.msg)
316           }
317           attribute[[i]] <- as.vector(x = attribute[[i]])
318         }
319       }
320       if (length(x = which(x = names(x = attribute) != '')) != length(x = attribute)) {
321         stop("Not all attributes had names provided")
322       }
323       grp.name <- c('row_attrs', 'col_attrs')[MARGIN]
324       grp <- self[[grp.name]]
325       if (!overwrite) {
326         names.fail <- which(x = names(x = attribute) %in% names(x = grp))
327         if (length(x = names.fail) != 0) {
328           stop(paste0(
329             "The following names are already used for ",
330             switch(EXPR = MARGIN, '1' = 'row', '2' = 'column'),
331             " attributes: '",
332             paste(names(x = attribute)[names.fail], collapse = ''),
333             "'"
334           ))
335         }
336       }
337       # Add the attributes as datasets for our MARGIN's group
338       for (i in 1:length(x = attribute)) {
339         try(expr = grp$link_delete(name = names(x = attribute)[i]), silent = TRUE)
340         grp[[names(x = attribute)[i]]] <- attribute[[i]]
341       }
342       self$flush()
343       gc(verbose = FALSE)
344       # Load the attributes for this margin
345       private$load_attributes(MARGIN = MARGIN)
346       invisible(x = self)
347     },
348     add.row.attribute = function(attribute, overwrite = FALSE) {
349       self$add.attribute(attribute = attribute, MARGIN = 1, overwrite = overwrite)
350       invisible(x = self)
351     },
352     add.col.attribute = function(attribute, overwrite = FALSE) {
353       self$add.attribute(attribute = attribute, MARGIN = 2, overwrite = overwrite)
354       invisible(x = self)
355     },
356     add.meta.data = function(meta.data, overwrite = FALSE) {
357       self$add.col.attribute(attribute = meta.data, overwrite = overwrite)
358       invisible(x = self)
359     },
360     # Get attribute information
361     get.attribute.df = function(
362       attribute.layer = c("row", "col"),
363       attribute.names = NULL,
364       row.names = "gene_names",
365       col.names = "cell_names"
366     ) {
367       # takes multiple row_attrs or col_attrs and combines them into a data.frame,
368       # removing rows or columns that are entirely NAs.
369       #
370       # attribute.layer either "row" for row_attrs or "col" col_attrs
371       # attribute.names name of rows or columns to combine into matrix
372       # row.names either a character vector or name of an element in row_attrs
373       # col.names either a character vector or name of an element in col_attrs
374       if (!attribute.layer %in% c("row", "col")) {
375         stop("Invalid attribute.layer. Please select either 'row' or 'col'.")
376       }
377       attribute.layer <- paste0(attribute.layer, "_attrs")
378       # check that attribute names are present
379       if (!all(attribute.names %in% self[[attribute.layer]]$names)) {
380         invalid.names <- attribute.names[which(!attribute.names %in% self[[attribute.layer]]$names)]
381         stop(paste0("Invalid attribute.names: ", paste0(invalid.names, collapse = ", ")))
382       }
383       if (attribute.layer == "row_attrs") {
384         combined.df <- data.frame(
385           self[[paste0(attribute.layer, "/", attribute.names[1])]][],
386           row.names = self[[paste0(attribute.layer, "/", row.names)]][]
387         )
388       } else {
389         combined.df <- data.frame(
390           self[[paste0(attribute.layer, "/", attribute.names[1])]][],
391           row.names = self[[paste0(attribute.layer, "/", col.names)]][]
392         )
393       }
394       if (length(x = attribute.names) > 1) {
395         for (i in 2:length(x = attribute.names)) {
396           combined.df[, attribute.names[i]] <- self[[paste0(attribute.layer, "/", attribute.names[i])]][]
397         }
398       }
399       colnames(x = combined.df) <- attribute.names
400       # check if any row is all NAs
401       rows.to.remove <- unname(obj = which(x = apply(
402         X = combined.df,
403         MARGIN = 1,
404         FUN = function(x) {
405           return(all(is.na(x = x)))
406         }
407       )))
408       if (length(x = rows.to.remove) > 1) {
409         combined.df <- combined.df[-rows.to.remove, ]
410       }
411       return(combined.df)
412     },
413     # Chunking functions
414     batch.scan = function(
415       chunk.size = NULL,
416       MARGIN = 2,
417       index.use = NULL,
418       dataset.use = 'matrix',
419       force.reset = FALSE
420     ) {
421       if (is.null(x = private$it) || !grepl(pattern = dataset.use, x = private$iter.dataset) || force.reset) {
422         # Check the existence of the dataset
423         private$iter.dataset <- grep(
424           pattern = dataset.use,
425           x = list.datasets(object = self),
426           value = TRUE
427         )
428         if (length(x = private$iter.dataset) != 1) {
429           private$reset_batch()
430           stop(paste0("Cannot find dataset '", dataset.use, "' in the loom file"))
431         }
432         if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
433           MARGIN <- 1
434         }
435         else if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
436           MARGIN <- 2
437         }
438         # Check the margin
439         if (!(MARGIN %in% c(1, 2))) {
440           private$reset_batch()
441           stop("MARGIN must be 1 (genes) or 2 (cells)")
442         } else {
443           private$iter.margin <- MARGIN
444         }
445         if (is.null(x = chunk.size)) {
446           chunk.size <- rev(x = self$chunksize)[private$iter.margin]
447         }
448         private$iter.chunksize <- chunk.size
449         # Set the indices to use
450         index.use <- private$iter_range(index.use = index.use)
451         # Setup our iterator
452         private$it <- ihasNext(iterable = ichunk(
453           iterable = index.use[1]:index.use[2],
454           chunkSize = chunk.size
455         ))
456         private$iter.index <- index.use
457       }
458       # Return the times we iterate, this isn't important, we only need the length of this vector
459       return(1:ceiling((private$iter.index[2] - private$iter.index[1]) / private$iter.chunksize))
460     },
461     batch.next = function(return.data = TRUE) {
462       # Ensure that we have a proper iterator
463       if (!'hasNext.ihasNext' %in% suppressWarnings(expr = methods(class = class(x = private$it)))) {
464         stop("Please setup the iterator with self$batch.scan")
465       }
466       # Do the iterating
467       if (hasNext(obj = private$it)) {
468         # Get the indices for this chunk
469         chunk.indices <- unlist(x = nextElem(obj = private$it))
470         if (return.data) {
471           # If we're working with a matrix dataset, ensure chunking on proper axis
472           if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
473             to.return <- switch(
474               EXPR = private$iter.margin,
475               '1' = self[[private$iter.dataset]][, chunk.indices], # Genes
476               '2' = self[[private$iter.dataset]][chunk.indices, ] # Cells
477             )
478           } else {
479             # Otherwise, iterating over an attribute (1 dimensional)
480             to.return <- self[[private$iter.dataset]][chunk.indices]
481           }
482         } else {
483           to.return <- chunk.indices
484         }
485         # Determine if we reset the iterator
486         if (!hasNext(obj = private$it)) {
487           private$reset_batch()
488         }
489         return(to.return)
490       } else {
491         # Just in case
492         private$reset_batch()
493         return(NULL)
494       }
495     },
496     apply = function(
497       name,
498       FUN,
499       MARGIN = 2,
500       index.use = NULL,
501       chunk.size = NULL,
502       dataset.use = 'matrix',
503       overwrite = FALSE,
504       display.progress = TRUE,
505       ...
506     ) {
507       if (self$mode == 'r') {
508         stop(private$err_mode)
509       }
510       if (!inherits(x = FUN, what = 'function')) {
511         stop("FUN must be a function")
512       }
513       # Check that we're storing our results properly
514       results.basename <- basename(path = name) # Basename of our results
515       results.dirname <- gsub(pattern = '/', replacement = '', x = dirname(path = name)) # Groupname of our results
516       dirnames <- c('row_attrs', 'col_attrs', 'layers') # Allowed group names
517       if (name %in% list.datasets(object = self)) {
518         if (overwrite) {
519           self$link_delete(name = name)
520         } else {
521           stop(paste("A dataset with the name", name, "already exists!"))
522         }
523       }
524       # Checks datset, index, and MARGIN
525       # Sets full dataset path in private$iter.dataset
526       # Sets proper MARGIN in private$iter.margin
527       batch <- self$batch.scan(
528         chunk.size = chunk.size,
529         MARGIN = MARGIN,
530         dataset.use = dataset.use,
531         force.reset = TRUE
532       )
533       MARGIN <- private$iter.margin
534       dataset.use <- private$iter.dataset
535       # Ensure that our group name is allowed
536       name.check <- which(x = dirnames == results.dirname)
537       if (!any(name.check)) {
538         private$reset_batch()
539         stop(paste(
540           "The dataset must go into one of:",
541           paste(dirnames, collapse = ', ')
542         ))
543       }
544       # Check that our group matches our iteration
545       # ie. To store in col_attrs, MARGIN must be 1
546       if (name.check %in% c(1, 2) && name.check != private$iter.margin) {
547         private$reset_batch()
548         stop(paste(
549           "Iteration must be over",
550           c('genes', 'cells')[name.check],
551           paste0("(MARGIN = ", name.check, ")"),
552           "to store results in",
553           paste0("'", dirnames[name.check], "'")
554         ))
555       }
556       # Check how we store our results
557       dataset.matrix <- ('matrix' %in% private$iter.dataset || grepl(pattern = 'layers', x = private$iter.dataset))
558       results.matrix <- name.check == 3
559       # Get a connection to the group we're iterating over
560       group <- self[[results.dirname]]
561       # Ensure index.use is integers within the bounds of [1, self$shape[MARGIN]]
562       if (!is.null(x = index.use)) {
563         # Filter index.use to values between 1 and self$shape[MARGIN]
564         index.use <- as.integer(x = index.use)
565         index.use[index.use >= 1 & index.use <= self$shape[MARGIN]]
566         index.use <- as.vector(x = na.omit(object = index.use))
567         # If we still have values, figure out NAs, otherwise set index.use to NULL
568         if (length(x = index.use) > 0) {
569           # Do a trial run to figure out the class of NAs
570           na.use <- NA
571           if (display.progress) {
572             catn("Running trial to determine class of NAs")
573           }
574           trial <- switch(
575             EXPR = MARGIN,
576             '1' = self[[dataset.use]][, 1],
577             '2' = self[[dataset.use]][1, ]
578           )
579           trial <- FUN(trial, ...)
580           if (is.list(x = trial)) {
581             trial <- unlist(x = trial)
582           }
583           class(x = na.use) <- class(x = trial)
584         } else {
585           warning("No values passed to 'index.use' fall within the data, using all values")
586           index.use <- NULL
587         }
588       }
589       if (display.progress) {
590         pb <- txtProgressBar(char = '=', style = 3)
591       }
592       # Have to initialize the dataset differently than
593       # appending to it
594       first <- TRUE
595       for (i in 1:length(x = batch)) {
596         # Get the indices we're iterating over
597         these.indices <- self$batch.next(return.data = FALSE)
598         if (is.null(x = index.use)) {
599           chunk.indices <- these.indices
600         } else {
601           chunk.indices <- index.use[index.use %in% these.indices]
602           chunk.na <- these.indices[!(these.indices %in% chunk.indices)]
603         }
604         # Get the data and apply FUN
605         chunk.data <- if (dataset.matrix) {
606           switch(
607             EXPR = MARGIN,
608             '1' = self[[dataset.use]][, chunk.indices], # Chunk genes
609             '2' = self[[dataset.use]][chunk.indices, ] # Chunk cells
610           )
611         } else {
612           self[[private$iter.datset]][chunk.indices]
613         }
614         chunk.data <- FUN(chunk.data, ...)
615         # If this is the first iteration
616         # Initialize the dataset within group, set first to FALSE
617         if (first) {
618           if (!is.null(x = index.use)) {
619             # If we had indices to chunk on, create a holding matrix the size
620             # of what the results should be, add the
621             holding <- switch(
622               EXPR = MARGIN,
623               '1' = matrix(nrow = nrow(x = chunk.data), ncol = length(x = these.indices)),
624               '2' = matrix(nrow = length(x = these.indices), ncol = ncol(x = chunk.data))
625             )
626             switch (
627               EXPR = MARGIN,
628               '1' = holding[, chunk.indices] <- chunk.data,
629               '2' = holding[chunk.indices, ] <- chunk.data
630             )
631             chunk.data <- holding
632           }
633           group[[results.basename]] <- chunk.data
634           first <- FALSE
635         } else {
636           # If we're writign to a matrix
637           # Figure out which way we're writing the data
638           if (results.matrix) {
639             switch(
640               EXPR = MARGIN,
641               '1' = group[[results.basename]][, chunk.indices] <- chunk.data,
642               '2' = group[[results.basename]][chunk.indices, ] <- chunk.data
643             )
644             if (!is.null(x = index.use)) {
645               switch(
646                 EXPR = MARGIN,
647                 '1' = group[results.basename][, chunk.na] <- na.use,
648                 '2' = group[results.basename][chunk.na, ] <- na.use,
649               )
650             }
651           } else {
652             # Just write to the vector
653             group[[results.basename]][chunk.indices] <- chunk.data
654             if (!is.null(x = index.use)) {
655               group[[results.basename]][chunk.na] <- na.use
656             }
657           }
658         }
659         if (display.progress) {
660           setTxtProgressBar(pb = pb, value = i / length(x = batch))
661         }
662       }
663       if (display.progress) {
664         close(con = pb)
665       }
666       # Clean up and allow chaining
667       private$reset_batch()
668       # Load layers and attributes
669       private$load_layers()
670       private$load_attributes(MARGIN = 1) # Genes (row_attrs)
671       private$load_attributes(MARGIN = 2) # Cells (col_attrs)
672       invisible(x = self)
673     },
674     map = function(
675       FUN,
676       MARGIN = 2,
677       chunk.size = NULL,
678       index.use = NULL,
679       dataset.use = 'matrix',
680       display.progress = TRUE,
681       expected = NULL,
682       ...
683     ) {
684       if (!inherits(x = FUN, what = 'function')) {
685         stop("FUN must be a function")
686       }
687       # Checks datset, index, and MARGIN
688       # Sets full dataset path in private$iter.dataset
689       # Sets proper MARGIN in private$iter.margin
690       batch <- self$batch.scan(
691         chunk.size = chunk.size,
692         MARGIN = MARGIN,
693         index.use = index.use,
694         dataset.use = dataset.use,
695         force.reset = TRUE
696       )
697       MARGIN <- private$iter.margin
698       dataset.use <- private$iter.dataset
699       # Check how we store our results
700       # And what the shape of our dataset is
701       results.matrix <- TRUE
702       dataset.matrix <- TRUE
703       if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
704         results.matrix <- FALSE
705         dataset.matrix <- FALSE
706       } else if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
707         results.matrix <- FALSE
708         dataset.matrix <- FALSE
709       }
710       if (!is.null(x = expected)) {
711         results.matrix <- switch(
712           EXPR = expected,
713           'vector' = FALSE,
714           'matrix' = TRUE,
715           stop("'expected' must be one of 'matrix', 'vector', or NULL")
716         )
717       }
718       # Create our results holding object
719       results <- vector(mode = "list", length = length(x = batch))
720       if (display.progress) {
721         pb <- txtProgressBar(char = '=', style = 3)
722       }
723       for (i in 1:length(x = batch)) {
724         chunk.indices <- self$batch.next(return.data = FALSE)
725         chunk.data <- if (dataset.matrix) {
726           switch(
727             EXPR = MARGIN,
728             '1' = self[[dataset.use]][, chunk.indices], # Chunk genes
729             '2' = self[[dataset.use]][chunk.indices, ] # Chunk cells
730           )
731         } else {
732           self[[dataset.use]][chunk.indices]
733         }
734         results[[i]] <- FUN(chunk.data, ...)
735         if (display.progress) {
736           setTxtProgressBar(pb = pb, value = i / length(x = batch))
737         }
738       }
739       # Bring result list into vector or matrix format
740       if (class(results[[1]]) == 'numeric') {
741         results <- unlist(x = results, use.names = FALSE)
742       } else {
743         if (MARGIN == 1) {
744           results <- Reduce(f = cbind, x = results)
745         } else if (MARGIN == 2) {
746           results <- Reduce(f = rbind, x = results)
747         }
748       }
749       if (display.progress) {
750         close(con = pb)
751       }
752       private$reset_batch()
753       return(results)
754     },
755     # Functions that modify `/matrix'
756     add.cells = function(
757       matrix.data,
758       attributes.data = NULL,
759       layers.data = NULL,
760       display.progress = TRUE
761     ) {
762       if (self$mode == 'r') {
763         stop(private$err_mode)
764       }
765       # Check inputs
766       n <- self[['matrix']]$dims[2]
767       if (display.progress) {
768         cate("Checking inputs...")
769       }
770       matrix.data <- check.matrix_data(x = matrix.data, n = n)
771       layers.data <- check.layers(
772         x = layers.data,
773         n = n,
774         layers.names = names(x = self[['layers']])
775       )
776       attributes.data <- check.col_attrs(
777         x = attributes.data,
778         attrs.names = names(x = self[['col_attrs']])
779       )
780       # Get the number of cells we're adding
781       num.cells <- c(
782         nCells.matrix_data(x = matrix.data),
783         nCells.layers(x = layers.data),
784         nCells.col_attrs(x = attributes.data)
785       )
786       num.cells <- max(num.cells)
787       # Flesh out the input data to add
788       if (display.progress) {
789         cate(paste("Adding", num.cells, "to this loom file"))
790       }
791       matrix.data <- addCells.matrix_data(x = matrix.data, n = n, m2 = num.cells)
792       layers.data <- addCells.layers(x = layers.data, n = n, m2 = num.cells)
793       attributes.data <- addCells.col_attrs(x = attributes.data, m2 = num.cells)
794       # Add the input to the loom file
795       dims.fill <- self[['matrix']]$dims[1]
796       dims.fill <- (dims.fill + 1L):(dims.fill + num.cells)
797       # Matrix data
798       if (display.progress) {
799         cate("Adding data to /matrix")
800       }
801       matrix.data <- t(x = as.matrix(x = data.frame(matrix.data)))
802       self[['matrix']][dims.fill, ] <- matrix.data
803       # Layer data
804       if (display.progress) {
805         cate("Adding data to /layers")
806         pb <- new.pb()
807         counter <- 0
808       }
809       layers.names <- names(x = self[['layers']])
810       for (i in layers.names) {
811         self[['layers']][[i]][dims.fill, ] <- t(x = layers.data[[i]])
812         if (display.progress) {
813           counter <- counter + 1
814           setTxtProgressBar(pb = pb, value = counter / length(x = layers.names))
815         }
816       }
817       # Column attributes
818       if (display.progress) {
819         cate("Adding data to /col_attrs")
820         pb <- new.pb()
821         counter <- 0
822       }
823       attrs.names <- names(x = self[['col_attrs']])
824       for (i in attrs.names) {
825         self[['col_attrs']][[i]][dims.fill] <- attributes.data[[i]]
826         if (display.progress) {
827           counter <- counter + 1
828           setTxtProgressBar(pb = pb, value = counter / length(x = attrs.names))
829         }
830       }
831       # # Load layers and attributes
832       # private$load_layers()
833       # private$load_attributes(MARGIN = 1)
834       # private$load_attributes(MARGIN = 2)
835       # self$shape <- self[['matrix']]$dims
836     },
837     add.loom = function(
838       other,
839       other.key = NULL,
840       self.key = NULL,
841       ...
842     ) {
843       if (self$mode == 'r') {
844         stop(private$err_mode)
845       }
846       # Connect to the other loom file
847       if (inherits(x = other, what = 'loom')) {
848         ofile <- other
849       } else if (is.character(x = other)) {
850         ofile <- connect(filename = other)
851       } else {
852         stop("'other' must be either a loom object or a path to a loom file")
853       }
854       # If we have row keys to use
855       if (!is.null(x = other.key) && !is.null(x = self.key)) {
856         other.key <- basename(path = other.key)
857         self.key <- basename(path = self.key)
858         tryCatch(
859           expr = other.key <- other[['row_attrs']][[other.key]][],
860           error = function(e) {
861             if (is.character(x = other)) {
862               ofile$close
863             }
864             stop("Failed to find the gene names dataset in the other loom file")
865           }
866         )
867         tryCatch(
868           expr = self.key <- self[['row_attrs']][[self.key]][],
869           error = function(e) {
870             if (is.character(x = other)) {
871               ofile$close
872             }
873             stop("Failed to find the gene names dataset in this loom file")
874           }
875         )
876         # Match rows
877         rows.use <- match(x = other.key, table = self.key, nomatch = 0)
878         rows.use <- rows.use[rows.use > 0]
879       } else {
880         cate("Adding the loom file as-is, assuming in the correct order")
881         Sys.sleep(time = 3)
882         rows.use <- 1:other[['matrix']]$dims[2]
883       }
884       if (max(rows.use) > self[['matrix']]$dims[2]) {
885         stop("More genes in the other loom file than present in this one")
886       } else if (max(rows.use) < self[['matrix']]$dims[2]) {
887         cate("...")
888       }
889       # Clean up
890       if (is.character(x = other)) {
891         ofile$close_all()
892       }
893     }
894   ),
895   # Private fields and methods
896   # @field err_init An error message for if this object hasn't been created with loomR::create or loomR::connect
897   # @field err_mode An error message for if this object is in read-only mode
898   # @field it Iterator object for batch.scan and batch.next
899   # @field iter.dataset Dataset for iterating on
900   # @field iter.margin Margin to iterate over
901   # @field iter.index Index values for iteration
902   # @field skipped.validation Was validation skipped?
903   # \describe{
904   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
905   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
906   #   \item{\code{reset_batch()}}{Reset the batch iterator fields}
907   #   \item{\code{iter_range(index.use)}}{Get the range of indices for a batch iteration}
908   # }
909   private = list(
910     # Fields
911     err_init = "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",
912     err_mode = "Cannot modify a loom file in read-only mode",
913     it = NULL,
914     iter.chunksize = NULL,
915     iter.dataset = NULL,
916     iter.margin = NULL,
917     iter.index = NULL,
918     skipped.validation = FALSE,
919     # Methods
920     load_attributes = function(MARGIN) {
921       attribute <- switch(
922         EXPR = MARGIN,
923         '1' = 'row_attrs',
924         '2' = 'col_attrs',
925         stop('Invalid attribute dimension')
926       )
927       group <- self[[attribute]]
928       attributes <- unlist(x = lapply(
929         X = names(x = group),
930         FUN = function(x) {
931           d <- list(group[[x]])
932           names(x = d) <- x
933           return(d)
934         }
935       ))
936       if (MARGIN == 1) {
937         self$row.attrs <- attributes
938       } else if (MARGIN == 2) {
939         self$col.attrs <- attributes
940       }
941     },
942     load_layers = function() {
943       if (private$skipped.validation) {
944         if (getOption(x = 'verbose')) {
945           warning("Not loading layers field")
946         }
947       } else {
948         self$layers <- unlist(x = lapply(
949           X = names(x = self[['layers']]),
950           FUN = function(n) {
951             d <- list(self[[paste('layers', n, sep = '/')]])
952             names(x = d) <- n
953             return(d)
954           }
955         ))
956       }
957     },
958     reset_batch = function() {
959       private$it <- NULL
960       private$iter.chunksize <- NULL
961       private$iter.dataset <- NULL
962       private$iter.margin <- NULL
963       private$iter.index <- NULL
964     },
965     iter_range = function(index.use) {
966       if (is.null(private$iter.margin)) {
967         stop("Batch processing hasn't been set up")
968       }
969       if (is.null(x = index.use)) {
970         # If no index was provided, use entire range for this margin
971         index.use <- c(1, self$shape[private$iter.margin])
972       } else if (length(x = index.use) == 1) {
973         # If one index was provided, start at one and go to index
974         index.use <- c(1, index.use)
975       } else {
976         # Use index.use[1] and index.use[2]
977         index.use <- c(index.use[1], index.use[2])
978       }
979       # Ensure the indices provided fit within the range of the dataset
980       index.use[1] <- max(1, index.use[1])
981       index.use[2] <- min(index.use[2], self$shape[private$iter.margin])
982       # Ensure that index.use[1] is greater than index.use[2]
983       if (index.use[1] > index.use[2]) {
984         stop(paste0(
985           "Starting index (",
986           index.use[1],
987           ") must be lower than the ending index (",
988           index.use[2],
989           ")"
990         ))
991       }
992       return(index.use)
993     }
994   )
995 )
996
997 #' Create a loom object
998 #'
999 #' @param filename The name of the new loom file
1000 #' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
1001 #' @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}
1002 #' @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}
1003 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
1004 #' @param chunk.size How many rows of \code{data} should we stream to the loom file at any given time?
1005 #' @param overwrite Overwrite an already existing loom file?
1006 #'
1007 #' @return A connection to a loom file
1008 #'
1009 #' @importFrom utils packageVersion txtProgressBar setTxtProgressBar
1010 #'
1011 #' @seealso \code{\link{loom-class}}
1012 #'
1013 #' @export
1014 #'
1015 create <- function(
1016   filename,
1017   data,
1018   gene.attrs = NULL,
1019   cell.attrs = NULL,
1020   layers = NULL,
1021   chunk.dims = 'auto',
1022   chunk.size = 1000,
1023   overwrite = FALSE,
1024   display.progress = TRUE
1025 ) {
1026   mode <- ifelse(test = overwrite, yes = 'w', no = 'w-')
1027   if (file.exists(filename) && !overwrite) {
1028     stop(paste('File', filename, 'already exists!'))
1029   }
1030   if (!inherits(x = data, what = c('matrix', 'Matrix'))) {
1031     data <- as.matrix(x = data)
1032   }
1033   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
1034     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
1035   } else if (length(x = chunk.dims) == 1) {
1036     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
1037       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
1038     }
1039   } else {
1040     chunk.dims <- as.integer(x = chunk.dims)
1041   }
1042   new.loom <- loom$new(filename = filename, mode = mode)
1043   # Create the matrix
1044   # new.loom$create_dataset(
1045   #   name = 'matrix',
1046   #   robj = data,
1047   #   chunk_dims = chunk.dims
1048   # )
1049   dtype <- getDtype(x = data[1, 1])
1050   new.loom$create_dataset(
1051     name = 'matrix',
1052     dtype = dtype,
1053     dims = dim(x = data)
1054   )
1055   chunk.points <- chunkPoints(
1056     data.size = dim(x = data)[1],
1057     chunk.size = chunk.size
1058   )
1059   if (display.progress) {
1060     pb <- txtProgressBar(char = '=', style = 3)
1061   }
1062   for (col in 1:ncol(x = chunk.points)) {
1063     row.start <- chunk.points[1, col]
1064     row.end <- chunk.points[2, col]
1065     new.loom[['matrix']][row.start:row.end, ] <- as.matrix(x = data[row.start:row.end, ])
1066     if (display.progress) {
1067       setTxtProgressBar(pb = pb, value = col / ncol(x = chunk.points))
1068     }
1069   }
1070   new.loom$matrix <- new.loom[['matrix']]
1071   new.loom$shape <- rev(x = new.loom[['matrix']]$dims)
1072   # Groups
1073   new.loom$create_group(name = 'layers')
1074   new.loom$create_group(name = 'row_attrs')
1075   new.loom$create_group(name = 'col_attrs')
1076   # Check for the existance of gene or cell names
1077   if (!is.null(x = colnames(x = data))) {
1078     new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
1079   }
1080   if (!is.null(x = rownames(x = data))) {
1081     new.loom$add.col.attribute(attribute = list('cell_names' = rownames(x = data)))
1082   }
1083   # Store some constants as HDF5 attributes
1084   h5attr(x = new.loom, which = 'version') <- new.loom$version
1085   h5attr(x = new.loom, which = 'chunks') <- paste0(
1086     '(',
1087     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
1088     ')'
1089   )
1090   # Add layers
1091   if (!is.null(x = layers)) {
1092     new.loom$add.layer(layer = layers)
1093   }
1094   if (!is.null(x = gene.attrs)) {
1095     new.loom$add.row.attribute(attribute = gene.attrs)
1096   }
1097   if (!is.null(x = cell.attrs)) {
1098     new.loom$add.col.attribute(attribute = cell.attrs)
1099   }
1100   # Set last bit of information
1101   chunks <- new.loom[['matrix']]$chunk_dims
1102   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
1103   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
1104   chunks <- unlist(x = strsplit(x = chunks, split = ','))
1105   new.loom$chunksize <- as.integer(x = chunks)
1106   # Return the connection
1107   return(new.loom)
1108 }
1109
1110 #' Connect to a loom file
1111 #'
1112 #' @param filename The loom file to connect to
1113 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write.
1114 #' If \code{mode} is 'r+', loomR will automatically add missing required groups during validation
1115 #' @param skip.validate Skip the validation steps, use only for extremely large loom files
1116 #'
1117 #' @return A loom file connection
1118 #'
1119 #' @seealso \code{\link{loom-class}}
1120 #'
1121 #' @export
1122 #'
1123 connect <- function(filename, mode = "r", skip.validate = FALSE) {
1124   if (!(mode %in% c('r', 'r+'))) {
1125     stop("'mode' must be one of 'r' or 'r+'")
1126   }
1127   new.loom <- loom$new(filename = filename, mode = mode, skip.validate = skip.validate)
1128   return(new.loom)
1129 }
1130
1131 #' Subset a loom file
1132 #'
1133 #' @param x A loom object
1134 #' @param m Rows (genes) to subset, defaults to all rows
1135 #' @param n Columns (cells) to subset, defaults to all columns
1136 #' @param filename Filename for new loom object, defaults to ...
1137 #' @param overwrite Overwrite \code{filename} if already exists?
1138 #' @param display.progress Display progress as we're copying over data
1139 #' @param ... Ignored for now
1140 #'
1141 #' @return A loom object connected to \code{filename}
1142 #'
1143 #' @importFrom utils setTxtProgressBar
1144 #'
1145 #' @export subset.loom
1146 #' @method subset loom
1147 #'
1148 subset.loom <- function(
1149   x,
1150   m = NULL,
1151   n = NULL,
1152   filename = NULL,
1153   overwrite = FALSE,
1154   display.progress = TRUE,
1155   ...
1156 ) {
1157   # Set some defaults
1158   if (is.null(x = m)) {
1159     m <- 1:x$shape[1]
1160   }
1161   if (is.null(x = n)) {
1162     n <- 1:x$shape[2]
1163   }
1164   if (is.null(x = filename)) {
1165     filename <- paste(
1166       unlist(x = strsplit(x = x$filename, split = '.', fixed = TRUE)),
1167       collapse = '_subset.'
1168     )
1169   }
1170   if (length(x = m) == 1) {
1171     m <- 1:m
1172   }
1173   if (length(x = n) == 1) {
1174     n <- 1:n
1175   }
1176   # Ensure that m and n are within the bounds of the loom file
1177   if (max(m) > x$shape[1] || max(n) > x$shape[2]) {
1178     stop(paste(
1179       "'m' and 'n' must be less than",
1180       x$shape[1],
1181       "and",
1182       x$shape[2],
1183       "respectively"
1184     ))
1185   }
1186   extension <- rev(x = unlist(x = strsplit(x = filename, split = '.', fixed = TRUE)))[1]
1187   # Ensure that we call our new file a loom file
1188   if (extension != 'loom') {
1189     filename <- paste0(filename, '.loom')
1190   }
1191   if (display.progress) {
1192     catn("Writing new loom file to", filename)
1193   }
1194   # Make the loom file
1195   new.loom <- create(
1196     filename = filename,
1197     # data = t(x = x[['matrix']][n, m]),
1198     data = x[['matrix']][n, m],
1199     overwrite = overwrite
1200   )
1201   # Add row attributes
1202   row.attrs <- list.datasets(object = x, path = 'row_attrs', full.names = TRUE)
1203   if (length(x = row.attrs) > 0) {
1204     if (display.progress) {
1205       catn("\nAdding", length(x = row.attrs), "row attributes")
1206       pb <- new.pb()
1207       counter <- 0
1208     }
1209     for (row.attr in row.attrs) {
1210       row.list <- list(x[[row.attr]][m])
1211       names(x = row.list) <- basename(path = row.attr)
1212       new.loom$add.row.attribute(attribute = row.list)
1213       if (display.progress) {
1214         counter <- counter + 1
1215         setTxtProgressBar(pb = pb, value = counter / length(x = row.attrs))
1216       }
1217     }
1218   } else {
1219     warning("No row attributes found")
1220   }
1221   # Add col attributes
1222   col.attrs <- list.datasets(object = x, path = 'col_attrs', full.names = TRUE)
1223   if (length(x = col.attrs) > 0) {
1224     if (display.progress) {
1225       catn("\nAdding", length(x = col.attrs), "column attributes")
1226       pb <- new.pb()
1227       counter <- 0
1228     }
1229     for (col.attr in col.attrs) {
1230       col.list <- list(x[[col.attr]][n])
1231       names(x = col.list) <- basename(path = col.attr)
1232       new.loom$add.col.attribute(attribute = col.list)
1233       if (display.progress) {
1234         counter <- counter + 1
1235         setTxtProgressBar(pb = pb, value = counter / length(x = col.attrs))
1236       }
1237     }
1238   } else {
1239     warning("No column attributes found")
1240   }
1241   # Add layers
1242   layers <- list.datasets(object = x, path = 'layers', full.names = TRUE)
1243   if (length(x = layers) > 0) {
1244     if (display.progress) {
1245       catn("\nAdding", length(x = layers), "layers")
1246       pb <- new.pb()
1247       counter <- 0
1248     }
1249     for (layer in layers) {
1250       layer.list <- list(x[[layer]][n, m])
1251       names(x = layer.list) <- basename(path = layer)
1252       new.loom$add.layer(layers = layer.list)
1253       if (display.progress) {
1254         counter <- counter + 1
1255         setTxtProgressBar(pb = pb, value = counter / length(x = layers))
1256       }
1257     }
1258   } else {
1259     warning("No layers found")
1260   }
1261   new.loom$flush()
1262   return(new.loom)
1263 }
1264
1265 #' Combine loom files
1266 #'
1267 #' @param x A vector of loom files or filenames
1268 #' @param filename Name for resultant vector
1269 #' @param chunk.size How many rows form each input loom should we stream to the merged loom file at any given time?
1270 #' @param order.by Optional row attribute to order each input loom by, must be one dimensional
1271 #' @param overwrite Overwrite \code{filename} if already exists?
1272 #' @param display.progress Display progress as we're copying over data
1273 #'
1274 #' @return A loom object connected to \code{filename}
1275 #'
1276 #' @importFrom utils setTxtProgressBar
1277 #'
1278 #' @export Combine
1279 #'
1280 combine <- function(
1281   looms,
1282   filename,
1283   chunk.size = 1000,
1284   order.by = NULL,
1285   overwrite = FALSE,
1286   display.progress = TRUE,
1287   ...
1288 ) {
1289   # Basic checking of input arguments
1290   looms <- looms[vapply(
1291     X = looms,
1292     FUN = inherits,
1293     FUN.VALUE = logical(length = 1L),
1294     what = c('loom', 'character')
1295   )]
1296   if (length(x = looms) < 2) {
1297     stop("Need at least two loom objects or files to merge")
1298   }
1299   # Check the existance of loom files
1300   loom.names <- looms[is.character(x = looms)]
1301   if (length(x = loom.names) > 0) {
1302     if (!all(file.exists(loom.names))) {
1303       stop(paste0(
1304         "Cannot find the following loom files: '",
1305         paste(loom.names[!file.exists(loom.names)], collapse = "', '"),
1306         "'"
1307       ))
1308     }
1309   }
1310   # Set mode and provide more useful error
1311   mode <- ifelse(test = overwrite, yes = 'w', no = 'w-')
1312   if (file.exists(filename) && !overwrite) {
1313     stop(paste('File', filename, 'already exists!'))
1314   }
1315   # Check loom contents
1316   # Every loom must have same number of genes (rows, MARGIN = 2)
1317   # and same datasets in the groups
1318   row.attrs <- vector(mode = 'list', length = length(x = looms))
1319   row.types <- list()
1320   col.attrs <- vector(mode = 'list', length = length(x = looms))
1321   col.types <- list()
1322   layers <- vector(mode = 'list', length = length(x = looms))
1323   layers.types <- list()
1324   nrows <- vector(mode = 'list', length = length(x = looms))
1325   ncols <- vector(mode = 'list', length = length(x = looms))
1326   for (i in 1:length(x = looms)) {
1327     this <- if (is.character(x = looms[i])) {
1328       connect(filename = looms[i])
1329     } else {
1330       looms[i]
1331     }
1332     row.attrs[[i]] <- sort(x = list.datasets(
1333       object = this,
1334       path = 'row_attrs',
1335       full.names = TRUE
1336     ))
1337     for (attr in row.attrs[[i]]) {
1338       if (length(x = attr) > 0) {
1339         row.types[[attr]] <- c(
1340           row.types[[attr]],
1341           as.character(x = this[[attr]]$get_type()$get_class())
1342         )
1343       }
1344     }
1345     col.attrs[[i]] <- sort(x = list.datasets(
1346       object = this,
1347       path = 'col_attrs',
1348       full.names = TRUE
1349     ))
1350     for (attr in col.attrs[[i]]) {
1351       if (length(x = attr) > 0) {
1352         col.types[[attr]] <- c(
1353           col.types[[attr]],
1354           as.character(x = this[[attr]]$get_type()$get_class())
1355         )
1356       }
1357     }
1358     layers[[i]] <- sort(x = list.datasets(
1359       object = this,
1360       path = 'layers',
1361       full.names = TRUE
1362     ))
1363     for (lay in layers) {
1364       if (length(x = lay) > 0) {
1365         layers.types[[lay]] <- c(
1366           layers.types[[attr]],
1367           as.character(x = this[[lay]]$get_type()$get_class())
1368         )
1369       }
1370     }
1371     nrows[[i]] <- this[['matrix']]$dims[2]
1372     ncols[[i]] <- this[['matrix']]$dims[1]
1373     if (is.character(x = looms[i])) {
1374       this$close_all()
1375     }
1376   }
1377   row.attrs <- unique(x = row.attrs)
1378   col.attrs <- unique(x = col.attrs)
1379   layers <- unique(x = layers)
1380   nrows <- unique(x = nrows)
1381   ncells <- sum(unlist(x = ncols))
1382   if (length(x = row.attrs) != 1) {
1383     stop("Not all loom objects have the same row attributes")
1384   }
1385   if (length(x = col.attrs) != 1) {
1386     stop("Not all loom objects have the same column attributes")
1387   }
1388   if (length(x = layers) != 1) {
1389     stop("Not all loom objects have the same layers")
1390   }
1391   if (length(x = nrows) != 1) {
1392     stop("Not all loom objects have the number of rows (MARGIN = 2)")
1393   }
1394   # Check for the row attribute to order by
1395   if (!is.null(x = order.by)) {
1396     if (!grepl(pattern = order.by, x = row.attrs)) {
1397       stop(paste0("Cannot find '", order.by, "' in the row attributes for the loom files provided"))
1398     } else {
1399       temp <- if (is.character(x = looms[1])) {
1400         connect(filename = looms[1])
1401       } else {
1402         looms[1]
1403       }
1404       order.dat <- temp[['row_attrs']][[basename(path = order.by)]]
1405       if (length(x = order.dat$dims) != 1) {
1406         if (is.character(x = looms[1])) {
1407           temp$close_all()
1408         }
1409         stop("'order.by' must reference a one dimensional attribute")
1410       }
1411       order.use <- order.dat[]
1412       if (is.character(x = looms[1])) {
1413         temp$close_all()
1414       }
1415     }
1416   }
1417   # Check data types:
1418   row.types <- lapply(X = row.types, FUN = unique)
1419   row.types.counts <- vapply(X = row.types, FUN = length, FUN.VALUE = integer(length = 1L))
1420   col.types <- lapply(X = col.types, FUN = unique)
1421   col.types.counts <- vapply(X = col.types, FUN = length, FUN.VALUE = integer(length = 1L))
1422   layers.types <- lapply(X = layers.types, FUN = unique)
1423   layers.types.counts <- vapply(X = layers.types, FUN = length, FUN.VALUE = integer(length = 1L))
1424   if (any(row.types.counts > 1)) {
1425     stop(paste0(
1426       "The following row attributes have multiple types across the input loom files: '",
1427       paste(names(x = row.types.counts[row.types.counts > 1]), collapse = "', '"),
1428       "'; cannot combine"
1429     ))
1430   }
1431   if (any(col.types.counts > 1)) {
1432     stop(paste0(
1433       "The following column attributes have multiple types across the input loom files: '",
1434       paste(names(x = col.types.counts[col.types.counts > 1]), collapse = "', '"),
1435       "'; cannot combine"
1436     ))
1437   }
1438   if (any(layers.types.counts > 1)) {
1439     stop(paste0(
1440       "The following layers have multiple types across the input loom files: '",
1441       paste(names(x = layers.types.counts[layers.types.counts > 1]), collapse = "', '"),
1442       "'; cannot combine"
1443     ))
1444   }
1445   # # Create new HDF5 file
1446   # hfile <- h5file(filename = filename, mode = mode)
1447   # Start adding loom objects
1448   for (i in 1:length(x = looms)) {
1449     catn("Adding loom file", i ,"of", length(x = looms))
1450     chunk.points <- chunkPoints(data.size = ncols[[i]], chunk.size = chunk.size)
1451     next
1452   }
1453   # # Close and reopen as loom object for returning
1454   # hfile$close_all()
1455   # lfile <- connect(filename = filename, mode = 'r+')
1456   # return(lfile)
1457 }
1458
1459 # #need to comment
1460 # #need to add progress bar
1461 # #but otherwise, pretty cool
1462 # #for paul to try :
1463 # # f <- connect("~/Downloads/10X43_1.loom")
1464 # # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
1465 # # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
1466 # map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
1467 #   n_func = length(f_list)
1468 #   if (n_func == 1) {
1469 #     f_list <- list(f_list)
1470 #   }
1471 #   if (MARGIN == 1) {
1472 #     results <- list()
1473 #     for (j in 1:n_func) {
1474 #       results[[j]] <- numeric(0)
1475 #     }
1476 #     rows_per_chunk <- chunksize
1477 #     ix <- 1
1478 #     while (ix <= self@shape[1]) {
1479 #       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
1480 #       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
1481 #       for (j in 1:n_func) {
1482 #         new_results <- apply(chunk, 1, FUN = f_list[[j]])
1483 #         results[[j]] <- c(results[[j]], new_results)
1484 #       }
1485 #       ix <- ix + chunksize
1486 #     }
1487 #   }
1488 #   if (MARGIN == 2) {
1489 #     results <- list()
1490 #     for (j in 1:n_func) {
1491 #       results[[j]] <- numeric(0)
1492 #     }
1493 #     cols_per_chunk <- chunksize
1494 #     ix <- 1
1495 #     while (ix <= self@shape[2]) {
1496 #       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
1497 #       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
1498 #       for (j in 1:n_func) {
1499 #         new_results <- apply(chunk, 2, FUN = f_list[[j]])
1500 #         results[[j]] <- c(results[[j]], new_results)
1501 #       }
1502 #       ix <- ix + chunksize
1503 #     }
1504 #   }
1505 #   if (n_func == 1) {
1506 #     results <- results[[1]]
1507 #   }
1508 #   return(results)
1509 # }