set stringsAsFactors = FALSE in get.attribute.df
[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       # by default return all attributes
392       if (is.null(attribute.names)) {
393         attribute.names <- self[[attribute.layer]]$names
394       }
395       # check that attribute names are present
396       if (!all(attribute.names %in% self[[attribute.layer]]$names)) {
397         invalid.names <- attribute.names[which(!attribute.names %in% self[[attribute.layer]]$names)]
398         stop(paste0("Invalid attribute.names: ", paste0(invalid.names, collapse = ", ")))
399       }
400       attr.paths <- paste0(attribute.layer, "/", attribute.names)
401       # keep only the one-dimensional attributes
402       dims <- sapply(attr.paths, function(x) length(self[[x]]$dims))
403       data.lst <- lapply(attr.paths[dims == 1], function(x) data.frame(self[[x]][], stringsAsFactors = FALSE))
404       combined.df <- Reduce(cbind, data.lst)
405       colnames(combined.df) <- attribute.names[dims == 1]
406       if (attribute.layer == "col_attrs") {
407         rownames(combined.df) <- self[[paste0(attribute.layer, "/", col.names)]][]
408       } else {
409         rownames(combined.df) <- self[[paste0(attribute.layer, "/", row.names)]][]
410       }
411       # check if any row is all NAs
412       rows.to.remove <- unname(obj = which(x = apply(
413         X = combined.df,
414         MARGIN = 1,
415         FUN = function(x) {
416           return(all(is.na(x = x)))
417         }
418       )))
419       if (length(x = rows.to.remove) > 1) {
420         combined.df <- combined.df[-rows.to.remove, ]
421       }
422       return(combined.df)
423     },
424     # Chunking functions
425     batch.scan = function(
426       chunk.size = NULL,
427       MARGIN = 2,
428       index.use = NULL,
429       dataset.use = 'matrix',
430       force.reset = FALSE
431     ) {
432       if (is.null(x = private$it) || !grepl(pattern = dataset.use, x = private$iter.dataset) || force.reset) {
433         # Check the existence of the dataset
434         private$iter.dataset <- grep(
435           pattern = dataset.use,
436           x = list.datasets(object = self),
437           value = TRUE
438         )
439         if (length(x = private$iter.dataset) != 1) {
440           private$reset_batch()
441           stop(paste0("Cannot find dataset '", dataset.use, "' in the loom file"))
442         }
443         if (grepl(pattern = 'row_attrs', x = private$iter.dataset)) {
444           MARGIN <- 1
445         }
446         else if (grepl(pattern = 'col_attrs', x = private$iter.dataset)) {
447           MARGIN <- 2
448         }
449         # Check the margin
450         if (!(MARGIN %in% c(1, 2))) {
451           private$reset_batch()
452           stop("MARGIN must be 1 (genes) or 2 (cells)")
453         } else {
454           private$iter.margin <- MARGIN
455         }
456         if (is.null(x = chunk.size)) {
457           chunk.size <- rev(x = self$chunksize)[private$iter.margin]
458         }
459         private$iter.chunksize <- chunk.size
460         # Set the indices to use
461         index.use <- private$iter_range(index.use = index.use)
462         # Setup our iterator
463         private$it <- ihasNext(iterable = ichunk(
464           iterable = index.use[1]:index.use[2],
465           chunkSize = chunk.size
466         ))
467         private$iter.index <- index.use
468       }
469       # Return the times we iterate, this isn't important, we only need the length of this vector
470       return(1:ceiling((private$iter.index[2] - private$iter.index[1]) / private$iter.chunksize))
471     },
472     batch.next = function(return.data = TRUE) {
473       # Ensure that we have a proper iterator
474       if (!'hasNext.ihasNext' %in% suppressWarnings(expr = methods(class = class(x = private$it)))) {
475         stop("Please setup the iterator with self$batch.scan")
476       }
477       # Do the iterating
478       if (hasNext(obj = private$it)) {
479         # Get the indices for this chunk
480         chunk.indices <- unlist(x = nextElem(obj = private$it))
481         if (return.data) {
482           # If we're working with a matrix dataset, ensure chunking on proper axis
483           if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
484             to.return <- switch(
485               EXPR = private$iter.margin,
486               '1' = self[[private$iter.dataset]][, chunk.indices], # Genes
487               '2' = self[[private$iter.dataset]][chunk.indices, ] # Cells
488             )
489           } else {
490             # Otherwise, iterating over an attribute (1 dimensional)
491             to.return <- self[[private$iter.dataset]][chunk.indices]
492           }
493         } else {
494           to.return <- chunk.indices
495         }
496         # Determine if we reset the iterator
497         if (!hasNext(obj = private$it)) {
498           private$reset_batch()
499         }
500         return(to.return)
501       } else {
502         # Just in case
503         private$reset_batch()
504         return(NULL)
505       }
506     },
507     apply = function(
508       name,
509       FUN,
510       MARGIN = 2,
511       index.use = NULL,
512       chunk.size = NULL,
513       dataset.use = 'matrix',
514       overwrite = FALSE,
515       display.progress = TRUE,
516       ...
517     ) {
518       if (self$mode == 'r') {
519         stop(private$err_mode)
520       }
521       if (!inherits(x = FUN, what = 'function')) {
522         stop("FUN must be a function")
523       }
524       # Check that we're storing our results properly
525       results.basename <- basename(path = name) # Basename of our results
526       results.dirname <- gsub(pattern = '/', replacement = '', x = dirname(path = name)) # Groupname of our results
527       dirnames <- c('row_attrs', 'col_attrs', 'layers') # Allowed group names
528       if (name %in% list.datasets(object = self)) {
529         if (overwrite) {
530           self$link_delete(name = name)
531         } else {
532           stop(paste("A dataset with the name", name, "already exists!"))
533         }
534       }
535       # Checks datset, index, and MARGIN
536       # Sets full dataset path in private$iter.dataset
537       # Sets proper MARGIN in private$iter.margin
538       batch <- self$batch.scan(
539         chunk.size = chunk.size,
540         MARGIN = MARGIN,
541         dataset.use = dataset.use,
542         force.reset = TRUE
543       )
544       MARGIN <- private$iter.margin
545       dataset.use <- private$iter.dataset
546       # Ensure that our group name is allowed
547       name.check <- which(x = dirnames == results.dirname)
548       if (!any(name.check)) {
549         private$reset_batch()
550         stop(paste(
551           "The dataset must go into one of:",
552           paste(dirnames, collapse = ', ')
553         ))
554       }
555       # Check that our group matches our iteration
556       # ie. To store in col_attrs, MARGIN must be 1
557       if (name.check %in% c(1, 2) && name.check != private$iter.margin) {
558         private$reset_batch()
559         stop(paste(
560           "Iteration must be over",
561           c('genes', 'cells')[name.check],
562           paste0("(MARGIN = ", name.check, ")"),
563           "to store results in",
564           paste0("'", dirnames[name.check], "'")
565         ))
566       }
567       # Check how we store our results
568       dataset.matrix <- ('matrix' %in% private$iter.dataset || grepl(pattern = 'layers', x = private$iter.dataset))
569       results.matrix <- name.check == 3
570       # Ensure index.use is integers within the bounds of [1, self$shape[MARGIN]]
571       if (!is.null(x = index.use)) {
572         # Filter index.use to values between 1 and self$shape[MARGIN]
573         index.use <- as.integer(x = index.use)
574         index.use[index.use >= 1 & index.use <= self$shape[MARGIN]]
575         index.use <- as.vector(x = na.omit(object = index.use))
576         # If we still have values, figure out NAs, otherwise set index.use to NULL
577         if (length(x = index.use) > 0) {
578
579         } else {
580           warning("No values passed to 'index.use' fall within the data, using all values")
581           index.use <- NULL
582         }
583       }
584       # Trial to get class of new dataset
585       # Do a trial run to figure out the class of dataset
586       na.use <- NA
587       if (display.progress) {
588         catn("Running trial to determine class of dataset")
589       }
590       trial.use <- if (is.null(x = index.use)) {
591         sample(x = 1:self$shape[MARGIN], size = 3, replace = FALSE)
592       } else {
593         sample(x = index.use, size = 3, replace = FALSE)
594       }
595       trial.use <- sort(x = trial.use)
596       trial <- if (grepl(pattern = 'layers', x = dataset.use) || dataset.use == 'matrix') {
597         switch(
598           EXPR = MARGIN,
599           '1' = self[[dataset.use]][, trial.use],
600           '2' = self[[dataset.use]][trial.use, ]
601         )
602       } else {
603         self[[dataset.use]][trial.use]
604       }
605       trial <- FUN(trial, ...)
606       if (is.list(x = trial)) {
607         trial <- unlist(x = trial)
608       }
609       trial <- as.vector(x = trial)
610       class(x = na.use) <- class(x = trial)
611       # Get a connection to the group we're iterating over
612       group <- self[[results.dirname]]
613       # Make results dataset
614       dtype.use <- getDtype(x = na.use)
615       dims.use <- switch(
616         EXPR = results.dirname,
617         'layers' = self[['matrix']]$dims,
618         'row_attrs' = self[['matrix']]$dims[2],
619         'col_attrs' = self[['matrix']]$dims[1]
620       )
621       group$create_dataset(
622         name = results.basename,
623         dtype = dtype.use,
624         dims = dims.use
625       )
626       # Start the iteration
627       if (display.progress) {
628         catn("Writing results to", name)
629         pb <- txtProgressBar(char = '=', style = 3)
630       }
631       # first <- TRUE
632       for (i in 1:length(x = batch)) {
633         # Get the indices we're iterating over
634         these.indices <- self$batch.next(return.data = FALSE)
635         if (is.null(x = index.use)) {
636           chunk.indices <- these.indices
637         } else {
638           chunk.indices <- index.use[index.use %in% these.indices]
639           chunk.na <- these.indices[!(these.indices %in% chunk.indices)]
640         }
641         # Get the data and apply FUN
642         chunk.data <- if (dataset.matrix) {
643           switch(
644             EXPR = MARGIN,
645             '1' = self[[dataset.use]][, chunk.indices], # Chunk genes
646             '2' = self[[dataset.use]][chunk.indices, ] # Chunk cells
647           )
648         } else {
649           self[[private$iter.datset]][chunk.indices]
650         }
651         chunk.data <- FUN(chunk.data, ...)
652         if (results.matrix) {
653           # If we're writign to a matrix
654           # Figure out which way we're writing the data
655           switch(
656             EXPR = MARGIN,
657             '1' = group[[results.basename]][, chunk.indices] <- chunk.data,
658             '2' = group[[results.basename]][chunk.indices, ] <- chunk.data
659           )
660           if (!is.null(x = index.use)) {
661             switch(
662               EXPR = MARGIN,
663               '1' = group[[results.basename]][, chunk.na] <- na.use,
664               '2' = group[[results.basename]][chunk.na, ] <- na.use
665             )
666           }
667         } else {
668           # Just write to the vector
669           group[[results.basename]][chunk.indices] <- chunk.data
670           if (!is.null(x = index.use)) {
671             group[[results.basename]][chunk.na] <- na.use
672           }
673         }
674         if (display.progress) {
675           setTxtProgressBar(pb = pb, value = i / length(x = batch))
676         }
677       }
678       if (display.progress) {
679         close(con = pb)
680       }
681       # Clean up and allow chaining
682       private$reset_batch()
683       # Load layers and attributes
684       private$load_layers()
685       private$load_attributes(MARGIN = 1) # Genes (row_attrs)
686       private$load_attributes(MARGIN = 2) # Cells (col_attrs)
687       invisible(x = self)
688     },
689     map = function(
690       FUN,
691       MARGIN = 2,
692       chunk.size = NULL,
693       index.use = NULL,
694       dataset.use = 'matrix',
695       display.progress = TRUE,
696       expected = NULL,
697       ...
698     ) {
699       if (!inherits(x = FUN, what = 'function')) {
700         stop("FUN must be a function")
701       }
702       # Checks datset, index, and MARGIN
703       # Sets full dataset path in private$iter.dataset
704       # Sets proper MARGIN in private$iter.margin
705       batch <- self$batch.scan(
706         chunk.size = chunk.size,
707         MARGIN = MARGIN,
708         index.use = index.use,
709         dataset.use = dataset.use,
710         force.reset = TRUE
711       )
712       MARGIN <- private$iter.margin
713       dataset.use <- private$iter.dataset
714       # Check how we store our results
715       # And what the shape of our dataset is
716       results.matrix <- TRUE
717       dataset.matrix <- TRUE
718       if (grepl(pattern = 'col_attrs', x = dataset.use) || grepl(pattern = 'row_attrs', x = dataset.use)) {
719         results.matrix <- FALSE
720         dataset.matrix <- FALSE
721       }
722       if (!is.null(x = expected)) {
723         results.matrix <- switch(
724           EXPR = expected,
725           'vector' = FALSE,
726           'matrix' = TRUE,
727           stop("'expected' must be one of 'matrix', 'vector', or NULL")
728         )
729       }
730       # Create our results holding object
731       results <- vector(mode = "list", length = length(x = batch))
732       if (display.progress) {
733         pb <- txtProgressBar(char = '=', style = 3)
734       }
735       for (i in 1:length(x = batch)) {
736         chunk.indices <- self$batch.next(return.data = FALSE)
737         chunk.data <- if (dataset.matrix) {
738           switch(
739             EXPR = MARGIN,
740             '1' = self[[dataset.use]][, chunk.indices], # Chunk genes
741             '2' = self[[dataset.use]][chunk.indices, ] # Chunk cells
742           )
743         } else {
744           self[[dataset.use]][chunk.indices]
745         }
746         results[[i]] <- FUN(chunk.data, ...)
747         if (display.progress) {
748           setTxtProgressBar(pb = pb, value = i / length(x = batch))
749         }
750       }
751       # Bring result list into vector or matrix format
752       if (class(results[[1]]) == 'numeric') {
753         results <- unlist(x = results, use.names = FALSE)
754       } else {
755         if (MARGIN == 1) {
756           results <- Reduce(f = cbind, x = results)
757         } else if (MARGIN == 2) {
758           results <- Reduce(f = rbind, x = results)
759         }
760       }
761       if (display.progress) {
762         close(con = pb)
763       }
764       private$reset_batch()
765       return(results)
766     },
767     # Functions that modify `/matrix'
768     add.cells = function(
769       matrix.data,
770       attributes.data = NULL,
771       layers.data = NULL,
772       display.progress = TRUE
773     ) {
774       if (self$mode == 'r') {
775         stop(private$err_mode)
776       }
777       # Check inputs
778       n <- self[['matrix']]$dims[2]
779       if (display.progress) {
780         cate("Checking inputs...")
781       }
782       matrix.data <- check.matrix_data(x = matrix.data, n = n)
783       layers.data <- check.layers(
784         x = layers.data,
785         n = n,
786         layers.names = names(x = self[['layers']])
787       )
788       attributes.data <- check.col_attrs(
789         x = attributes.data,
790         attrs.names = names(x = self[['col_attrs']])
791       )
792       # Get the number of cells we're adding
793       num.cells <- c(
794         nCells.matrix_data(x = matrix.data),
795         nCells.layers(x = layers.data),
796         nCells.col_attrs(x = attributes.data)
797       )
798       num.cells <- max(num.cells)
799       # Flesh out the input data to add
800       if (display.progress) {
801         cate(paste("Adding", num.cells, "to this loom file"))
802       }
803       matrix.data <- addCells.matrix_data(x = matrix.data, n = n, m2 = num.cells)
804       layers.data <- addCells.layers(x = layers.data, n = n, m2 = num.cells)
805       attributes.data <- addCells.col_attrs(x = attributes.data, m2 = num.cells)
806       # Add the input to the loom file
807       dims.fill <- self[['matrix']]$dims[1]
808       dims.fill <- (dims.fill + 1L):(dims.fill + num.cells)
809       # Matrix data
810       if (display.progress) {
811         cate("Adding data to /matrix")
812       }
813       matrix.data <- t(x = as.matrix(x = data.frame(matrix.data)))
814       self[['matrix']][dims.fill, ] <- matrix.data
815       # Layer data
816       if (display.progress) {
817         cate("Adding data to /layers")
818         pb <- new.pb()
819         counter <- 0
820       }
821       layers.names <- names(x = self[['layers']])
822       for (i in layers.names) {
823         self[['layers']][[i]][dims.fill, ] <- t(x = layers.data[[i]])
824         if (display.progress) {
825           counter <- counter + 1
826           setTxtProgressBar(pb = pb, value = counter / length(x = layers.names))
827         }
828       }
829       # Column attributes
830       if (display.progress) {
831         cate("Adding data to /col_attrs")
832         pb <- new.pb()
833         counter <- 0
834       }
835       attrs.names <- names(x = self[['col_attrs']])
836       for (i in attrs.names) {
837         self[['col_attrs']][[i]][dims.fill] <- attributes.data[[i]]
838         if (display.progress) {
839           counter <- counter + 1
840           setTxtProgressBar(pb = pb, value = counter / length(x = attrs.names))
841         }
842       }
843       # # Load layers and attributes
844       # private$load_layers()
845       # private$load_attributes(MARGIN = 1)
846       # private$load_attributes(MARGIN = 2)
847       # self$shape <- self[['matrix']]$dims
848     },
849     add.loom = function(
850       other,
851       other.key = NULL,
852       self.key = NULL,
853       ...
854     ) {
855       if (self$mode == 'r') {
856         stop(private$err_mode)
857       }
858       # Connect to the other loom file
859       if (inherits(x = other, what = 'loom')) {
860         ofile <- other
861       } else if (is.character(x = other)) {
862         ofile <- connect(filename = other)
863       } else {
864         stop("'other' must be either a loom object or a path to a loom file")
865       }
866       # If we have row keys to use
867       if (!is.null(x = other.key) && !is.null(x = self.key)) {
868         other.key <- basename(path = other.key)
869         self.key <- basename(path = self.key)
870         tryCatch(
871           expr = other.key <- other[['row_attrs']][[other.key]][],
872           error = function(e) {
873             if (is.character(x = other)) {
874               ofile$close
875             }
876             stop("Failed to find the gene names dataset in the other loom file")
877           }
878         )
879         tryCatch(
880           expr = self.key <- self[['row_attrs']][[self.key]][],
881           error = function(e) {
882             if (is.character(x = other)) {
883               ofile$close
884             }
885             stop("Failed to find the gene names dataset in this loom file")
886           }
887         )
888         # Match rows
889         rows.use <- match(x = other.key, table = self.key, nomatch = 0)
890         rows.use <- rows.use[rows.use > 0]
891       } else {
892         cate("Adding the loom file as-is, assuming in the correct order")
893         Sys.sleep(time = 3)
894         rows.use <- 1:other[['matrix']]$dims[2]
895       }
896       if (max(rows.use) > self[['matrix']]$dims[2]) {
897         stop("More genes in the other loom file than present in this one")
898       } else if (max(rows.use) < self[['matrix']]$dims[2]) {
899         cate("...")
900       }
901       # Clean up
902       if (is.character(x = other)) {
903         ofile$close_all()
904       }
905     }
906   ),
907   # Private fields and methods
908   # @field err_init An error message for if this object hasn't been created with loomR::create or loomR::connect
909   # @field err_mode An error message for if this object is in read-only mode
910   # @field it Iterator object for batch.scan and batch.next
911   # @field iter.dataset Dataset for iterating on
912   # @field iter.margin Margin to iterate over
913   # @field iter.index Index values for iteration
914   # @field skipped.validation Was validation skipped?
915   # \describe{
916   #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
917   #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
918   #   \item{\code{reset_batch()}}{Reset the batch iterator fields}
919   #   \item{\code{iter_range(index.use)}}{Get the range of indices for a batch iteration}
920   # }
921   private = list(
922     # Fields
923     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",
924     err_mode = "Cannot modify a loom file in read-only mode",
925     it = NULL,
926     iter.chunksize = NULL,
927     iter.dataset = NULL,
928     iter.margin = NULL,
929     iter.index = NULL,
930     skipped.validation = FALSE,
931     # Methods
932     load_attributes = function(MARGIN) {
933       attribute <- switch(
934         EXPR = MARGIN,
935         '1' = 'row_attrs',
936         '2' = 'col_attrs',
937         stop('Invalid attribute dimension')
938       )
939       group <- self[[attribute]]
940       if (length(x = names(x = group)) > 0) {
941         attributes <- unlist(x = lapply(
942           X = names(x = group),
943           FUN = function(x) {
944             d <- list(group[[x]])
945             names(x = d) <- x
946             return(d)
947           }
948         ))
949       } else {
950         attributes <- NULL
951       }
952       if (MARGIN == 1) {
953         self$row.attrs <- attributes
954       } else if (MARGIN == 2) {
955         self$col.attrs <- attributes
956       }
957     },
958     load_layers = function() {
959       if (private$skipped.validation) {
960         if (getOption(x = 'verbose')) {
961           warning("Not loading layers field")
962         }
963       } else if (length(x = names(x = self[['layers']])) > 0) {
964         self$layers <- unlist(x = lapply(
965           X = names(x = self[['layers']]),
966           FUN = function(n) {
967             d <- list(self[[paste('layers', n, sep = '/')]])
968             names(x = d) <- n
969             return(d)
970           }
971         ))
972       }
973     },
974     reset_batch = function() {
975       private$it <- NULL
976       private$iter.chunksize <- NULL
977       private$iter.dataset <- NULL
978       private$iter.margin <- NULL
979       private$iter.index <- NULL
980     },
981     iter_range = function(index.use) {
982       if (is.null(private$iter.margin)) {
983         stop("Batch processing hasn't been set up")
984       }
985       if (is.null(x = index.use)) {
986         # If no index was provided, use entire range for this margin
987         index.use <- c(1, self$shape[private$iter.margin])
988       } else if (length(x = index.use) == 1) {
989         # If one index was provided, start at one and go to index
990         index.use <- c(1, index.use)
991       } else {
992         # Use index.use[1] and index.use[2]
993         index.use <- c(index.use[1], index.use[2])
994       }
995       # Ensure the indices provided fit within the range of the dataset
996       index.use[1] <- max(1, index.use[1])
997       index.use[2] <- min(index.use[2], self$shape[private$iter.margin])
998       # Ensure that index.use[1] is greater than index.use[2]
999       if (index.use[1] > index.use[2]) {
1000         stop(paste0(
1001           "Starting index (",
1002           index.use[1],
1003           ") must be lower than the ending index (",
1004           index.use[2],
1005           ")"
1006         ))
1007       }
1008       return(index.use)
1009     }
1010   )
1011 )
1012
1013 #' Create a loom object
1014 #'
1015 #' @param filename The name of the new loom file
1016 #' @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}
1017 #' @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}
1018 #' @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}
1019 #' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
1020 #' @param do.transpose Transpose the input? Should be \code{TRUE} if \code{data} has genes as rows and cells as columns
1021 #' @param chunk.size How many rows of \code{data} should we stream to the loom file at any given time?
1022 #' @param overwrite Overwrite an already existing loom file?
1023 #'
1024 #' @return A connection to a loom file
1025 #'
1026 #' @importFrom utils packageVersion txtProgressBar setTxtProgressBar
1027 #'
1028 #' @seealso \code{\link{loom-class}}
1029 #'
1030 #' @export
1031 #'
1032 create <- function(
1033   filename,
1034   data,
1035   gene.attrs = NULL,
1036   cell.attrs = NULL,
1037   layers = NULL,
1038   chunk.dims = 'auto',
1039   chunk.size = 1000,
1040   do.transpose = TRUE,
1041   overwrite = FALSE,
1042   display.progress = TRUE
1043 ) {
1044   mode <- ifelse(test = overwrite, yes = 'w', no = 'w-')
1045   if (file.exists(filename) && !overwrite) {
1046     stop(paste('File', filename, 'already exists!'))
1047   }
1048   if (!inherits(x = data, what = c('matrix', 'Matrix'))) {
1049     data <- as.matrix(x = data)
1050   }
1051   if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
1052     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
1053   } else if (length(x = chunk.dims) == 1) {
1054     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
1055       chunk.dims <- rep.int(x = as.integer(x = chunk.dims), times = 2)
1056     }
1057   } else {
1058     chunk.dims <- as.integer(x = chunk.dims)
1059   }
1060   new.loom <- loom$new(filename = filename, mode = mode)
1061   dtype <- getDtype(x = data[1, 1])
1062   matrix.shape <- dim(x = data)
1063   if (do.transpose) {
1064     cate("Transposing input data: loom file will show input columns (cells) as rows and input rows (genes) as columns")
1065     cate("This is to maintain compatibility with other loom tools")
1066     matrix.shape <- rev(x = matrix.shape)
1067   } else {
1068     cate("Not tranposing data: loom file will show data exactly like input")
1069     cate("Please note, other loom tools will show this flipped")
1070   }
1071   new.loom$create_dataset(
1072     name = 'matrix',
1073     dtype = dtype,
1074     dims = matrix.shape
1075   )
1076   chunk.points <- chunkPoints(
1077     data.size = matrix.shape[1],
1078     chunk.size = chunk.size
1079   )
1080   if (display.progress) {
1081     pb <- txtProgressBar(char = '=', style = 3)
1082   }
1083   for (col in 1:ncol(x = chunk.points)) {
1084     row.start <- chunk.points[1, col]
1085     row.end <- chunk.points[2, col]
1086     data.add <- if (do.transpose) {
1087       t(x = as.matrix(x = data[, row.start:row.end]))
1088     } else {
1089       as.matrix(x = data[row.start:row.end, ])
1090     }
1091     # new.loom[['matrix']][row.start:row.end, ] <- as.matrix(x = data[row.start:row.end, ])
1092     new.loom[['matrix']][row.start:row.end, ] <- data.add
1093     if (display.progress) {
1094       setTxtProgressBar(pb = pb, value = col / ncol(x = chunk.points))
1095     }
1096   }
1097   new.loom$matrix <- new.loom[['matrix']]
1098   new.loom$shape <- rev(x = new.loom[['matrix']]$dims)
1099   # Groups
1100   groups <- c('layers', 'row_attrs', 'col_attrs', 'row_graphs', 'col_graphs')
1101   for (group in groups) {
1102     new.loom$create_group(name = group)
1103   }
1104   # Check for the existance of gene or cell names
1105   if (!is.null(x = colnames(x = data))) {
1106     if (do.transpose) {
1107       new.loom$add.col.attribute(attribute = list('cell_names' = colnames(x = data)))
1108     } else {
1109       new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
1110     }
1111   }
1112   if (!is.null(x = rownames(x = data))) {
1113     if (do.transpose) {
1114       new.loom$add.row.attribute(attribute = list('gene_names' = rownames(x = data)))
1115     } else {
1116       new.loom$add.col.attribute(attribute = list('cell_names' = rownames(x = data)))
1117     }
1118   }
1119   # Store some constants as HDF5 attributes
1120   h5attr(x = new.loom, which = 'version') <- new.loom$version
1121   h5attr(x = new.loom, which = 'chunks') <- paste0(
1122     '(',
1123     paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
1124     ')'
1125   )
1126   # Add layers
1127   if (!is.null(x = layers)) {
1128     new.loom$add.layer(layer = layers)
1129   }
1130   if (!is.null(x = gene.attrs)) {
1131     new.loom$add.row.attribute(attribute = gene.attrs)
1132   }
1133   if (!is.null(x = cell.attrs)) {
1134     new.loom$add.col.attribute(attribute = cell.attrs)
1135   }
1136   # Set last bit of information
1137   chunks <- new.loom[['matrix']]$chunk_dims
1138   chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
1139   chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
1140   chunks <- unlist(x = strsplit(x = chunks, split = ','))
1141   new.loom$chunksize <- as.integer(x = chunks)
1142   # Return the connection
1143   return(new.loom)
1144 }
1145
1146 #' Connect to a loom file
1147 #'
1148 #' @param filename The loom file to connect to
1149 #' @param mode How do we connect to it? Pass 'r' for read-only or 'r+' for read/write.
1150 #' If \code{mode} is 'r+', loomR will automatically add missing required groups during validation
1151 #' @param skip.validate Skip the validation steps, use only for extremely large loom files
1152 #'
1153 #' @return A loom file connection
1154 #'
1155 #' @seealso \code{\link{loom-class}}
1156 #'
1157 #' @export
1158 #'
1159 connect <- function(filename, mode = "r", skip.validate = FALSE) {
1160   if (!(mode %in% c('r', 'r+'))) {
1161     stop("'mode' must be one of 'r' or 'r+'")
1162   }
1163   new.loom <- loom$new(filename = filename, mode = mode, skip.validate = skip.validate)
1164   return(new.loom)
1165 }
1166
1167 #' Subset a loom file
1168 #'
1169 #' @param x A loom object
1170 #' @param m Rows (cells) to subset, defaults to all rows
1171 #' @param n Columns (genes) to subset, defaults to all columns
1172 #' @param filename Filename for new loom object, defaults to ...
1173 #' @param chunk.m Chunk size for m
1174 #' @param chunk.n Chunk size for n
1175 #' @param overwrite Overwrite \code{filename} if already exists?
1176 #' @param display.progress Display progress bars?
1177 #' @param ... Ignored for now
1178 #'
1179 #' @return A loom object connected to \code{filename}
1180 #'
1181 #' @aliases subset
1182 #'
1183 #' @importFrom utils setTxtProgressBar
1184 #'
1185 #' @export subset.loom
1186 #' @method subset loom
1187 #'
1188 subset.loom <- function(
1189   x,
1190   m = NULL,
1191   n = NULL,
1192   filename = NULL,
1193   chunk.m = 1000,
1194   chunk.n = 1000,
1195   overwrite = FALSE,
1196   display.progress = TRUE,
1197   ...
1198 ) {
1199   m.all <- is.null(x = m)
1200   n.all <- is.null(x = n)
1201   if (m.all && n.all) {
1202     stop("Subset is set to all data, not subsetting")
1203   }
1204   # Set some defaults
1205   if (m.all) { # Pull all cells
1206     m <- 1:x[['matrix']]$dims[1]
1207   }
1208   if (n.all) { # Pull all genes
1209     n <- 1:x[['matrix']]$dims[2]
1210   }
1211   if (is.null(x = filename)) { # Set default filename
1212     filename <- paste(
1213       unlist(x = strsplit(x = x$filename, split = '.', fixed = TRUE)),
1214       collapse = '_subset.'
1215     )
1216   }
1217   # Ensure that m and n are within the bounds of the loom file
1218   if (max(m) > x[['matrix']]$dims[1] || max(n) > x[['matrix']]$dims[2]) {
1219     stop(paste(
1220       "'m' and 'n' must be less than",
1221       x[['matrix']]$dims[1],
1222       "and",
1223       x[['matrix']]$dims[2],
1224       "respectively"
1225     ))
1226   }
1227   extension <- rev(x = unlist(x = strsplit(x = filename, split = '.', fixed = TRUE)))[1]
1228   # Ensure that we call our new file a loom file
1229   if (extension != 'loom') {
1230     filename <- paste0(filename, '.loom')
1231   }
1232   # Make the loom file
1233   mode <- ifelse(test = overwrite, yes = 'w', no = 'w-')
1234   if (file.exists(filename) && !overwrite) {
1235     stop(paste('File', filename, 'already exists!'))
1236   } else if (display.progress) {
1237     catn("Writing new loom file to", filename)
1238   }
1239   new.loom <- loom$new(filename = filename, mode = mode)
1240   # Add /matrix
1241   matrix.dims <- c(length(x = m), length(x = n))
1242   new.loom$create_dataset(
1243     name = 'matrix',
1244     dtype = getDtype2(x = class(x = x[['matrix']]$get_type())[1]),
1245     dims = matrix.dims
1246   )
1247   new.loom$shape <- rev(x = matrix.dims)
1248   batch <- new.loom$batch.scan(chunk.size = chunk.m)
1249   if (display.progress) {
1250     catn("Adding data for /matrix")
1251     pb <- new.pb()
1252   }
1253   for (i in 1:length(x = batch)) {
1254     these.indices <- new.loom$batch.next(return.data = FALSE)
1255     chunk.indices <- m[these.indices]
1256     new.loom[['matrix']][these.indices, ] <- x[['matrix']][chunk.indices, n]
1257     if (display.progress) {
1258       setTxtProgressBar(pb = pb, value = i / length(x = batch))
1259     }
1260   }
1261   if (display.progress) {
1262     close(con = pb)
1263   }
1264   # Add groups
1265   for (group in c('row_attrs', 'row_graphs', 'col_attrs', 'col_graphs', 'layers')) {
1266     new.loom$create_group(name = group)
1267   }
1268   # Add row attributes
1269   row.attrs <- list.datasets(object = x, path = 'row_attrs', full.names = TRUE)
1270   if (length(x = row.attrs) > 0) {
1271     if (display.progress) {
1272       catn("Adding", length(x = row.attrs), "row attributes")
1273       pb <- new.pb()
1274       counter <- 0
1275     }
1276     for (row.attr in row.attrs) {
1277       if (length(x = x[[row.attr]]$dims) == 2) {
1278         new.loom[[row.attr]] <- x[[row.attr]][, n]
1279       } else {
1280         new.loom[[row.attr]] <- x[[row.attr]][n]
1281       }
1282       if (display.progress) {
1283         counter <- counter + 1
1284         setTxtProgressBar(pb = pb, value = counter / length(x = row.attrs))
1285       }
1286     }
1287     if (display.progress) {
1288       close(con = pb)
1289     }
1290   } else {
1291     cate("No row attributes found")
1292   }
1293   # Add col attributes
1294   col.attrs <- list.datasets(object = x, path = 'col_attrs', full.names = TRUE)
1295   if (length(x = col.attrs) > 0) {
1296     if (display.progress) {
1297       catn("Adding", length(x = col.attrs), "column attributes")
1298       pb <- new.pb()
1299       counter <- 0
1300     }
1301     for (col.attr in col.attrs) {
1302       if (length(x = x[[col.attr]]$dims) == 2) {
1303         new.loom[[col.attr]] <- x[[col.attr]][, m]
1304       } else {
1305         new.loom[[col.attr]] <- x[[col.attr]][m]
1306       }
1307       if (display.progress) {
1308         counter <- counter + 1
1309         setTxtProgressBar(pb = pb, value = counter / length(x = col.attrs))
1310       }
1311     }
1312     if (display.progress) {
1313       close(con = pb)
1314     }
1315   } else {
1316     cate("No column attributes found")
1317   }
1318   # Add layers
1319   layers <- list.datasets(object = x, path = 'layers', full.names = TRUE)
1320   if (length(x = layers) > 0) {
1321     if (display.progress) {
1322       catn("Adding", length(x = layers), "layers")
1323       pb <- new.pb()
1324     }
1325     # Initialize datasets
1326     for (layer in layers) {
1327       new.loom[['layers']]$create_dataset(
1328         name = basename(path = layer),
1329         dtype = getDtype2(x = class(x = x[[layer]]$get_type())[1]),
1330         dims = matrix.dims
1331       )
1332     }
1333     # Start adding in batches
1334     batch <- new.loom$batch.scan(chunk.size = chunk.m)
1335     for (i in 1:length(x = batch)) {
1336       these.indices <- new.loom$batch.next(return.data = FALSE)
1337       chunk.indices <- m[these.indices]
1338       # Add the data for each loom for this batch
1339       for (layer in layers) {
1340         new.loom[[layer]][these.indices, ] <- x[[layer]][chunk.indices, n]
1341       }
1342       if (display.progress) {
1343         setTxtProgressBar(pb = pb, value = i / length(x = batch))
1344       }
1345     }
1346     if (display.progress) {
1347       close(con = pb)
1348     }
1349   } else {
1350     cate("No layers found")
1351   }
1352   new.loom$flush()
1353   new.loom$load.fields()
1354   return(new.loom)
1355 }
1356
1357 #' Combine loom files
1358 #'
1359 #' @param looms A vector of loom files or filenames
1360 #' @param filename Name for resultant vector
1361 #' @param chunk.size How many rows form each input loom should we stream to the merged loom file at any given time?
1362 #' @param order.by Optional row attribute to order each input loom by, must be one dimensional
1363 #' @param overwrite Overwrite \code{filename} if already exists?
1364 #' @param display.progress Display progress as we're copying over data
1365 #'
1366 #' @return A loom object connected to \code{filename}
1367 #'
1368 #' @importFrom utils setTxtProgressBar
1369 #'
1370 #' @seealso \code{\link{loom-class}}
1371 #'
1372 #' @export
1373 #'
1374 combine <- function(
1375   looms,
1376   filename,
1377   chunk.size = 1000,
1378   order.by = NULL,
1379   overwrite = FALSE,
1380   display.progress = TRUE,
1381   ...
1382 ) {
1383   # Basic checking of input arguments
1384   looms <- looms[vapply(
1385     X = looms,
1386     FUN = inherits,
1387     FUN.VALUE = logical(length = 1L),
1388     what = c('loom', 'character')
1389   )]
1390   if (length(x = looms) < 2) {
1391     stop("Need at least two loom objects or files to merge")
1392   }
1393   # Check the existance of loom files
1394   loom.names <- looms[is.character(x = looms)]
1395   if (length(x = loom.names) > 0) {
1396     if (!all(file.exists(loom.names))) {
1397       stop(paste0(
1398         "Cannot find the following loom files: '",
1399         paste(loom.names[!file.exists(loom.names)], collapse = "', '"),
1400         "'"
1401       ))
1402     }
1403   }
1404   # Set mode and provide more useful error
1405   mode <- ifelse(test = overwrite, yes = 'w', no = 'w-')
1406   if (file.exists(filename) && !overwrite) {
1407     stop(paste('File', filename, 'already exists!'))
1408   }
1409   # Check loom contents
1410   # Every loom must have same number of genes (rows, MARGIN = 2)
1411   # and same datasets in the groups
1412   row.attrs <- vector(mode = 'list', length = length(x = looms))
1413   row.types <- list()
1414   col.attrs <- vector(mode = 'list', length = length(x = looms))
1415   col.dims <- list()
1416   col.ndims <- list()
1417   col.types <- list()
1418   layers <- vector(mode = 'list', length = length(x = looms))
1419   layers.types <- list()
1420   nrows <- vector(mode = 'list', length = length(x = looms))
1421   ncols <- vector(mode = 'list', length = length(x = looms))
1422   matrix.type <- list()
1423   if (display.progress) {
1424     catn("Validating", length(x = looms), "input loom files")
1425     pb <- new.pb()
1426   }
1427   for (i in 1:length(x = looms)) {
1428     this <- if (is.character(x = looms[i])) {
1429       connect(filename = looms[i])
1430     } else {
1431       looms[i]
1432     }
1433     row.attrs[[i]] <- sort(x = list.datasets(
1434       object = this,
1435       path = 'row_attrs',
1436       full.names = TRUE
1437     ))
1438     for (attr in row.attrs[[i]]) {
1439       if (length(x = attr) > 0) {
1440         row.types[[attr]] <- c(
1441           row.types[[attr]],
1442           class(x = this[[attr]]$get_type())[1]
1443         )
1444       }
1445     }
1446     col.attrs[[i]] <- sort(x = list.datasets(
1447       object = this,
1448       path = 'col_attrs',
1449       full.names = TRUE
1450     ))
1451     for (attr in col.attrs[[i]]) {
1452       if (length(x = attr) > 0) {
1453         col.types[[attr]] <- c(
1454           col.types[[attr]],
1455           class(x = this[[attr]]$get_type())[1]
1456         )
1457         col.dims[[attr]] <- sum(col.dims[[attr]], this[[attr]]$dims[1])
1458         col.ndims[[attr]] <- c(
1459           col.ndims[[attr]],
1460           length(x = this[[attr]]$dims)
1461         )
1462       }
1463     }
1464     layers[[i]] <- sort(x = list.datasets(
1465       object = this,
1466       path = 'layers',
1467       full.names = TRUE
1468     ))
1469     for (lay in layers) {
1470       if (length(x = lay) > 0) {
1471         layers.types[[lay]] <- c(
1472           layers.types[[lay]],
1473           class(x = this[[lay]]$get_type())[1]
1474         )
1475       }
1476     }
1477     nrows[[i]] <- this[['matrix']]$dims[2]
1478     ncols[[i]] <- this[['matrix']]$dims[1]
1479     matrix.type[[i]] <- class(x = this[['matrix']]$get_type())[1]
1480     if (is.character(x = looms[i])) {
1481       this$close_all()
1482     }
1483     if (display.progress) {
1484       setTxtProgressBar(pb = pb, value = i / length(x = looms))
1485     }
1486   }
1487   row.attrs <- unique(x = row.attrs)
1488   col.attrs <- unique(x = col.attrs)
1489   layers <- unique(x = layers)
1490   nrows <- unlist(x = unique(x = nrows))
1491   ncols <- unlist(x = ncols)
1492   ncells <- sum(ncols)
1493   if (length(x = row.attrs) != 1) {
1494     stop("Not all loom objects have the same row attributes")
1495   }
1496   if (length(x = col.attrs) != 1) {
1497     stop("Not all loom objects have the same column attributes")
1498   }
1499   if (length(x = layers) != 1) {
1500     stop("Not all loom objects have the same layers")
1501   }
1502   if (length(x = nrows) != 1) {
1503     stop("Not all loom objects have the number of rows (MARGIN = 2)")
1504   }
1505   # Check for the row attribute to order by
1506   if (!is.null(x = order.by)) { # We have something to order by
1507     if (!grepl(pattern = order.by, x = row.attrs)) { # Check to ensure that the order.by is in each of our row attributes
1508       stop(paste0("Cannot find '", order.by, "' in the row attributes for the loom files provided"))
1509     } else {
1510       # If it is, get the values to order by
1511       order.by <- basename(path = order.by) # Basename of the order.by attribute name
1512       temp <- if (is.character(x = looms[1])) { # Connect to the loom first loom file
1513         connect(filename = looms[1])
1514       } else {
1515         looms[1]
1516       }
1517       order.dat <- temp[['row_attrs']][[order.by]] # Ensure that our order.by attribute is one dimmensional
1518       if (length(x = order.dat$dims) != 1) {
1519         if (is.character(x = looms[1])) {
1520           temp$close_all()
1521         }
1522         stop("'order.by' must reference a one dimensional attribute")
1523       }
1524       # If so, pull the data
1525       order.use <- order.dat[]
1526       # If the first loom was initially closed, close it again
1527       if (is.character(x = looms[1])) {
1528         temp$close_all()
1529       }
1530     }
1531   }
1532   # Check data types:
1533   matrix.type <- unlist(x = unique(x = matrix.type))
1534   row.types <- lapply(X = row.types, FUN = unique)
1535   row.types.counts <- vapply(X = row.types, FUN = length, FUN.VALUE = integer(length = 1L))
1536   col.types <- lapply(X = col.types, FUN = unique)
1537   col.types.counts <- vapply(X = col.types, FUN = length, FUN.VALUE = integer(length = 1L))
1538   layers.types <- lapply(X = layers.types, FUN = unique)
1539   layers.types.counts <- vapply(X = layers.types, FUN = length, FUN.VALUE = integer(length = 1L))
1540   if (any(row.types.counts > 1)) {
1541     stop(paste0(
1542       "The following row attributes have multiple types across the input loom files: '",
1543       paste(names(x = row.types.counts[row.types.counts > 1]), collapse = "', '"),
1544       "'; cannot combine"
1545     ))
1546   }
1547   if (any(col.types.counts > 1)) {
1548     stop(paste0(
1549       "The following column attributes have multiple types across the input loom files: '",
1550       paste(names(x = col.types.counts[col.types.counts > 1]), collapse = "', '"),
1551       "'; cannot combine"
1552     ))
1553   }
1554   if (any(layers.types.counts > 1)) {
1555     stop(paste0(
1556       "The following layers have multiple types across the input loom files: '",
1557       paste(names(x = layers.types.counts[layers.types.counts > 1]), collapse = "', '"),
1558       "'; cannot combine"
1559     ))
1560   }
1561   if (length(x = matrix.type) != 1) {
1562     stop("Cannot combine multiple datatypes for '/matrix'")
1563   }
1564   # Check dimmensions for col_attrs
1565   col.ndims <- lapply(X = col.ndims, FUN = unique)
1566   col.ndims.counts <- vapply(X = col.ndims, FUN = length, FUN.VALUE = integer(length = 1L))
1567   if (any(col.ndims.counts > 1)) {
1568     stop(paste0(
1569       "The following column attributes have varying dimmensions across the input loom files: '",
1570       paste(names(x = col.ndims.counts[col.ndims.counts > 1]), collapse = "', '"),
1571       "'; cannot combine"
1572     ))
1573   }
1574   # Create the new HDF5 file and the required groups
1575   new.loom <- loom$new(filename = filename, mode = mode)
1576   groups <- c('layers', "row_attrs", "col_attrs", "row_graphs", "col_graphs")
1577   for (group in groups) {
1578     new.loom$create_group(name = group)
1579   }
1580   # Initialize the '/matrix' dataset as well as datasets in 'col_attrs' and 'layers'
1581   new.loom$create_dataset(
1582     name = 'matrix', # Name is '/matrix'
1583     dtype = getDtype2(x = matrix.type), # Use the single type that we got from above
1584     dims = c(ncells, nrows) # Use the number of cells from the sum of ncols above, nrows should be the same for everyone
1585   )
1586   for (attr in col.attrs) {
1587     if (length(x = attr) > 0) {
1588       dims.use <- switch(
1589         EXPR = col.ndims[[attr]],
1590         '1' = ncells,
1591         '2' = c(col.dims[[attr]], ncells),
1592         stop("loomR supports only one- and two-dimmensional attribute datasets")
1593       )
1594       new.loom$create_dataset(
1595         name = attr,
1596         dtype = getDtype2(x = col.types[[attr]]),
1597         dims = dims.use
1598       )
1599     }
1600   }
1601   for (lay in layers) {
1602     if (length(x = lay) > 1) {
1603       new.loom$create_dataset(
1604         name = lay,
1605         dtype = getDtype2(x = layers.types[[lay]]),
1606         dims = c(ncells, nrows)
1607       )
1608     }
1609   }
1610   # Start adding loom objects
1611   matrix.previous <- 0
1612   col.previous <- vector(mode = 'integer', length = length(x = col.attrs))
1613   names(x = col.previous) <- unlist(x = col.attrs)
1614   for (i in 1:length(x = looms)) {
1615     if (display.progress) {
1616       catn("\nAdding loom file", i ,"of", length(x = looms))
1617     }
1618     # Open the loom file
1619     this <- if (is.character(x = looms[i])) {
1620       connect(filename = looms[i])
1621     } else {
1622       looms[i]
1623     }
1624     # Get the chunk points
1625     chunk.points <- chunkPoints(data.size = ncols[[i]], chunk.size = chunk.size)
1626     # print(chunk.points)
1627     # Pull ordering information
1628     order.genes <- if (is.null(x = order.by)) {
1629       1:nrows # If order.by wasn't provided, just use 1:number of genes
1630     } else {
1631       # Otherwise, match the data order from our order.use
1632       order(match(x = this[[row.attrs]][[order.by]][], table = order.use))
1633     }
1634     # Add data to /matrix, /col_attrs, and /layers
1635     if (display.progress) {
1636       catn("Adding data to /matrix and /layers")
1637       pb <- new.pb()
1638     }
1639     for (col in 1:ncol(x = chunk.points)) {
1640       cells.use <- chunk.points[1, col]:chunk.points[2, col]
1641       # Add /matrix for this chunk
1642       matrix.add <- this[['matrix']][cells.use, ]
1643       new.loom[['matrix']][cells.use + matrix.previous, ] <- matrix.add[, order.genes]
1644       # Add layers for this chunk
1645       for (lay in list.datasets(object = this[['layers']])) {
1646         lay.add <- this[['layers']][[lay]][cells.use, ]
1647         new.loom[['layers']][[lay]][cells.use + matrix.previous, ] <- lay.add[, order.genes]
1648       }
1649       if (display.progress) {
1650         setTxtProgressBar(pb = pb, value = col / ncol(x = chunk.points))
1651       }
1652       gc()
1653     }
1654     matrix.previous <- matrix.previous + max(chunk.points)
1655     # Add col_attrs for this chunk
1656     if (display.progress) {
1657       catn("\nAdding data to /col_attrs")
1658       pb <- new.pb()
1659     }
1660     for (attr in list.datasets(object = this[['col_attrs']], full.names = TRUE)) {
1661       start <- col.previous[attr]
1662       end <- start + this[[attr]]$dims[length(x = this[[attr]]$dims)]
1663       if (col.ndims[[attr]] == 1) {
1664         new.loom[[attr]][(start + 1):end] <- this[[attr]][]
1665       } else {
1666         for (col in 1:ncol(x = chunk.points)) {
1667           col.use <- chunk.points[1, col]:chunk.points[2, col]
1668           new.loom[[attr]][col.use, (start + 1):end] <- this[[attr]][col.use, ]
1669         }
1670       }
1671       col.previous[attr] <- end
1672       if (display.progress) {
1673         setTxtProgressBar(
1674           pb = pb,
1675           value = grep(
1676             pattern = attr,
1677             x = list.datasets(object = this[['col_attrs']], full.names = TRUE)
1678           ) / length(x = list.datasets(object = this[['col_attrs']], full.names = TRUE))
1679         )
1680       }
1681     }
1682     # Copy row attributes from the first loom object into the merged one
1683     if (i == 1) {
1684       if (display.progress) {
1685         catn("\nCopying data for /row_attrs")
1686         pb <- new.pb()
1687       }
1688       for (attr in list.datasets(object = this, path = 'row_attrs')) {
1689         new.loom[['row_attrs']]$obj_copy_from(
1690           src_loc = this[['row_attrs']],
1691           src_name = attr,
1692           dst_name = attr
1693         )
1694         if (display.progress) {
1695           setTxtProgressBar(
1696             pb = pb,
1697             value = grep(
1698               pattern = attr,
1699               x = list.datasets(object = this[['row_attrs']])
1700             ) / length(x = list.datasets(object = this[['row_attrs']]))
1701           )
1702         }
1703       }
1704     }
1705     new.loom$flush()
1706     # Close current loom file, if not open previously
1707     if (is.character(x = looms[i])) {
1708       this$close_all()
1709     }
1710   }
1711   return(new.loom)
1712 }
1713
1714 # #need to comment
1715 # #need to add progress bar
1716 # #but otherwise, pretty cool
1717 # #for paul to try :
1718 # # f <- connect("~/Downloads/10X43_1.loom")
1719 # # mean_var = map(f,f_list = c(mean,var),chunksize = 5000)
1720 # # nGene <- map(f, f_list = function(x) length(which(x>0)), MARGIN = 2)
1721 # map <- function(self, f_list = list(mean, var), MARGIN=1, chunksize=1000, selection) {
1722 #   n_func = length(f_list)
1723 #   if (n_func == 1) {
1724 #     f_list <- list(f_list)
1725 #   }
1726 #   if (MARGIN == 1) {
1727 #     results <- list()
1728 #     for (j in 1:n_func) {
1729 #       results[[j]] <- numeric(0)
1730 #     }
1731 #     rows_per_chunk <- chunksize
1732 #     ix <- 1
1733 #     while (ix <= self@shape[1]) {
1734 #       rows_per_chunk <- min(rows_per_chunk, self@shape[1] - ix + 1)
1735 #       chunk <- self["matrix"][ix:(ix + rows_per_chunk - 1), ]
1736 #       for (j in 1:n_func) {
1737 #         new_results <- apply(chunk, 1, FUN = f_list[[j]])
1738 #         results[[j]] <- c(results[[j]], new_results)
1739 #       }
1740 #       ix <- ix + chunksize
1741 #     }
1742 #   }
1743 #   if (MARGIN == 2) {
1744 #     results <- list()
1745 #     for (j in 1:n_func) {
1746 #       results[[j]] <- numeric(0)
1747 #     }
1748 #     cols_per_chunk <- chunksize
1749 #     ix <- 1
1750 #     while (ix <= self@shape[2]) {
1751 #       cols_per_chunk <- min(cols_per_chunk, self@shape[2] - ix + 1)
1752 #       chunk <- self["matrix"][, ix:(ix + cols_per_chunk - 1)]
1753 #       for (j in 1:n_func) {
1754 #         new_results <- apply(chunk, 2, FUN = f_list[[j]])
1755 #         results[[j]] <- c(results[[j]], new_results)
1756 #       }
1757 #       ix <- ix + chunksize
1758 #     }
1759 #   }
1760 #   if (n_func == 1) {
1761 #     results <- results[[1]]
1762 #   }
1763 #   return(results)
1764 # }