Fix issues with create, further progress on batching
[loomr.git] / R / loom.R
index f1f316138303e36f70890a8f00daf9bc0b44cd31..758a19f114aa8a1526f32d331f08797ea68a5dbe 100644 (file)
--- a/R/loom.R
+++ b/R/loom.R
@@ -2,24 +2,63 @@
 #' @importFrom R6 R6Class
 NULL
 
-#' A class for loom
+#' A class for loom files
 #'
 #' @docType class
 #' @name loom-class
 #' @rdname loom-class
-#' @return Object of class \code{\link{loom}}
+#' @return Object of \code{\link{R6::R6Class}} to generate \code{loom} objects
+#' @format An \code{\link{R6::R6Class}} object
 #' @seealso \code{\link{hdf5r::H5File}}
 #'
+#' @field version Version of loomR object was created under
+#' @field shape Shape of \code{/matrix} in columns (cells) by rows (genes)
+#' @field chunksize Chunks set for this dataset in columns (cells) by rows (genes)
+#' @field matrix The main data matrix, stored as columns (cells) by rows (genes)
+#' @field layers Additional data matricies, the same shape as \code{/matrix}
+#' @field col.attrs Extra information about cells
+#' @field row.attrs Extra information about genes
+#'
+#' @section Methods:
+#' \describe{
+#'   \item{\code{add.layer(layer)}}{Add a data layer to this loom file, must be in column (cells) by row (genes) orientation}
+#'   \item{\code{add.attribute(attribute, MARGIN)}}{
+#'     Add extra information to this loom file where
+#'     \code{attribute} is a named list where each element is a vector that is as long as one dimension of \code{/matrix} and
+#'     \code{MARGIN} is either 1 for cells or 2 for genes
+#'   }
+#'   \item{\code{add.row.attribute(attribute)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 2)}}
+#'   \item{\code{add.col.attribute(attribute)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 1)}}
+#'   \item{\code{add.meta.data(meta.data)}}{A wrapper for \code{add.attribute(attribute, MARGIN = 1)}}
+#'   \item{\code{batch.scan(chunk.size, MARGIN, index.use, dataset.use)}, \code{batch.next}}{
+#'     Scan a dataset in the loom file from \code{index.use[1]} to \code{index.use[2]}, iterating by \code{chunk.size}.
+#'     \code{dataset.use} can be the name, not \code{group/name}, unless the name is present in multiple groups.
+#'     If the dataset is in col_attrs, pass \code{MARGIN = 1}; if in row_attrs, pass \code{MARGIN = 2}.
+#'     Otherwise, pass \code{MARGIN = 1} to iterate on cells or \code{MARGIN = 2} to iterate on genes.
+#'     \code{chunk.size} defaults to \code{self$chunksize}, \code{MARGIN} defaults to 1,
+#'     \code{index.use} defaults to \code{1:self$shape[MARGIN]}, \code{dataset.use} defaults to 'matrix'
+#'   }
+#' }
+#'
+#' @importFrom iterators nextElem
 #' @importFrom utils packageVersion
