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