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