+#' @importFrom itertools ihasNext ichunk hasNext
 #'
 #' @export
 #'
 loom <- R6Class(
+  # The loom class
+  # Based on the H5File class from hdf5r
+  # Not clonable (no loom$clone method), doesn't make sense since all data is on disk, not in memory
+  # Yes to portability, other packages can subclass the loom class
+  # Class is locked, other fields and methods cannot be added
   classname = 'loom',
   inherit = hdf5r::H5File,
   cloneable = FALSE,
   portable = TRUE,
   lock_class = TRUE,
+  # Public fields and methods
+  # See above for documentation
   public = list(
     # Fields
     version = NULL,
@@ -30,21 +69,22 @@ loom <- R6Class(
     col.attrs = NULL,
     row.attrs = NULL,
     # Methods
-    initialize = function(filename = NULL, mode = c('a', 'r', 'r+', 'w', 'w+'), ...) {
+    initialize = function(filename = NULL, mode = c('a', 'r', 'r+', 'w', 'w-'), ...) {
       # If the file exists, run validation steps
       do.validate <- file.exists(filename) && !(mode %in% c('w', 'w+'))
       super$initialize(filename = filename, mode = mode, ...)
       if (do.validate) {
         # Run the validation steps
         validateLoom(object = self)
-        # Store the shape of /matrix
+        # Store /matrix and the shape of /matrix
+        self$matrix <- self[['matrix']]
         self$shape <- self[['matrix']]$dims
         # Store the chunk size
         chunks <- h5attr(x = self, which = 'chunks')
         chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
         chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
         chunks <- unlist(x = strsplit(x = chunks, split = ','))
-        self$chunks <- as.integer(x = chunks)
+        self$chunksize <- as.integer(x = chunks)
         # Store version information
         self$version <- as.character(x = tryCatch(
           # Try getting a version
@@ -65,35 +105,225 @@ loom <- R6Class(
         # Load layers
         private$load_layers()
         # Load attributes
-        private$load_attirubtes(MARGIN = 1)
-        private$load_attributes(MARGIN = 2)
+        private$load_attributes(MARGIN = 1) # Cells (col_attrs)
+        private$load_attributes(MARGIN = 2) # Genes (row_attrs)
       } else {
         # Assume new HDF5 file
-        self$version <- packageVersion(pkg = 'loomR')
+        self$version <- as.character(x = packageVersion(pkg = 'loomR'))
       }
     },
-    add.layer = function(layer) {
-      invisible(x = NULL)
+    add.layer = function(layers) {
+      # Value checking
+      if (!is.list(x = layers) || is.null(x = names(x = layers))) {
+        stop("'layers' must be a named list")
+      }
+      if (is.null(x = self$shape)) {
+        stop(private$err_msg)
+      }
+      # Add layers
+      for (i in 1:length(x = layers)) {
+        if (!is.matrix(x = layers[[i]])) {
+          layers[[i]] <- as.matrix(x = layers[[i]])
+        }
+        do.transpose <- FALSE
+        if (any(dim(x = layers[[i]]) != self$shape)) {
+          if (all(rev(x = dim(x = layers[[i]])) == self$shape)) {
+            do.transpose <- TRUE
+          } else {
+            stop(paste(
+              "All layers must have",
+              self$shape[1],
+              "rows for cells and",
+              self$shape[2],
+              "columns for genes"
+            ))
+          }
+        }
+        if (do.transpose) {
+          self[['layers', names(x = layers)[i]]] <- t(x = layers[[i]])
+        } else {
+          self[['layers', names(x = layers)[i]]] <- layers[[i]]
+        }
+      }
+      self$flush()
+      private$load_layers()
+      invisible(x = self)
     },
-    add.attribute = function(attribute, MARGIN = 1) {
-      invisible(x = NULL)
+    add.attribute = function(attribute, MARGIN) {
+      # Value checking
+      if (!is.list(x = attribute) || is.null(x = names(x = attribute))) {
+        stop("'attribute' must be a named list")
+      }
+      if (length(x = attribute) > 1) {
+        for (i in attribute) {
+          if (!is.vector(x = attribute)) {
+            stop("All attributes must be one-dimensional vectors")
+          }
+        }
+      }
+      if (length(x = which(x = names(x = attribute) != '')) != length(x = attribute)) {
+        stop("Not all attributes had names provided")
+      }
+      if (!MARGIN %in% c(1, 2)) {
+        stop("'MARGIN' must be 1 or 2")
+      }
+      # Add the attributes as datasets for our MARGIN's group
+      if (is.null(x = self$shape)) {
+        stop(private$err_msg)
+      }
+      grp.name <- c('col_attrs', 'row_attrs')[MARGIN]
+      grp <- self[[grp.name]]
+      for (i in 1:length(x = attribute)) {
+        if (length(attribute[[i]]) != self$shape[MARGIN])
+          stop(paste(
+            "All",
+            switch(EXPR = MARGIN, '1' = 'cell', '2' = 'gene'),
+            "attributes must be of length",
+            self$shape[MARGIN]
+          ))
+        grp[[names(x = attribute)[i]]] <- attribute[[i]]
+      }
+      self$flush()
+      gc(verbose = FALSE)
+      # Load the attributes for this margin
+      private$load_attributes(MARGIN = MARGIN)
+      invisible(x = self)
     },
+    # Add attributes for genes
     add.row.attribute = function(attribute) {
-      invisible(x = NULL)
+      self$add.attribute(attribute = attribute, MARGIN = 2)
+      invisible(x = self)
     },
+    # Add attributes for cells
     add.col.attribute = function(attribute) {
-      invisible(x = NULL)
+      self$add.attribute(attribute = attribute, MARGIN = 1)
+      invisible(x = self)
+    },
+    # Add metadata, follows cells
+    add.meta.data = function(meta.data) {
+      self$add.col.attribute(attribute = meta.data)
+      invisible(x = self)
+    },
+    # Batch scan
+    batch.scan = function(
+      chunk.size = NULL,
+      MARGIN = 1,
+      index.use = NULL,
+      dataset.use = 'matrix'
+    ) {
+      if (is.null(x = private$it) || !grepl(pattern = dataset.use, x = private$iter.dataset)) {
+        # Check the existence of the dataset
+        private$iter.dataset <- grep(
+          pattern = dataset.use,
+          x = list.datasets(object = self),
+          value = TRUE
+        )
+        if (length(x = private$iter.dataset) != 1) {
+          stop(paste0("Cannot find dataset '", dataset.use, "' in the loom file"))
+        }
+        # Check the margin
+        if (!(MARGIN %in% c(1, 2))) {
+          stop("MARGIN must be 1 (cells) or 2 (genes)")
+        } else {
+          private$iter.margin <- MARGIN
+        }
+        if (is.null(x = chunk.size)) {
+          chunk.size <- self$chunksize[MARGIN]
+        }
+        # Set the indices to use
+        if (is.null(x = index.use)) {
+          index.use <- c(1, self$shape[MARGIN])
+        } else if (length(x = index.use) == 1) {
+          index.use <- index.use <- 1:index.use
+        } else {
+          index.use <- c(index.use[1], index.use[2])
+        }
+        index.use[1] <- max(1, index.use[1])
+        index.use[2] <- min(index.use[2], self$shape[MARGIN])
+        if (index.use[1] > index.use[2]) {
+          stop(paste0(
+            "Starting index (",
+            index.use[1],
+            ") must be lower than the ending index (",
+            index.use[2],
+            ")"
+          ))
+        }
+        # Setup our iterator
+        private$it <- ihasNext(iterable = ichunk(
+          iterable = index.use[1]:index.use[2],
+          chunkSize = chunk.size
+        ))
+        private$iter.index <- c(index.use[1], ceiling(x = index.use[2] / chunk.size))
+      }
+      return(private$iter.index[1]:private$iter.index[2])
+      # # Do the iterating
+      # if (hasNext(obj = private$it)) {
+      #   chunk.indices <- unlist(x = nextElem(obj = private$it))
+      #   if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
+      #     return(switch(
+      #       EXPR = private$iter.margin,
+      #       '1' = self[[private$iter.dataset]][chunk.indices, ],
+      #       '2' = self[[private$iter.dataset]][, chunk.indices]
+      #     ))
+      #   } else {
+      #     return(self[[private$iter.dataset]][chunk.indices])
+      #   }
+      # } else {
+      #   private$it <- NULL
+      #   return(NULL)
+      # }
+    },
+    batch.next = function() {
+      if (!'hasNext.ihasNext' %in% suppressWarnings(expr = methods(class = class(x = private$it)))) {
+        stop("Please setup the iterator with self$batch.scan")
+      }
+      # Do the iterating
+      if (hasNext(obj = private$it)) {
+        chunk.indices <- unlist(x = nextElem(obj = private$it))
+        if (private$iter.dataset == 'matrix' || grepl(pattern = 'layers', x = private$iter.dataset)) {
+          to.return <- switch(
+            EXPR = private$iter.margin,
+            '1' = self[[private$iter.dataset]][chunk.indices, ],
+            '2' = self[[private$iter.dataset]][, chunk.indices]
+          )
+        } else {
+          to.return <- self[[private$iter.dataset]][chunk.indices]
+        }
+        if (!hasNext(obj = private$it)) {
+          private$reset_batch()
+        }
+        return(to.return)
+      } else {
+        private$reset_batch()
+        return(NULL)
+      }
     }
   ),
+  # Private fields and methods
+  # @field err_msg A simple error message if this object hasn't been created with loomR::create or loomR::connect
+  # @field it
+  # @field iter.dataset
+  # @field iter.margin
+  # @field iter.index
+  # \describe{
+  #   \item{\code{load_attributes(MARGIN)}}{Load attributes of a given MARGIN into \code{self$col.attrs} or \code{self$row.attrs}}
+  #   \item{\code{load_layers()}}{Load layers into \code{self$layers}}
+  #   \item{\code{reset_batch()}}{Reset the batch iterator fields}
+  # }
   private = list(
-    add_attribute = function(attribute, MARGIN) {
-      invisible(x = NULL)
-    },
+    # Fields
+    err_msg = "This loom object has not been created with either loomR::create or loomR::connect, please use these function to create or connect to a loom file",
+    it = NULL,
+    iter.dataset = NULL,
+    iter.margin = NULL,
+    iter.index = NULL,
+    # Methods
     load_attributes = function(MARGIN) {
       attribute <- switch(
         EXPR = MARGIN,
-        '1' = 'row_attrs',
-        '2' = 'col_attrs',
+        '1' = 'col_attrs',
+        '2' = 'row_attrs',
         stop('Invalid attribute dimension')
       )
       group <- self[[attribute]]
@@ -105,25 +335,38 @@ loom <- R6Class(
           return(d)
         }
       ))
-      switch(
-        EXPR = MARGIN,
-        '1' = self$row.attrs <- attributes,
-        '2' = self$col.attrs <- attributes
-      )
+      if (MARGIN == 1) {
+        self$col.attrs <- attributes
+      } else if (MARGIN == 2) {
+        self$row.attrs <- attributes
+      }
     },
     load_layers = function() {
-      invisible(x = NULL)
+      self$layers <- unlist(x = lapply(
+        X = names(x = self[['layers']]),
+        FUN = function(n) {
+          d <- c(self[['layers', n]])
+          names(x = d) <- n
+          return(d)
+        }
+      ))
+    },
+    reset_batch = function() {
+      private$it <- NULL
+      private$iter.dataset <- NULL
+      private$iter.margin <- NULL
+      private$iter.index <- NULL
     }
   )
 )
 
 #' Create a loom object
 #'
-#' @param filename ...
-#' @param data ...
-#' @param row.attrs ...
-#' @param col.attrs ...
-#' @param chunk.dims ...
+#' @param filename The name of the new loom file
+#' @param data The data for \code{/matrix}, should be cells as rows and genes as columns
+#' @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}
+#' @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}
+#' @param chunk.dims A one- or two-length integer vector of chunksizes for \code{/matrix}, defaults to 'auto' to automatically determine chunksize
 #'
 #' @return A connection to a loom file
 #'
@@ -131,11 +374,13 @@ loom <- R6Class(
 #'
 #' @seealso \code{\link{loom-class}}
 #'
+#' @export
+#'
 create <- function(
   filename,
   data,
-  row.attrs = NULL,
-  col.attrs = NULL,
+  gene.attrs = NULL,
+  cell.attrs = NULL,
   layers = NULL,
   chunk.dims = 'auto'
 ) {
@@ -145,7 +390,7 @@ create <- function(
   if (!is.matrix(x = data)) {
     data <- as.matrix(x = data)
   }
-  if (length(x = chunk.dims) > 2 || length(x = chunk.dims < 1)) {
+  if (length(x = chunk.dims) > 2 || length(x = chunk.dims) < 1) {
     stop("'chunk.dims' must be a one- or two-length integer vector or 'auto'")
   } else if (length(x = chunk.dims == 1)) {
     if (!grepl(pattern = '^auto$', x = chunk.dims, perl = TRUE)) {
@@ -154,57 +399,60 @@ create <- function(
   } else {
     chunk.dims <- as.integer(x = chunk.dims)
   }
-  new.loom <- loom$new(filename = filename, mode = 'r')
-  h5attr(x = new.loom, which = 'version') <- as.character(x = packageVersion(pkg = 'loomR'))
+  new.loom <- loom$new(filename = filename, mode = 'w-')
+  # Create the matrix
   new.loom$create_dataset(
     name = 'matrix',
     robj = data,
     chunk_dims = chunk.dims
   )
+  new.loom$matrix <- new.loom[['matrix']]
+  new.loom$shape <- new.loom[['matrix']]$dims
   # Groups
   new.loom$create_group(name = 'layers')
   new.loom$create_group(name = 'row_attrs')
   new.loom$create_group(name = 'col_attrs')
+  # Check for the existance of gene or cell names
+  if (!is.null(x = colnames(x = data))) {
+    new.loom$add.row.attribute(attribute = list('gene_names' = colnames(x = data)))
+  }
+  if (!is.null(x = rownames(x = data))) {
+    new.loom$add.col.attribute(attribute = list('cell_names' = rownames(x = data)))
+  }
+  # Store some constants as HDF5 attributes
+  h5attr(x = new.loom, which = 'version') <- new.loom$version
+  h5attr(x = new.loom, which = 'chunks') <- paste0(
+    '(',
+    paste(new.loom[['matrix']]$chunk_dims, collapse = ', '),
+    ')'
+  )
   # Add layers
-  for (ly in layers) {
-    new.loom$add.layer(layer = ly)
+  if (!is.null(x = layers)) {
+    new.loom$add.layer(layer = layers)
   }
-  for (rw in row.attrs) {
-    new.loom$add.row.attribute(attribute = rw)
+  if (!is.null(x = gene.attrs)) {
+    new.loom$add.row.attribute(attribute = gene.attrs)
   }
-  for (cl in col.attrs) {
-    new.loom$add.col.attribute(attribute = cl)
+  if (!is.null(x = cell.attrs)) {
+    new.loom$add.col.attribute(attribute = cell.attrs)
   }
   # Set last bit of information
-  new.loom$shape <- ''
-  new.loom$chunksize <- ''
+  chunks <- new.loom[['matrix']]$chunk_dims
+  chunks <- gsub(pattern = '(', replacement = '', x = chunks, fixed = TRUE)
+  chunks <- gsub(pattern = ')', replacement = '', x = chunks, fixed = TRUE)
+  chunks <- unlist(x = strsplit(x = chunks, split = ','))
+  new.loom$chunksize <- as.integer(x = chunks)
+  # Return the connection
+  return(new.loom)
 }
 
-# #' @importFrom utils packageVersion
-# #'
-# setMethod(
-#   f = 'initialize',
-#   signature = 'loom',
-#   definition = function(.Object, name, mode = 'a') {
-#     .Object <- callNextMethod(
-#       .Object,
-#       name = name,
-#       mode = mode
-#     )
-#     validateLoom(object = .Object)
-#     #.Object@version <- packageVersion(pkg = 'loom')
-#     # .Object@filename <- name
-#     .Object@shape <- dim(.Object['/matrix'])
-#     return(.Object)
-#   }
-# )
-
-
 #' Validate a loom object
 #'
 #' @param object A loom object
 #'
-#' @return None, errors if object is an invalid loom object
+#' @return None, errors out if object is an invalid loom connection
+#'
+#' @seealso \code{\link{loom-class}}
 #'
 #' @export
 #'
@@ -220,8 +468,8 @@ validateLoom <- function(object) {
   }
   # There must be groups called '/col_attrs', '/row_attrs', and '/layers'
   required.groups <- c('row_attrs', 'col_attrs', 'layers')
-  dim.matrix <- object[[root.datasets]]$dims # Rows x Columns
-  names(dim.matrix) <- required.groups[c(2, 1)]
+  dim.matrix <- object[[root.datasets]]$dims # Columns x Rows
+  names(x = dim.matrix) <- required.groups[c(2, 1)]
   root.groups <- list.groups(object = object, path = '/', recursive = FALSE)
   group.msg <- paste0(
     "There can only be three groups in the loom file: '",
@@ -264,13 +512,32 @@ validateLoom <- function(object) {
 #'
 #' @return A loom file connection
 #'
+#' @seealso \code{\link{loom-class}}
+#'
 #' @export
 #'
 connect <- function(filename, mode = "r") {
+  if (!(mode %in% c('r', 'r+'))) {
+    stop("'mode' must be one of 'r' or 'r+'")
+  }
   new.loom <- loom$new(filename = filename, mode = mode)
   return(new.loom)
 }
 
+CreateLoomFromSeurat <- function(object, filename) {
+  object.data=t(object@raw.data[rownames(object@data),object@cell.names])
+  object.meta.data=object@meta.data
+  row_attrs=list(); col_attrs=list()
+  gene.names=colnames(object.data)
+  object.meta.data$ident = object@ident
+  object.meta.data$CellID = object@cell.names
+  for(i in 1:ncol(object.meta.data)) {
+    col_attrs[[colnames(object.meta.data)[i]]]=object.meta.data[,i]
+  }
+  row_attrs[["Gene"]]=gene.names
+  create(filename,object.data,gene.attrs = row_attrs, cell.attrs = col_attrs)
+}
+
 #need to comment
 #need to add progress bar
 #but otherwise, pretty cool