diff --git a/.Rbuildignore b/.Rbuildignore
index f7bd2a92..a31bf5a3 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -8,6 +8,7 @@
 ^doc$
 ^Meta$
 ^README\.Rmd$
+^README\.md$
 ^cran-comments\.md$
 ^.*\.loom$
 ^.*\.h5ad$
diff --git a/DESCRIPTION b/DESCRIPTION
index 71ba5959..7dcbefba 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,12 @@
 Package: SeuratObject
 Type: Package
 Title: Data Structures for Single Cell Data
-Version: 4.0.0
-Date: 2021-01-07
+Version: 4.0.1
+Date: 2021-05-07
 Authors@R: c(
   person(given = 'Rahul', family = 'Satija', email = 'rsatija@nygenome.org', role = 'aut', comment = c(ORCID = '0000-0001-9448-8833')),
   person(given = 'Andrew', family = 'Butler', email = 'abutler@nygenome.org', role = 'aut', comment = c(ORCID = '0000-0003-3608-0463')),
-  person(given = 'Paul', family = 'Hoffman', email = 'nygcSatijalab@nygenome.org', role = c('aut', 'cre'), comment = c(ORCID = '0000-0002-7693-8957')),
+  person(given = 'Paul', family = 'Hoffman', email = 'seurat@nygenome.org', role = c('aut', 'cre'), comment = c(ORCID = '0000-0002-7693-8957')),
   person(given = 'Tim', family = 'Stuart', email = 'tstuart@nygenome.org', role = 'aut', comment = c(ORCID = '0000-0002-3044-0897')),
   person(given = 'Jeff', family = 'Farrell', email = 'jfarrell@g.harvard.edu', role = 'ctb'),
   person(given = 'Shiwei', family = 'Zheng', email = 'szheng@nygenome.org', role = 'ctb', comment = c(ORCID = '0000-0001-6682-6743')),
@@ -17,10 +17,10 @@ Authors@R: c(
 Description: Defines S4 classes for single-cell genomic data and associated
   information, such as dimensionality reduction embeddings, nearest-neighbor
   graphs, and spatially-resolved coordinates. Provides data access methods and
-  R-native hooks to ensure the Seurat object is familiar to other R users. See 
-  Satija R, Farrell J, Gennert D, et al (2015) <doi:10.1038/nbt.3192>, 
-  Macosko E, Basu A, Satija R, et al (2015) <doi:10.1016/j.cell.2015.05.002>, 
-  and Stuart T, Butler A, et al (2019) <doi:10.1016/j.cell.2019.05.031> for 
+  R-native hooks to ensure the Seurat object is familiar to other R users. See
+  Satija R, Farrell J, Gennert D, et al (2015) <doi:10.1038/nbt.3192>,
+  Macosko E, Basu A, Satija R, et al (2015) <doi:10.1016/j.cell.2015.05.002>,
+  and Stuart T, Butler A, et al (2019) <doi:10.1016/j.cell.2019.05.031> for
   more details.
 URL: https://satijalab.org/seurat,
       https://github.com/mojaveazure/seurat-object
@@ -30,21 +30,20 @@ License: GPL-3
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.1.1
-Depends: 
-  R (>= 3.5.0)
-Imports: 
+Depends:
+  R (>= 4.0.0)
+Imports:
     grDevices,
     grid,
-    Matrix (>= 1.2.18),
+    Matrix (>= 1.3.3),
     methods,
     Rcpp (>= 1.0.5),
     rlang (>= 0.4.7),
-    sctransform,
     stats,
     tools,
     utils
-Suggests: 
-    tinytest
+Suggests:
+    testthat
 Collate: 'RcppExports.R' 'utils.R' 'zzz.R' 'generics.R'
       'assay.R' 'command.R' 'data.R' 'default.R'
       'jackstraw.R' 'dimreduc.R' 'graph.R' 'neighbor.R'
diff --git a/NAMESPACE b/NAMESPACE
index 77fc3697..e13d1596 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -103,6 +103,7 @@ S3method(Stdev,Seurat)
 S3method(Tool,Seurat)
 S3method(VariableFeatures,Assay)
 S3method(VariableFeatures,Seurat)
+S3method(Version,Seurat)
 S3method(WhichCells,Assay)
 S3method(WhichCells,Seurat)
 S3method(as.Graph,Matrix)
@@ -123,6 +124,8 @@ S3method(dimnames,Assay)
 S3method(dimnames,DimReduc)
 S3method(dimnames,Seurat)
 S3method(droplevels,Seurat)
+S3method(head,Assay)
+S3method(head,Seurat)
 S3method(length,DimReduc)
 S3method(levels,Seurat)
 S3method(merge,Assay)
@@ -135,6 +138,8 @@ S3method(subset,Assay)
 S3method(subset,DimReduc)
 S3method(subset,Seurat)
 S3method(subset,SpatialImage)
+S3method(tail,Assay)
+S3method(tail,Seurat)
 export("%iff%")
 export("%||%")
 export("DefaultAssay<-")
@@ -149,6 +154,7 @@ export("Tool<-")
 export("VariableFeatures<-")
 export(AddMetaData)
 export(Assays)
+export(AttachDeps)
 export(Cells)
 export(CellsByIdentities)
 export(CheckGC)
@@ -170,6 +176,7 @@ export(Images)
 export(Index)
 export(Indices)
 export(IsGlobal)
+export(IsMatrixEmpty)
 export(IsS4List)
 export(JS)
 export(Key)
@@ -196,6 +203,7 @@ export(Stdev)
 export(Tool)
 export(UpdateSeuratObject)
 export(VariableFeatures)
+export(Version)
 export(WhichCells)
 export(as.Graph)
 export(as.Neighbor)
@@ -245,7 +253,6 @@ importFrom(rlang,enquo)
 importFrom(rlang,eval_tidy)
 importFrom(rlang,is_null)
 importFrom(rlang,is_quosure)
-importFrom(sctransform,get_residuals)
 importFrom(stats,na.omit)
 importFrom(tools,file_path_sans_ext)
 importFrom(utils,.DollarNames)
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 00000000..b6093dcc
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +1,23 @@
+# SeuratObject 4.0.1
+
+## Added
+- `head` and `tail` methods for `Seurat` and `Assay` objects (#5)
+- New utility functions (#6):
+  - `AttachDeps` to attach required imported dependencies on package attachment
+  - `IsMatrixEmpty` to test if a matrix is empty or not
+
+## Changed
+- Allow super classes to replace child classes (#1). For example, allows `Assay`
+  objects to replace `Seurat::SCTAssay` or `Signac::ChromatinAssay` objects of
+  the same name
+- Better support for creating sparse matrices from `data.table`/`tibble` 
+  objects (#4)
+- Improved error messages for clashing object names (#7)
+- Allow returning a `NULL` if a subset results in zero cells (#9)
+
+## Removed
+- SCT-specific code (#2)
+
+# SeuratObject 4.0.0
+
+- Initial release of SeuratObject
diff --git a/R/assay.R b/R/assay.R
index cd87c3c1..25da9e28 100644
--- a/R/assay.R
+++ b/R/assay.R
@@ -67,6 +67,7 @@ Assay <- setClass(
 #' new object with a lower cutoff.
 #' @param min.features Include cells where at least this many features are
 #' detected.
+#' @param ... Arguments passed to \code{\link{as.sparse}}
 #'
 #' @return A \code{\link{Assay}} object
 #'
@@ -91,7 +92,8 @@ CreateAssayObject <- function(
   counts,
   data,
   min.cells = 0,
-  min.features = 0
+  min.features = 0,
+  ...
 ) {
   if (missing(x = counts) && missing(x = data)) {
     stop("Must provide either 'counts' or 'data'")
@@ -125,7 +127,7 @@ CreateAssayObject <- function(
       stop("No feature names (rownames) names present in the input matrix")
     }
     if (!inherits(x = counts, what = 'dgCMatrix')) {
-      counts <- as(object = as.matrix(x = counts), Class = 'dgCMatrix')
+      counts <- as.sparse(x = counts, ...)
     }
     # Filter based on min.features
     if (min.features > 0) {
@@ -387,24 +389,6 @@ Misc.Assay <- .Misc
 #'
 RenameCells.Assay <- function(object, new.names = NULL, ...) {
   CheckDots(...)
-  if (IsSCT(assay = object)) {
-    if (is.null(x = Misc(object = object, slot = 'vst.set'))) {
-      suppressWarnings(Misc(object = object, slot = "vst.out")$cells_step1 <- new.names)
-      suppressWarnings(rownames(x = Misc(object = object, slot = "vst.out")$cell_attr) <- new.names)
-    } else{
-      suppressWarnings(
-        Misc(object, slot = "vst.set") <- lapply(
-          X = Misc(object = object, slot = "vst.set"),
-          FUN = function(x) {
-            new.names.vst <- new.names[which(x = x$cells_step1 %in% Cells(x = object))]
-            x$cells_step1 <- new.names.vst
-            rownames(x = x$cell_attr) <- new.names.vst
-            return(x)
-          }
-        )
-      )
-    }
-  }
   for (data.slot in c("counts", "data", "scale.data")) {
     old.data <- GetAssayData(object = object, slot = data.slot)
     if (ncol(x = old.data) <= 1) {
@@ -777,6 +761,17 @@ dimnames.Assay <- function(x) {
   return(dimnames(x = GetAssayData(object = x)))
 }
 
+#' @describeIn Assay-methods Get the first rows of feature-level metadata
+#'
+#' @inheritParams utils::head
+#'
+#' @return \code{head}: The first \code{n} rows of feature-level metadata
+#'
+#' @export
+#' @method head Assay
+#'
+head.Assay <- .head
+
 #' @describeIn Assay-methods Merge \code{Assay} objects
 #'
 #' @param y A vector or list of one or more objects to merge
@@ -839,43 +834,6 @@ merge.Assay <- function(
       new.data = merged.data
     )
   }
-  # merge SCT assay misc vst info and scale.data
-  if (all(IsSCT(assay = assays))) {
-    vst.set.new <- list()
-    idx <- 1
-    umi.assay.new <- list()
-    for (i in 1:length(x = assays)) {
-      vst.set.old <- Misc(object = assays[[i]], slot = "vst.set")
-      umi.assay.old <- Misc(object = assays[[i]], slot = "umi.assay")
-      if (!is.null(x = vst.set.old) && length(x = vst.set.old) > 1) {
-        for (j in 1:length(x = vst.set.old)) {
-          vst.set.new[[idx]] <- vst.set.old[[j]]
-          umi.assay.new[[idx]] <- umi.assay.old[[j]]
-          idx <- idx + 1
-        }
-      } else if (!is.null(x = Misc(object = assays[[i]], slot = "vst.out"))) {
-        vst.set.new[[idx]] <- Misc(object = assays[[i]], slot = "vst.out")
-        umi.assay.new[[idx]] <- Misc(object = assays[[i]], slot = "umi.assay")
-        idx <- idx + 1
-      }
-    }
-    Misc(object = combined.assay, slot = "vst.set") <- vst.set.new
-    Misc(object = combined.assay, slot = "umi.assay") <- umi.assay.new
-    scale.data <- do.call(
-      what = cbind,
-      args = lapply(
-        X = assays,
-        FUN = GetAssayData,
-        slot = 'scale.data'
-        # FUN = function(x) GetAssayData(object = x, slot = "scale.data")
-      )
-    )
-    combined.assay <- SetAssayData(
-      object = combined.assay,
-      slot = "scale.data",
-      new.data = scale.data
-    )
-  }
   return(combined.assay)
 }
 
@@ -935,33 +893,20 @@ subset.Assay <- function(x, cells = NULL, features = NULL, ...) {
   }
   VariableFeatures(object = x) <- VariableFeatures(object = x)[VariableFeatures(object = x) %in% features]
   slot(object = x, name = 'meta.features') <- x[[]][features, , drop = FALSE]
-  if (IsSCT(assay = x)) {
-    # subset cells and genes in the SCT assay
-    obj.misc <- Misc(object = x)
-    if ("vst.set" %in% names(x = obj.misc)) {
-      # set of vst.out objects
-      vst.info <- obj.misc[["vst.set"]]
-      for (i in seq_along(along.with = vst.info)) {
-        vst.info[[i]] <- SubsetVST(
-          sct.info = vst.info[[i]],
-          cells = cells,
-          features = features
-        )
-      }
-      obj.misc[["vst.set"]] <- vst.info
-    } else {
-      # just one vst.out
-      obj.misc[["vst.out"]] <- SubsetVST(
-        sct.info = obj.misc[["vst.out"]],
-        cells = cells,
-        features = features
-      )
-    }
-    slot(object = x, name = "misc") <- obj.misc
-  }
   return(x)
 }
 
+#' @describeIn Assay-methods Get the last rows of feature-level metadata
+#'
+#' @return \code{tail}: The last \code{n} rows of feature-level metadata
+#'
+#' @importFrom utils tail
+#'
+#' @export
+#' @method tail Assay
+#'
+tail.Assay <- .tail
+
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # S4 methods
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -1188,31 +1133,11 @@ CalcN <- function(object) {
     return(NULL)
   }
   return(list(
-    nCount = colSums(x = object, slot = 'counts'),
-    nFeature = colSums(x = GetAssayData(object = object, slot = 'counts') > 0)
+    nCount = Matrix::colSums(x = object, slot = 'counts'),
+    nFeature = Matrix::colSums(x = GetAssayData(object = object, slot = 'counts') > 0)
   ))
 }
 
-#' Check whether an assay has been processed by sctransform
-#'
-#' @param assay assay to check
-#'
-#' @return Boolean
-#'
-#' @keywords internal
-#'
-#' @noRd
-#'
-IsSCT <- function(assay) {
-  if (is.list(x = assay)) {
-    sct.check <- lapply(X = assay, FUN = function(x) {
-      return(!is.null(x = Misc(object = x, slot = 'vst.out')) | !is.null(x = Misc(object = x, slot = 'vst.set')))
-    })
-    return(unlist(x = sct.check))
-  }
-  return(!is.null(x = Misc(object = assay, slot = 'vst.out')) | !is.null(x = Misc(object = assay, slot = 'vst.set')))
-}
-
 #' Subset cells in vst data
 #'
 #' @param sct.info A vst.out list
diff --git a/R/generics.R b/R/generics.R
index fe38023c..42081a58 100644
--- a/R/generics.R
+++ b/R/generics.R
@@ -827,6 +827,23 @@ VariableFeatures <- function(object, selection.method = NULL, ...) {
   UseMethod(generic = 'VariableFeatures<-', object = object)
 }
 
+#' Get Version Information
+#'
+#' @param object An object
+#' @param ... Arguments passed to other methods
+#'
+#' @rdname Version
+#' @export Version
+#'
+#' @concept data-access
+#'
+#' @examples
+#' Version(pbmc_small)
+#'
+Version <- function(object, ...) {
+  UseMethod(generic = "Version", object = object)
+}
+
 #' Identify cells matching certain criteria
 #'
 #' Returns a list of cells that match a particular set of criteria such as
diff --git a/R/seurat.R b/R/seurat.R
index c8d492ef..142c2eab 100644
--- a/R/seurat.R
+++ b/R/seurat.R
@@ -189,6 +189,8 @@ Assays <- function(object, slot = NULL) {
 #' @param idents A vector of identity class levels to limit resulting list to;
 #' defaults to all identity class levels
 #' @param cells A vector of cells to grouping to
+#' @param return.null If no cells are request, return a \code{NULL};
+#' by default, throws an error
 #'
 #' @return A named list where names are identity classes and values are vectors
 #' of cells belonging to that class
@@ -200,10 +202,18 @@ Assays <- function(object, slot = NULL) {
 #' @examples
 #' CellsByIdentities(object = pbmc_small)
 #'
-CellsByIdentities <- function(object, idents = NULL, cells = NULL) {
+CellsByIdentities <- function(
+  object,
+  idents = NULL,
+  cells = NULL,
+  return.null = FALSE
+) {
   cells <- cells %||% colnames(x = object)
   cells <- intersect(x = cells, y = colnames(x = object))
   if (length(x = cells) == 0) {
+    if (isTRUE(x = return.null)) {
+      return(NULL)
+    }
     stop("Cannot find cells provided")
   }
   idents <- idents %||% levels(x = object)
@@ -347,6 +357,17 @@ FetchData <- function(object, vars, cells = NULL, slot = 'data') {
   # Pull vars from object metadata
   meta.vars <- vars[vars %in% colnames(x = object[[]]) & !(vars %in% names(x = data.fetched))]
   data.fetched <- c(data.fetched, object[[meta.vars]][cells, , drop = FALSE])
+  meta.default <- meta.vars[meta.vars %in% rownames(x = GetAssayData(object = object, slot = slot))]
+  if (length(x = meta.default)) {
+    warning(
+      "The following variables were found in both object metadata and the default assay: ",
+      paste0(meta.default, collapse = ", "),
+      "\nReturning metadata; if you want the feature, please use the assay's key (eg. ",
+      paste0(Key(object = object[[DefaultAssay(object = object)]]), meta.default[1]),
+      ")",
+      call. = FALSE
+    )
+  }
   # Pull vars from the default assay
   default.vars <- vars[vars %in% rownames(x = GetAssayData(object = object, slot = slot)) & !(vars %in% names(x = data.fetched))]
   data.fetched <- c(
@@ -909,6 +930,10 @@ Command.Seurat <- function(object, command = NULL, value = NULL, ...) {
   return(params[[value]])
 }
 
+#' @param row.names When \code{counts} is a \code{data.frame} or
+#' \code{data.frame}-derived object: an optional vector of feature names to be
+#' used
+#'
 #' @rdname CreateSeuratObject
 #' @method CreateSeuratObject default
 #' @export
@@ -922,6 +947,7 @@ CreateSeuratObject.default <- function(
   meta.data = NULL,
   min.cells = 0,
   min.features = 0,
+  row.names = NULL,
   ...
 ) {
   if (!is.null(x = meta.data)) {
@@ -932,7 +958,8 @@ CreateSeuratObject.default <- function(
   assay.data <- CreateAssayObject(
     counts = counts,
     min.cells = min.cells,
-    min.features = min.features
+    min.features = min.features,
+    row.names = row.names
   )
   if (!is.null(x = meta.data)) {
     common.cells <- intersect(
@@ -1763,10 +1790,11 @@ VariableFeatures.Seurat <- function(
 
 #' @param idents A vector of identity classes to keep
 #' @param slot Slot to pull feature data for
-#' @param downsample Maximum number of cells per identity class, default is \code{Inf};
-#' downsampling will happen after all other operations, including inverting the
-#' cell selection
+#' @param downsample Maximum number of cells per identity class, default is
+#' \code{Inf}; downsampling will happen after all other operations, including
+#' inverting the cell selection
 #' @param seed Random seed for downsampling. If NULL, does not set a seed
+#' @inheritDotParams CellsByIdentities
 #'
 #' @importFrom stats na.omit
 #' @importFrom rlang is_quosure enquo eval_tidy
@@ -1786,7 +1814,7 @@ WhichCells.Seurat <- function(
   seed = 1,
   ...
 ) {
-  CheckDots(...)
+  CheckDots(..., fxns = CellsByIdentities)
   if (!is.null(x = seed)) {
     set.seed(seed = seed)
   }
@@ -1861,11 +1889,11 @@ WhichCells.Seurat <- function(
     )
     cells <- rownames(x = data.subset)[eval_tidy(expr = expr, data = data.subset)]
   }
-  if (invert) {
+  if (isTRUE(x = invert)) {
     cell.order <- colnames(x = object)
     cells <- colnames(x = object)[!colnames(x = object) %in% cells]
   }
-  cells <- CellsByIdentities(object = object, cells = cells)
+  cells <- CellsByIdentities(object = object, cells = cells, ...)
   cells <- lapply(
     X = cells,
     FUN = function(x) {
@@ -1880,6 +1908,15 @@ WhichCells.Seurat <- function(
   return(cells)
 }
 
+#' @rdname Version
+#' @method Version Seurat
+#' @export
+#'
+Version.Seurat <- function(object, ...) {
+  CheckDots(...)
+  return(slot(object = object, name = 'version'))
+}
+
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # Methods for R-defined generics
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -1902,6 +1939,7 @@ WhichCells.Seurat <- function(
 #'  }
 #' }
 #' @param j,cells Cell names or indices
+#' @param n The number of rows of metadata to return
 #' @param ... Arguments passed to other methods
 #'
 #' @name Seurat-methods
@@ -2122,6 +2160,21 @@ droplevels.Seurat <- function(x, ...) {
   return(x)
 }
 
+#' @describeIn Seurat-methods Get the first rows of cell-level metadata
+#'
+#' @return \code{head}: The first \code{n} rows of cell-level metadata
+#'
+#' @importFrom utils head
+#'
+#' @export
+#' @method head Seurat
+#'
+#' @examples
+#' # Get the first 10 rows of cell-level metadata
+#' head(pbmc_small)
+#'
+head.Seurat <- .head
+
 #' @rdname Idents
 #' @export
 #' @method levels Seurat
@@ -2245,39 +2298,6 @@ merge.Seurat <- function(
         ))
       }
     )
-    if (all(IsSCT(assay = assays.merge))) {
-      scaled.features <- unique(x = unlist(x = lapply(
-        X = assays.merge,
-        FUN = function(x) rownames(x = GetAssayData(object = x, slot = "scale.data")))
-      ))
-      for (ob in 1:length(x = objects)) {
-        if (assay %in% FilterObjects(object = objects[[ob]], classes.keep = "Assay")) {
-          objects[[ob]] <- suppressWarnings(GetResidual(object = objects[[ob]], features = scaled.features, assay = assay, verbose = FALSE))
-          assays.merge[[ob]] <- objects[[ob]][[assay]]
-        }
-      }
-      # handle case where some features aren't in counts and can't be retrieved with
-      # GetResidual - take intersection
-      scaled.features <- names(x = which(x = table(x = unlist(x = lapply(
-        X = assays.merge,
-        FUN = function(x) rownames(x = GetAssayData(object = x, slot = "scale.data")))
-      )) == length(x = assays.merge)))
-      if (length(x = scaled.features) > 0) {
-        for (a in 1:length(x = assays.merge)) {
-          assays.merge[[a]] <- SetAssayData(
-            object = assays.merge[[a]],
-            slot = "scale.data",
-            new.data = GetAssayData(object = assays.merge[[a]], slot = "scale.data")[scaled.features, ])
-        }
-      } else {
-        for (a in 1:length(x = assays.merge)) {
-          assays.merge[[a]] <- SetAssayData(
-            object = assays.merge[[a]],
-            slot = "scale.data",
-            new.data = new(Class = "matrix"))
-        }
-      }
-    }
     merged.assay <- merge(
       x = assays.merge[[1]],
       y = assays.merge[2:length(x = assays.merge)],
@@ -2398,6 +2418,7 @@ names.Seurat <- function(x) {
 
 #' @describeIn Seurat-methods Subset a \code{\link{Seurat}} object
 #'
+#' @inheritParams CellsByIdentities
 #' @param subset Logical expression indicating features/variables to keep
 #' @param idents A vector of identity classes to keep
 #'
@@ -2425,6 +2446,7 @@ subset.Seurat <- function(
   cells = NULL,
   features = NULL,
   idents = NULL,
+  return.null = FALSE,
   ...
 ) {
   x <- UpdateSlots(object = x)
@@ -2436,9 +2458,13 @@ subset.Seurat <- function(
     cells = cells,
     idents = idents,
     expression = subset,
+    return.null = TRUE,
     ...
   )
   if (length(x = cells) == 0) {
+    if (isTRUE(x = return.null)) {
+      return(NULL)
+    }
     stop("No cells found", call. = FALSE)
   }
   if (all(cells %in% Cells(x = x)) && length(x = cells) == length(x = Cells(x = x)) && is.null(x = features)) {
@@ -2502,6 +2528,21 @@ subset.Seurat <- function(
   return(x)
 }
 
+#' @describeIn Seurat-methods Get the last rows of cell-level metadata
+#'
+#' @return \code{tail}: The last \code{n} rows of cell-level metadata
+#'
+#' @importFrom utils tail
+#'
+#' @export
+#' @method tail Seurat
+#'
+#' @examples
+#' # Get the last 10 rows of cell-level metadata
+#' tail(pbmc_small)
+#'
+tail.Seurat <- .tail
+
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # S4 methods
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -2697,7 +2738,10 @@ setMethod( # because R doesn't allow S3-style [[<- for S4 classes
         stop("All cells in the object being added must match the cells in this object", call. = FALSE)
       }
       # Ensure we're not duplicating object names
-      if (!is.null(x = FindObject(object = x, name = i)) && !inherits(x = value, what = c(class(x = x[[i]]), 'NULL'))) {
+      duplicate <- !is.null(x = FindObject(object = x, name = i)) &&
+        !inherits(x = value, what = c(class(x = x[[i]]), 'NULL')) &&
+        !inherits(x = x[[i]], what = class(x = value))
+      if (isTRUE(x = duplicate)) {
         stop(
           "This object already contains ",
           i,
@@ -2708,7 +2752,10 @@ setMethod( # because R doesn't allow S3-style [[<- for S4 classes
             no = ' '
           ),
           class(x = x[[i]]),
-          "; duplicate names are not allowed",
+          ", so ",
+          i,
+          " cannot be used for a ",
+          class(x = value),
           call. = FALSE
         )
       }
@@ -3118,259 +3165,6 @@ FindObject <- function(object, name) {
   return(NULL)
 }
 
-#' Calculate pearson residuals of features not in the scale.data
-#'
-#' This function calls sctransform::get_residuals.
-#'
-#' @param object A seurat object
-#' @param features Name of features to add into the scale.data
-#' @param assay Name of the assay of the seurat object generated by SCTransform
-#' @param umi.assay Name of the assay of the seurat object containing UMI matrix and the default is
-#' RNA
-#' @param clip.range Numeric of length two specifying the min and max values the Pearson residual
-#' will be clipped to
-#' @param replace.value Recalculate residuals for all features, even if they are already present.
-#' Useful if you want to change the clip.range.
-#' @param verbose Whether to print messages and progress bars
-#'
-#' @return Returns a Seurat object containing  pearson residuals of added features in its scale.data
-#'
-#' @importFrom sctransform get_residuals
-#'
-#' @keywords internal
-#'
-#' @seealso \code{\link[sctransform]{get_residuals}}
-#'
-#' @noRd
-#'
-#' @examples
-#' \dontrun{
-#' pbmc_small <- SCTransform(pbmc_small, variable.features.n = 20)
-#' pbmc_small <- GetResidual(pbmc_small, features = c('MS4A1', 'TCL1A'))
-#' }
-#'
-GetResidual <- function(
-  object,
-  features,
-  assay = "SCT",
-  umi.assay = NULL,
-  clip.range = NULL,
-  replace.value = FALSE,
-  verbose = TRUE
-) {
-  if (!IsSCT(assay = object[[assay]])) {
-    stop(assay, " assay was not generated by SCTransform")
-  }
-  model.name <- ifelse(
-    test = "vst.set" %in% names(x = Misc(object = object[[assay]])),
-    yes = "vst.set",
-    no = "vst.out"
-  )
-  if (length(x = Misc(object = object[[assay]])[[model.name]]) == 0) {
-    warning("SCT model not present in assay")
-    return(object)
-  }
-  umi.assay <- umi.assay %||% Misc(object = object[[assay]], slot = "umi.assay")
-  umi.assay <- umi.assay %||% "RNA" # for object created in 3.1.1 or earlier, default to RNA
-  if (replace.value) {
-    new_features <- features
-  } else {
-    new_features <- setdiff(
-      x = features,
-      y = rownames(x = GetAssayData(object = object, assay = assay, slot = "scale.data"))
-    )
-  }
-  if (length(x = new_features) == 0) {
-    if (verbose) {
-      message("Pearson residuals of input features exist already")
-    }
-  } else {
-    if (is.null(x = Misc(object = object[[assay]], slot = 'vst.set'))) {
-      vst_out <- Misc(object = object[[assay]], slot = 'vst.out')
-      # filter cells not in the object but in the SCT model
-      vst_out$cell_attr <- vst_out$cell_attr[Cells(x = object), ]
-      vst_out$cells_step1 <- intersect(x = vst_out$cells_step1, y = Cells(x = object))
-      object <- GetResidualVstOut(
-        object = object,
-        assay = assay,
-        umi.assay = umi.assay,
-        new_features = new_features,
-        vst_out = vst_out,
-        clip.range = clip.range,
-        verbose = verbose
-      )
-    } else {
-      # Calculate Pearson Residual from integrated object SCT assay
-      vst.set <- Misc(object = object[[assay]], slot = 'vst.set')
-      scale.data <- GetAssayData(
-        object = object,
-        assay = assay,
-        slot = "scale.data"
-      )
-      vst_set_genes <-  sapply(1:length(vst.set), function(x) rownames(vst.set[[x]]$model_pars_fit))
-      vst_set_genes <- Reduce(intersect, vst_set_genes)
-      diff_features <- setdiff(
-        x = new_features,
-        y = vst_set_genes
-      )
-      if (length(x = diff_features) != 0) {
-        warning(
-          "The following ", length(x = diff_features),
-          " features do not exist in all SCT models: ",
-          paste(diff_features, collapse = " ")
-        )
-      }
-      new_features <- intersect(
-        x = new_features,
-        y = vst_set_genes
-      )
-      if (length(new_features) != 0) {
-        object <- SetAssayData(
-          object = object,
-          assay = assay,
-          slot = "scale.data",
-          new.data = scale.data[!rownames(x = scale.data) %in% new_features, , drop = FALSE]
-        )
-        new.scale.data <- matrix(nrow = length(new_features), ncol = 0)
-        rownames(x = new.scale.data) <- new_features
-        for (v in 1:length(x = vst.set)) {
-          vst_out <- vst.set[[v]]
-          # confirm that cells from SCT model also exist in the integrated object
-          cells.v <- intersect(x = rownames(x = vst_out$cell_attr), y = Cells(x = object))
-          if (length(x = cells.v) != 0) {
-            vst_out$cell_attr <- vst_out$cell_attr[cells.v, ]
-            vst_out$cells_step1 <- intersect(x = vst_out$cells_step1, y = cells.v)
-            object.v <- subset(x = object, cells = cells.v)
-            object.v <- GetResidualVstOut(
-              object = object.v,
-              assay = assay,
-              umi.assay = umi.assay[[v]],
-              new_features = new_features,
-              vst_out = vst_out,
-              clip.range = clip.range,
-              verbose = verbose
-            )
-            new.scale.data <- cbind(
-              new.scale.data,
-              GetAssayData(
-                object = object.v,
-                assay = assay,
-                slot = "scale.data")[new_features, , drop = FALSE]
-            )
-          }
-        }
-        object <- SetAssayData(
-          object = object,
-          assay = assay,
-          slot = "scale.data",
-          new.data = rbind(
-            GetAssayData(object = object, slot = 'scale.data', assay = assay),
-            new.scale.data
-          )
-        )
-      }
-    }
-  }
-  return(object)
-}
-
-# Calculate pearson residuals of features not in the scale.data
-# This function is the secondary function under GetResidual
-#
-# @param object A seurat object
-# @param features Name of features to add into the scale.data
-# @param assay Name of the assay of the seurat object generated by SCTransform
-# @param vst_out The SCT parameter list
-# @param clip.range Numeric of length two specifying the min and max values the Pearson residual
-# will be clipped to
-# Useful if you want to change the clip.range.
-# @param verbose Whether to print messages and progress bars
-#
-# @return Returns a Seurat object containing  pearson residuals of added features in its scale.data
-#
-#' @importFrom sctransform get_residuals
-#
-GetResidualVstOut <- function(
-  object,
-  assay,
-  umi.assay,
-  new_features,
-  vst_out,
-  clip.range,
-  verbose
-) {
-  diff_features <- setdiff(
-    x = new_features,
-    y = rownames(x = vst_out$model_pars_fit)
-  )
-  intersect_feature <- intersect(
-    x = new_features,
-    y = rownames(x = vst_out$model_pars_fit)
-  )
-  if (length(x = diff_features) == 0) {
-    umi <- GetAssayData(object = object, assay = umi.assay, slot = "counts" )[new_features, , drop = FALSE]
-  } else {
-    warning(
-      "The following ", length(x = diff_features),
-      " features do not exist in the counts slot: ",
-      paste(diff_features, collapse = " ")
-    )
-    if (length(x = intersect_feature) == 0) {
-      return(object)
-    }
-    umi <- GetAssayData(object = object, assay = umi.assay, slot = "counts" )[intersect_feature, , drop = FALSE]
-  }
-  if (is.null(x = clip.range)) {
-    if (length(vst_out$arguments$sct.clip.range) != 0 ) {
-      clip.max <- max(vst_out$arguments$sct.clip.range)
-      clip.min <- min(vst_out$arguments$sct.clip.range)
-    } else {
-      clip.max <- max(vst_out$arguments$res_clip_range)
-      clip.min <- min(vst_out$arguments$res_clip_range)
-    }
-  } else {
-    clip.max <- max(clip.range)
-    clip.min <- min(clip.range)
-  }
-  new_residual <- get_residuals(
-    vst_out = vst_out,
-    umi = umi,
-    residual_type = "pearson",
-    res_clip_range = c(clip.min, clip.max),
-    verbosity = as.numeric(x = verbose) * 2
-  )
-  new_residual <- as.matrix(x = new_residual)
-  # centered data
-  new_residual <- new_residual - rowMeans(new_residual)
-  # remove genes from the scale.data if genes are part of new_features
-  scale.data <- GetAssayData(object = object, assay = assay, slot = "scale.data")
-  object <- SetAssayData(
-    object = object,
-    assay = assay,
-    slot = "scale.data",
-    new.data = scale.data[!rownames(x = scale.data) %in% new_features, , drop = FALSE]
-  )
-  if (nrow(x = GetAssayData(object = object, slot = 'scale.data', assay = assay)) == 0 ) {
-    object <- SetAssayData(
-      object = object,
-      slot = 'scale.data',
-      new.data = new_residual,
-      assay = assay
-    )
-  } else {
-    object <- SetAssayData(
-      object = object,
-      slot = 'scale.data',
-      new.data = rbind(
-        GetAssayData(object = object, slot = 'scale.data', assay = assay),
-        new_residual
-      ),
-      assay = assay
-    )
-  }
-  return(object)
-}
-
 #' Update Seurat v2 Internal Objects
 #'
 #' Helper functions to update old Seurat v2 objects to v3/v4 objects
diff --git a/R/utils.R b/R/utils.R
index 1b1c9082..198f0411 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -106,6 +106,37 @@ rlang::`%||%`
   return(x)
 }
 
+#' Attach Required Packages
+#'
+#' Helper function to attach required packages. Detects if a package is already
+#' attached and if so, skips it. Should be called in \code{\link[base]{.onAttach}}
+#'
+#' @param deps A character vector of packages to attach
+#'
+#' @return Invisibly returns \code{NULL}
+#'
+#' @export
+#'
+#' @concept utils
+#'
+#' @examples
+#' # Use in your .onAttach hook
+#' if (FALSE) {
+#'   .onAttach <- function(libname, pkgname) {
+#'     AttachDeps(c("SeuratObject", "rlang"))
+#'   }
+#' }
+#'
+AttachDeps <- function(deps) {
+  for (d in deps) {
+    if (!paste0('package:', d) %in% search()) {
+      packageStartupMessage("Attaching ", d)
+      attachNamespace(ns = d)
+    }
+  }
+  return(invisible(x = NULL))
+}
+
 #' Conditional Garbage Collection
 #'
 #' Call \code{gc} only when desired
@@ -125,6 +156,29 @@ CheckGC <- function(option = 'SeuratObject.memsafe') {
   return(invisible(x = NULL))
 }
 
+#' Check if a matrix is empty
+#'
+#' Takes a matrix and asks if it's empty (either 0x0 or 1x1 with a value of NA)
+#'
+#' @param x A matrix
+#'
+#' @return Whether or not \code{x} is empty
+#'
+#' @export
+#'
+#' @concept utils
+#'
+#' @examples
+#' IsMatrixEmpty(new("matrix"))
+#' IsMatrixEmpty(matrix())
+#' IsMatrixEmpty(matrix(1:3))
+#'
+IsMatrixEmpty <- function(x) {
+  matrix.dims <- dim(x = x)
+  matrix.na <- all(matrix.dims == 1) && all(is.na(x = x))
+  return(all(matrix.dims == 0) || matrix.na)
+}
+
 #' @name s4list
 #' @rdname s4list
 #'
@@ -232,12 +286,21 @@ RowMergeSparseMatrices <- function(mat1, mat2) {
 # Methods for Seurat-defined generics
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+#' @param row.names \code{NULL} or a character vector giving the row names for
+#' the data; missing values are not allowed
+#'
 #' @rdname as.sparse
 #' @export
 #' @method as.sparse data.frame
 #'
-as.sparse.data.frame <- function(x, ...) {
+as.sparse.data.frame <- function(x, row.names = NULL, ...) {
   CheckDots(...)
+  dnames <- list(row.names %||% rownames(x = x), colnames(x = x))
+  if (length(x = dnames[[1]]) != nrow(x = x)) {
+    stop("Differing numbers of rownames and rows", call. = FALSE)
+  }
+  x <- as.data.frame(x = x)
+  dimnames(x = x) <- dnames
   return(as.sparse(x = as.matrix(x = x)))
 }
 
@@ -522,31 +585,6 @@ ExtractField <- function(string, field = 1, delim = "_") {
   ))
 }
 
-#' Check if a matrix is empty
-#'
-#' Takes a matrix and asks if it's empty (either 0x0 or 1x1 with a value of NA)
-#'
-#' @param x A matrix
-#'
-#' @return Whether or not \code{x} is empty
-#'
-#' @keywords internal
-#'
-#' @noRd
-#'
-#' @examples
-#' \donttest{
-#' SeuratObject:::IsMatrixEmpty(new("matrix"))
-#' SeuratObject:::IsMatrixEmpty(matrix())
-#' SeuratObject:::IsMatrixEmpty(matrix(1:3))
-#' }
-#'
-IsMatrixEmpty <- function(x) {
-  matrix.dims <- dim(x = x)
-  matrix.na <- all(matrix.dims == 1) && all(is.na(x = x))
-  return(all(matrix.dims == 0) || matrix.na)
-}
-
 #' Test Null Pointers
 #'
 #' Check to see if a C++ pointer is a null pointer on the compiled side
diff --git a/R/zzz.R b/R/zzz.R
index 4c6c0d32..c3571909 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -48,7 +48,7 @@ setOldClass(Classes = 'package_version')
 # Internal
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-#' Object Metadata
+#' Add Object Metadata
 #'
 #' Internal \code{\link{AddMetaData}} definition
 #'
@@ -77,6 +77,28 @@ setOldClass(Classes = 'package_version')
   return(object)
 }
 
+#' Head and Tail Object Metadata
+#'
+#' Internal \code{\link[utils]{head}} and \code{\link[utils]{tail}} definitions
+#'
+#' @param x An object
+#' @param n Number of rows to return
+#' @inheritDotParams utils::head
+#'
+#' @return The first or last \code{n} rows of object metadata
+#'
+#' @keywords internal
+#'
+#' @noRd
+#'
+.head <- function(x, n = 10L, ...) {
+  return(head(x = x[[]], n = n, ...))
+}
+
+.tail <- function(x, n = 10L, ...) {
+  return(tail(x = x[[]], n = n, ...))
+}
+
 #' Miscellaneous Data
 #'
 #' Internal functions for getting and setting miscellaneous data
diff --git a/cran-comments.md b/cran-comments.md
index b42113b4..9c573c69 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,18 +1,26 @@
-# SeuratObject v4.0.0
+# SeuratObject v4.0.1
 
 ## Test environments
-* local ubuntu 20.04 install, R 4.0.3
-* local macOS 10.14.6 install, R 4.0.2
-* win-builder (oldrelease, release, devel)
+* local ubuntu 20.04 install, R 4.0.5
+* local macOS 10.14.6 install, R 4.0.5
+* win-builder (release, devel)
 
 ## R CMD check results
 
-There were no ERRORs or WARNINGs
+There were no ERRORs, WARNINGs
 
-There was 1 NOTE:
+There was one NOTE:
 
-* This is a new release.
+* checking CRAN incoming feasibility ... NOTE
+Maintainer: ‘Paul Hoffman <seurat@nygenome.org>’
+
+New maintainer:
+  Paul Hoffman <seurat@nygenome.org>
+Old maintainer(s):
+  Paul Hoffman <nygcSatijalab@nygenome.org>
+
+We are switching email addresses for the maintainer
 
 ## Downstream dependencies
 
-There are no packages that import or depend on SeuratObject
+There are two packages that import SeuratObject; this update does impact their functionality. There are two packages that suggest SeuratObject; this update does not impact their functionality.
diff --git a/man/Assay-methods.Rd b/man/Assay-methods.Rd
index 3d4192c4..66613da4 100644
--- a/man/Assay-methods.Rd
+++ b/man/Assay-methods.Rd
@@ -6,8 +6,10 @@
 \alias{[[.Assay}
 \alias{dim.Assay}
 \alias{dimnames.Assay}
+\alias{head.Assay}
 \alias{merge.Assay}
 \alias{subset.Assay}
+\alias{tail.Assay}
 \alias{[[<-,Assay-method}
 \alias{colMeans,Assay-method}
 \alias{colSums,Assay-method}
@@ -24,10 +26,14 @@
 
 \method{dimnames}{Assay}(x)
 
+\method{head}{Assay}(x, n = 10L, ...)
+
 \method{merge}{Assay}(x = NULL, y = NULL, add.cell.ids = NULL, merge.data = TRUE, ...)
 
 \method{subset}{Assay}(x, cells = NULL, features = NULL, ...)
 
+\method{tail}{Assay}(x, n = 10L, ...)
+
 \S4method{[[}{Assay}(x, i, j, ...) <- value
 
 \S4method{colMeans}{Assay}(x, na.rm = FALSE, dims = 1, ..., slot = "data")
@@ -52,6 +58,16 @@ feature names or indices}
 
 \item{drop}{See \code{\link[base]{drop}}}
 
+\item{n}{an integer vector of length up to \code{dim(x)} (or 1,
+    for non-dimensioned objects). Values specify the indices to be
+    selected in the corresponding dimension (or along the length) of the
+    object. A positive value of \code{n[i]} includes the first/last
+    \code{n[i]} indices in that dimension, while a negative value
+    excludes the last/first \code{abs(n[i])}, including all remaining
+    indices. \code{NA} or non-specified values (when \code{length(n) <
+      length(dim(x))}) select all indices in that dimension. Must
+    contain at least one non-missing value.}
+
 \item{y}{A vector or list of one or more objects to merge}
 
 \item{add.cell.ids}{A character vector of \code{length(x = c(x, y))};
@@ -82,10 +98,14 @@ means/sums on}
 
 \code{dimnames}: Feature (row) and cell (column) names
 
+\code{head}: The first \code{n} rows of feature-level metadata
+
 \code{merge}: Merged object
 
 \code{subset}: A subsetted \code{Assay}
 
+\code{tail}: The last \code{n} rows of feature-level metadata
+
 \code{[[<-}: \code{x} with metadata \code{value} added as \code{i}
 
 \code{colMeans}: The column (cell-wise) means of \code{slot}
@@ -113,10 +133,14 @@ other packages
 
 \item \code{dimnames.Assay}: Cell- and feature-names for an \code{Assay}
 
+\item \code{head.Assay}: Get the first rows of feature-level metadata
+
 \item \code{merge.Assay}: Merge \code{Assay} objects
 
 \item \code{subset.Assay}: Subset an \code{Assay}
 
+\item \code{tail.Assay}: Get the last rows of feature-level metadata
+
 \item \code{[[<-,Assay-method}: Add feature-level metadata
 
 \item \code{colMeans,Assay-method}: Calculate \code{\link[base]{colMeans}} on an
diff --git a/man/AttachDeps.Rd b/man/AttachDeps.Rd
new file mode 100644
index 00000000..65028afa
--- /dev/null
+++ b/man/AttachDeps.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{AttachDeps}
+\alias{AttachDeps}
+\title{Attach Required Packages}
+\usage{
+AttachDeps(deps)
+}
+\arguments{
+\item{deps}{A character vector of packages to attach}
+}
+\value{
+Invisibly returns \code{NULL}
+}
+\description{
+Helper function to attach required packages. Detects if a package is already
+attached and if so, skips it. Should be called in \code{\link[base]{.onAttach}}
+}
+\examples{
+# Use in your .onAttach hook
+if (FALSE) {
+  .onAttach <- function(libname, pkgname) {
+    AttachDeps(c("SeuratObject", "rlang"))
+  }
+}
+
+}
+\concept{utils}
diff --git a/man/CellsByIdentities.Rd b/man/CellsByIdentities.Rd
index 4223624a..3fdbad3d 100644
--- a/man/CellsByIdentities.Rd
+++ b/man/CellsByIdentities.Rd
@@ -4,7 +4,7 @@
 \alias{CellsByIdentities}
 \title{Get cell names grouped by identity class}
 \usage{
-CellsByIdentities(object, idents = NULL, cells = NULL)
+CellsByIdentities(object, idents = NULL, cells = NULL, return.null = FALSE)
 }
 \arguments{
 \item{object}{A Seurat object}
@@ -13,6 +13,9 @@ CellsByIdentities(object, idents = NULL, cells = NULL)
 defaults to all identity class levels}
 
 \item{cells}{A vector of cells to grouping to}
+
+\item{return.null}{If no cells are request, return a \code{NULL};
+by default, throws an error}
 }
 \value{
 A named list where names are identity classes and values are vectors
diff --git a/man/CreateAssayObject.Rd b/man/CreateAssayObject.Rd
index d792d1e5..ae6aadde 100644
--- a/man/CreateAssayObject.Rd
+++ b/man/CreateAssayObject.Rd
@@ -4,7 +4,7 @@
 \alias{CreateAssayObject}
 \title{Create an Assay object}
 \usage{
-CreateAssayObject(counts, data, min.cells = 0, min.features = 0)
+CreateAssayObject(counts, data, min.cells = 0, min.features = 0, ...)
 }
 \arguments{
 \item{counts}{Unnormalized data such as raw counts or TPMs}
@@ -17,6 +17,8 @@ new object with a lower cutoff.}
 
 \item{min.features}{Include cells where at least this many features are
 detected.}
+
+\item{...}{Arguments passed to \code{\link{as.sparse}}}
 }
 \value{
 A \code{\link{Assay}} object
diff --git a/man/CreateSeuratObject.Rd b/man/CreateSeuratObject.Rd
index e65df510..179a3393 100644
--- a/man/CreateSeuratObject.Rd
+++ b/man/CreateSeuratObject.Rd
@@ -25,6 +25,7 @@ CreateSeuratObject(
   meta.data = NULL,
   min.cells = 0,
   min.features = 0,
+  row.names = NULL,
   ...
 )
 
@@ -70,6 +71,10 @@ new object with a lower cutoff.}
 
 \item{min.features}{Include cells where at least this many features are
 detected.}
+
+\item{row.names}{When \code{counts} is a \code{data.frame} or
+\code{data.frame}-derived object: an optional vector of feature names to be
+used}
 }
 \value{
 A \code{\link{Seurat}} object
diff --git a/man/IsMatrixEmpty.Rd b/man/IsMatrixEmpty.Rd
new file mode 100644
index 00000000..ddfa4bf2
--- /dev/null
+++ b/man/IsMatrixEmpty.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{IsMatrixEmpty}
+\alias{IsMatrixEmpty}
+\title{Check if a matrix is empty}
+\usage{
+IsMatrixEmpty(x)
+}
+\arguments{
+\item{x}{A matrix}
+}
+\value{
+Whether or not \code{x} is empty
+}
+\description{
+Takes a matrix and asks if it's empty (either 0x0 or 1x1 with a value of NA)
+}
+\examples{
+IsMatrixEmpty(new("matrix"))
+IsMatrixEmpty(matrix())
+IsMatrixEmpty(matrix(1:3))
+
+}
+\concept{utils}
diff --git a/man/Seurat-methods.Rd b/man/Seurat-methods.Rd
index abee528d..fa352364 100644
--- a/man/Seurat-methods.Rd
+++ b/man/Seurat-methods.Rd
@@ -9,6 +9,7 @@
 \alias{[[.Seurat}
 \alias{dim.Seurat}
 \alias{dimnames.Seurat}
+\alias{head.Seurat}
 \alias{merge.Seurat}
 \alias{merge}
 \alias{MergeSeurat}
@@ -16,6 +17,7 @@
 \alias{names.Seurat}
 \alias{subset.Seurat}
 \alias{subset}
+\alias{tail.Seurat}
 \alias{[[<-,Seurat-method}
 \alias{colMeans,Seurat-method}
 \alias{colSums,Seurat-method}
@@ -38,6 +40,8 @@
 
 \method{dimnames}{Seurat}(x)
 
+\method{head}{Seurat}(x, n = 10L, ...)
+
 \method{merge}{Seurat}(
   x = NULL,
   y = NULL,
@@ -50,7 +54,17 @@
 
 \method{names}{Seurat}(x)
 
-\method{subset}{Seurat}(x, subset, cells = NULL, features = NULL, idents = NULL, ...)
+\method{subset}{Seurat}(
+  x,
+  subset,
+  cells = NULL,
+  features = NULL,
+  idents = NULL,
+  return.null = FALSE,
+  ...
+)
+
+\method{tail}{Seurat}(x, n = 10L, ...)
 
 \S4method{[[}{Seurat}(x, i, j, ...) <- value
 
@@ -92,6 +106,8 @@ can pass \code{NULL} to remove metadata or an associated object}
 
 \item{drop}{See \code{\link[base]{drop}}}
 
+\item{n}{The number of rows of metadata to return}
+
 \item{y}{A single \code{Seurat} object or a list of \code{Seurat} objects}
 
 \item{add.cell.ids}{A character vector of \code{length(x = c(x, y))};
@@ -111,6 +127,9 @@ are shared across all objects.}
 
 \item{idents}{A vector of identity classes to keep}
 
+\item{return.null}{If no cells are request, return a \code{NULL};
+by default, throws an error}
+
 \item{na.rm}{logical.  Should missing values (including \code{NaN})
     be omitted from the calculations?}
 
@@ -142,6 +161,8 @@ the same across all assays
 \strong{note}: while the features change depending on the active assay, the
 cell names remain the same across all assays
 
+\code{head}: The first \code{n} rows of cell-level metadata
+
 \code{merge}: Merged object
 
 \code{names}: The names of all \code{\link{Assay}},
@@ -150,6 +171,8 @@ objects in the \code{Seurat} object
 
 \code{subset}: A subsetted \code{Seurat} object
 
+\code{tail}: The last \code{n} rows of cell-level metadata
+
 \code{[[<-}: \code{x} with the metadata or associated objects added
 as \code{i}; if \code{value} is \code{NULL}, removes metadata or associated
 object \code{i} from object \code{x}
@@ -178,12 +201,16 @@ packages
 
 \item \code{dimnames.Seurat}: The cell and feature names for the active assay
 
+\item \code{head.Seurat}: Get the first rows of cell-level metadata
+
 \item \code{merge.Seurat}: Merge two or more \code{Seurat} objects together
 
 \item \code{names.Seurat}: Common associated objects
 
 \item \code{subset.Seurat}: Subset a \code{\link{Seurat}} object
 
+\item \code{tail.Seurat}: Get the last rows of cell-level metadata
+
 \item \code{[[<-,Seurat-method}: Add cell-level metadata or associated objects
 
 \item \code{colMeans,Seurat-method}: Calculate \code{\link[base]{colMeans}} on a
@@ -254,6 +281,9 @@ rownames(pbmc_small)
 # Get the cell names of an object
 colnames(pbmc_small)
 
+# Get the first 10 rows of cell-level metadata
+head(pbmc_small)
+
 # `merge' examples
 # merge two objects
 merge(pbmc_small, y = pbmc_small)
@@ -269,6 +299,9 @@ subset(pbmc_small, idents = '0', invert = TRUE)
 subset(pbmc_small, subset = MS4A1 > 3, slot = 'counts')
 subset(pbmc_small, features = VariableFeatures(object = pbmc_small))
 
+# Get the last 10 rows of cell-level metadata
+tail(pbmc_small)
+
 head(colMeans(pbmc_small))
 
 head(colSums(pbmc_small))
diff --git a/man/SeuratObject-package.Rd b/man/SeuratObject-package.Rd
index 75f8f7e8..7048921d 100644
--- a/man/SeuratObject-package.Rd
+++ b/man/SeuratObject-package.Rd
@@ -9,10 +9,10 @@
 Defines S4 classes for single-cell genomic data and associated
   information, such as dimensionality reduction embeddings, nearest-neighbor
   graphs, and spatially-resolved coordinates. Provides data access methods and
-  R-native hooks to ensure the Seurat object is familiar to other R users. See 
-  Satija R, Farrell J, Gennert D, et al (2015) <doi:10.1038/nbt.3192>, 
-  Macosko E, Basu A, Satija R, et al (2015) <doi:10.1016/j.cell.2015.05.002>, 
-  and Stuart T, Butler A, et al (2019) <doi:10.1016/j.cell.2019.05.031> for 
+  R-native hooks to ensure the Seurat object is familiar to other R users. See
+  Satija R, Farrell J, Gennert D, et al (2015) <doi:10.1038/nbt.3192>,
+  Macosko E, Basu A, Satija R, et al (2015) <doi:10.1016/j.cell.2015.05.002>,
+  and Stuart T, Butler A, et al (2019) <doi:10.1016/j.cell.2019.05.031> for
   more details.
 }
 \seealso{
diff --git a/man/Version.Rd b/man/Version.Rd
new file mode 100644
index 00000000..df106335
--- /dev/null
+++ b/man/Version.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/generics.R, R/seurat.R
+\name{Version}
+\alias{Version}
+\alias{Version.Seurat}
+\title{Get Version Information}
+\usage{
+Version(object, ...)
+
+\method{Version}{Seurat}(object, ...)
+}
+\arguments{
+\item{object}{An object}
+
+\item{...}{Arguments passed to other methods}
+}
+\description{
+Get Version Information
+}
+\examples{
+Version(pbmc_small)
+
+}
+\concept{data-access}
diff --git a/man/WhichCells.Rd b/man/WhichCells.Rd
index 7f131be9..f4d22552 100644
--- a/man/WhichCells.Rd
+++ b/man/WhichCells.Rd
@@ -25,7 +25,12 @@ WhichCells(object, ...)
 \arguments{
 \item{object}{An object}
 
-\item{...}{Arguments passed to other methods}
+\item{...}{
+  Arguments passed on to \code{\link[=CellsByIdentities]{CellsByIdentities}}
+  \describe{
+    \item{\code{return.null}}{If no cells are request, return a \code{NULL};
+by default, throws an error}
+  }}
 
 \item{cells}{Subset of cell names}
 
@@ -40,9 +45,9 @@ between numbers are present in the feature name}
 
 \item{slot}{Slot to pull feature data for}
 
-\item{downsample}{Maximum number of cells per identity class, default is \code{Inf};
-downsampling will happen after all other operations, including inverting the
-cell selection}
+\item{downsample}{Maximum number of cells per identity class, default is
+\code{Inf}; downsampling will happen after all other operations, including
+inverting the cell selection}
 
 \item{seed}{Random seed for downsampling. If NULL, does not set a seed}
 }
diff --git a/man/as.sparse.Rd b/man/as.sparse.Rd
index 4977f3a9..ab0fadc1 100644
--- a/man/as.sparse.Rd
+++ b/man/as.sparse.Rd
@@ -9,7 +9,7 @@
 \usage{
 as.sparse(x, ...)
 
-\method{as.sparse}{data.frame}(x, ...)
+\method{as.sparse}{data.frame}(x, row.names = NULL, ...)
 
 \method{as.sparse}{Matrix}(x, ...)
 
@@ -19,6 +19,9 @@ as.sparse(x, ...)
 \item{x}{An object}
 
 \item{...}{Arguments passed to other methods}
+
+\item{row.names}{\code{NULL} or a character vector giving the row names for
+the data; missing values are not allowed}
 }
 \value{
 A sparse representation of the input data
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 00000000..72a15ce4
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,4 @@
+library(testthat)
+library(Seurat)
+
+test_check("Seurat")
diff --git a/tests/testthat/test_objects.R b/tests/testthat/test_objects.R
new file mode 100644
index 00000000..150f3dbe
--- /dev/null
+++ b/tests/testthat/test_objects.R
@@ -0,0 +1,343 @@
+# Tests for functions in objects.R
+
+# Tests for interacting with the meta.data slot
+# ------------------------------------------------------------------------------
+context("Metadata")
+
+data("pbmc_small")
+
+pbmc_small <- suppressWarnings(suppressMessages(UpdateSeuratObject(pbmc_small)))
+cluster_letters <- LETTERS[Idents(object = pbmc_small)]
+names(cluster_letters) <- colnames(x = pbmc_small)
+cluster_letters_shuffled <- sample(x = cluster_letters)
+
+test_that("AddMetaData adds in cell-level vector properly ", {
+  pbmc_small <- AddMetaData(object = pbmc_small, metadata = cluster_letters, col.name = 'letter.idents')
+  expect_equal(pbmc_small$letter.idents, cluster_letters)
+  pbmc_small <- AddMetaData(object = pbmc_small, metadata = cluster_letters_shuffled, col.name = 'letter.idents.shuffled')
+  expect_equal(pbmc_small$letter.idents, pbmc_small$letter.idents.shuffled)
+})
+
+cluster_letters_df <- data.frame(A = cluster_letters, B = cluster_letters_shuffled)
+test_that("AddMetaData adds in data frame properly for cell-level metadata", {
+  pbmc_small <- AddMetaData(object = pbmc_small, metadata = cluster_letters_df)
+  expect_equal(pbmc_small[[c("A", "B")]], cluster_letters_df)
+})
+
+feature_letters <- sample(x = LETTERS, size = nrow(x = pbmc_small[["RNA"]]), replace = TRUE)
+names(feature_letters) <- rownames(x = pbmc_small[["RNA"]])
+feature_letters_shuffled <- sample(x = feature_letters)
+
+test_that("AddMetaData adds feature level metadata", {
+  pbmc_small[["RNA"]] <- AddMetaData(object = pbmc_small[["RNA"]], metadata = feature_letters, col.name = 'feature_letters')
+  expect_equal(pbmc_small[["RNA"]][["feature_letters", drop = TRUE]], feature_letters)
+  pbmc_small[["RNA"]] <- AddMetaData(object = pbmc_small[["RNA"]], metadata = feature_letters_shuffled, col.name = 'feature_letters_shuffled')
+  expect_equal(pbmc_small[["RNA"]][["feature_letters", drop = TRUE]], pbmc_small[["RNA"]][["feature_letters_shuffled", drop = TRUE]])
+})
+
+feature_letters_df <- data.frame(A = feature_letters, B = feature_letters_shuffled)
+test_that("AddMetaData adds in data frame properly for Assays", {
+  pbmc_small[["RNA"]] <- AddMetaData(object = pbmc_small[["RNA"]], metadata = feature_letters_df)
+  expect_equal(pbmc_small[["RNA"]][[c("A", "B")]], feature_letters_df)
+})
+
+test_that("AddMetaData errors", {
+  expect_error(AddMetaData(object = pbmc_small, metadata = cluster_letters, col.name = "RNA"))
+  expect_error(AddMetaData(object = pbmc_small, metadata = c(unname(cluster_letters), "A"), col.name = "letter.idents"))
+  expect_error(AddMetaData(object = pbmc_small, metadata = feature_letters, col.name = "letter.idents"))
+  expect_error(AddMetaData(object = pbmc_small[["RNA"]], metadata = cluster_letters, col.name = "letter.idents"))
+})
+
+# Tests for creating an Assay object
+# ------------------------------------------------------------------------------
+context("CreateAssayObject")
+
+pbmc.raw <- GetAssayData(object = pbmc_small[["RNA"]], slot = "counts")
+rna.assay <- CreateAssayObject(counts = pbmc.raw)
+rna.assay2 <- CreateAssayObject(data = pbmc.raw)
+
+test_that("CreateAssayObject works as expected", {
+  expect_equal(dim(x = rna.assay), c(230, 80))
+  expect_equal(rownames(x = rna.assay), rownames(x = pbmc.raw))
+  expect_equal(colnames(x = rna.assay), colnames(x = pbmc.raw))
+  expect_equal(GetAssayData(object = rna.assay, slot = "counts"), pbmc.raw)
+  expect_equal(GetAssayData(object = rna.assay, slot = "data"), pbmc.raw)
+  expect_equal(GetAssayData(object = rna.assay, slot = "scale.data"), new(Class = "matrix"))
+  expect_equal(dim(rna.assay[[]]), c(230, 0))
+  expect_equal(rownames(x = rna.assay[[]]), rownames(x = rna.assay))
+  expect_equal(VariableFeatures(object = rna.assay), vector())
+  expect_equal(rna.assay@misc, list())
+  expect_equal(GetAssayData(object = rna.assay2, slot = "counts"), new(Class = "matrix"))
+})
+
+rna.assay2 <- CreateAssayObject(counts = pbmc.raw, min.cells = 10, min.features = 30)
+test_that("CreateAssayObject filtering works", {
+  expect_equal(dim(x = rna.assay2), c(163, 77))
+  expect_true(all(rowSums(GetAssayData(object = rna.assay2, slot = "counts")) >= 10))
+  expect_true(all(colSums(GetAssayData(object = rna.assay2, slot = "counts")) >= 30))
+})
+
+test_that("CreateAssayObject catches improper input", {
+  expect_error(CreateAssayObject())
+  expect_error(CreateAssayObject(counts = pbmc.raw, data = pbmc.raw))
+  pbmc.raw2 <- cbind(pbmc.raw[, 1:10], pbmc.raw[, 1:10])
+  expect_warning(CreateAssayObject(counts = pbmc.raw2))
+  expect_warning(CreateAssayObject(data = pbmc.raw2))
+  pbmc.raw2 <- rbind(pbmc.raw[1:10, ], pbmc.raw[1:10, ])
+  expect_warning(CreateAssayObject(counts = pbmc.raw2))
+  expect_warning(CreateAssayObject(data = pbmc.raw2))
+  pbmc.raw2 <- pbmc.raw
+  colnames(x = pbmc.raw2) <- c()
+  expect_error(CreateAssayObject(counts = pbmc.raw2))
+  expect_error(CreateAssayObject(data = pbmc.raw2))
+  pbmc.raw2 <- pbmc.raw
+  rownames(x = pbmc.raw2) <- c()
+  expect_error(CreateAssayObject(counts = pbmc.raw2))
+  expect_error(CreateAssayObject(data = pbmc.raw2))
+  pbmc.raw.mat <- as.matrix(x = pbmc.raw)
+  pbmc.raw.df <- as.data.frame(x = pbmc.raw.mat)
+  rna.assay3 <- CreateAssayObject(counts = pbmc.raw.df)
+  rna.assay4 <- CreateAssayObject(counts = pbmc.raw.mat)
+  expect_is(object = GetAssayData(object = rna.assay3, slot = "counts"), class = "dgCMatrix")
+  expect_is(object = GetAssayData(object = rna.assay4, slot = "counts"), class = "dgCMatrix")
+  pbmc.raw.underscores <- pbmc.raw
+  rownames(pbmc.raw.underscores) <- gsub(pattern = "-", replacement = "_", x = rownames(pbmc.raw.underscores))
+  expect_warning(CreateAssayObject(counts = pbmc.raw.underscores))
+})
+
+# Tests for creating an DimReduc object
+# ------------------------------------------------------------------------------
+context("CreateDimReducObject")
+
+pca <- pbmc_small[["pca"]]
+Key(object = pca) <- 'PC_'
+
+test_that("CreateDimReducObject works", {
+  pca.dr <- CreateDimReducObject(
+    embeddings = Embeddings(object = pca),
+    loadings = Loadings(object = pca),
+    projected = Loadings(object = pca, projected = TRUE),
+    assay = "RNA"
+  )
+  expect_equal(Embeddings(object = pca.dr), Embeddings(object = pca))
+  expect_equal(Loadings(object = pca.dr), Loadings(object = pca))
+  expect_equal(Loadings(object = pca.dr, projected = TRUE), Loadings(object = pca, projected = TRUE))
+  expect_equal(Key(object = pca.dr), "PC_")
+  expect_equal(pca.dr@assay.used, "RNA")
+})
+
+test_that("CreateDimReducObject catches improper input", {
+  bad.embeddings <- Embeddings(object = pca)
+  colnames(x = bad.embeddings) <- paste0("PCA", 1:ncol(x = bad.embeddings))
+  expect_warning(CreateDimReducObject(embeddings = bad.embeddings, key = "PC"))
+  colnames(x = bad.embeddings) <- paste0("PC", 1:ncol(x = bad.embeddings), "X")
+  suppressWarnings(expect_error(CreateDimReducObject(embeddings = bad.embeddings, key = "PC")))
+  suppressWarnings(expect_error(CreateDimReducObject(embeddings = bad.embeddings)))
+})
+
+# Tests for creating a Seurat object
+# ------------------------------------------------------------------------------
+context("CreateSeuratObject")
+
+colnames(x = pbmc.raw) <- paste0(colnames(x = pbmc.raw), "-", pbmc_small$groups)
+metadata.test <- pbmc_small[[]][, 5:7]
+rownames(x = metadata.test) <- colnames(x = pbmc.raw)
+
+test_that("CreateSeuratObject works", {
+  seurat.object <- CreateSeuratObject(
+    counts = pbmc.raw,
+    project = "TESTING",
+    assay = "RNA.TEST",
+    names.field = 2,
+    names.delim = "-",
+    meta.data = metadata.test
+  )
+  expect_equal(seurat.object[[]][, 4:6], metadata.test)
+  expect_equal(seurat.object@project.name, "TESTING")
+  expect_equal(names(x = seurat.object), "RNA.TEST")
+  expect_equal(as.vector(x = unname(obj = Idents(object = seurat.object))), unname(pbmc_small$groups))
+})
+
+test_that("CreateSeuratObject handles bad names.field/names.delim", {
+  expect_warning(seurat.object <- CreateSeuratObject(
+    counts = pbmc.raw[1:5,1:5],
+    names.field = 3,
+    names.delim = ":",
+    meta.data = metadata.test
+  ))
+})
+
+# Tests for creating a Seurat object
+# ------------------------------------------------------------------------------
+context("Merging")
+
+pbmc.assay <- pbmc_small[["RNA"]]
+x <- merge(x = pbmc.assay, y = pbmc.assay)
+
+test_that("Merging Assays works properly", {
+  expect_equal(dim(GetAssayData(object = x, slot = "counts")), c(230, 160))
+  expect_equal(dim(GetAssayData(object = x, slot = "data")), c(230, 160))
+  expect_equal(GetAssayData(object = x, slot = "scale.data"), new(Class = "matrix"))
+  expect_equal(Key(object = x), "rna_")
+  expect_equal(VariableFeatures(object = x), vector())
+  expect_equal(x[[]], data.frame(row.names = rownames(x = pbmc.assay)))
+})
+
+pbmc.assay2 <- pbmc.assay
+pbmc.assay2@counts <- new("dgCMatrix")
+test_that("Merging Assays handles case when counts not present", {
+  y <- merge(x = pbmc.assay2, y = pbmc.assay)
+  expect_equal(unname(colSums(x = GetAssayData(object = y, slot = "counts"))[1:80]), rep.int(x = 0, times = 80))
+  z <- merge(x = pbmc.assay2, pbmc.assay2)
+  expect_equal(Matrix::nnzero(x = GetAssayData(object = z, slot = "counts")), 0)
+})
+
+pbmc.assay2 <- pbmc.assay
+pbmc.assay2@data <- new("dgCMatrix")
+test_that("Merging Assays handles case when data not present", {
+  y <- merge(x = pbmc.assay2, y = pbmc.assay, merge.data = TRUE)
+  expect_equal(unname(colSums(x = GetAssayData(object = y, slot = "data"))[1:80]), rep.int(x = 0, times = 80))
+  z <- merge(x = pbmc.assay2, y = pbmc.assay2, merge.data = TRUE)
+  expect_equal(Matrix::nnzero(x = GetAssayData(object = z, slot = "data")), 0)
+})
+
+# # Tests for Neighbor object
+# # ------------------------------------------------------------------------------
+# context("Neighbor")
+#
+# # converting to Graph and back
+#
+# n.rann.ob <- NNHelper(
+#   data = Embeddings(object = pbmc_small[["pca"]]),
+#   query = Embeddings(object = pbmc_small[["pca"]]),
+#   k = 10,
+#   method = "rann")
+#
+# test_that("Neighbor object methods work", {
+#   expect_equal(dim(x = Indices(object = n.rann.ob)), c(80, 10))
+#   expect_equal(dim(x = n.rann.ob), c(80, 10))
+#   expect_equal(as.numeric(Indices(object = n.rann.ob)[1, 7]), 45, )
+#   expect_equal(dim(x = Distances(object = n.rann.ob)), c(80, 10))
+#   expect_equal(as.numeric(Distances(object = n.rann.ob)[2, 2]), 2.643759, tolerance = 1e-6)
+#   expect_equal(length(x = Cells(x = n.rann.ob)), 80)
+#   expect_equal(Cells(x = n.rann.ob)[c(1, 20, 80)], c("ATGCCAGAACGACT", "TACATCACGCTAAC", "CTTGATTGATCTTC"))
+#   pbmc_small[["n.ob"]] <- n.rann.ob
+#   pbmc_small <- RenameCells(object = pbmc_small, add.cell.id = "test")
+#   expect_equal(Cells(x = pbmc_small[['n.ob']])[1], c("test_ATGCCAGAACGACT"))
+#   expect_equal(TopNeighbors(object = n.rann.ob, cell = "ATGCCAGAACGACT", n = 5)[5], "GATATAACACGCAT")
+#   expect_equal(length(TopNeighbors(object = n.rann.ob, cell = "ATGCCAGAACGACT", n = 7)), 7)
+#   nrg <- as.Graph(x = n.rann.ob)
+#   expect_true(inherits(x = nrg, what = "Graph"))
+#   expect_equal(as.numeric(Distances(object = n.rann.ob)[2, 3]), nrg[2, Indices(object = n.rann.ob)[2, 3]])
+#   nro2 <- as.Neighbor(x = nrg)
+#   expect_true(inherits(x = nro2, what = "Neighbor"))
+#   expect_equal(Distances(object = n.rann.ob)[2, 3], Distances(object = nro2)[2, 3])
+#   expect_equal(Indices(object = n.rann.ob)[1, 6], Indices(object = nro2)[1, 6])
+# })
+#
+# n.annoy.ob <- NNHelper(
+#   data = Embeddings(object = pbmc_small[["pca"]]),
+#   query = Embeddings(object = pbmc_small[["pca"]]),
+#   k = 10,
+#   method = "annoy",
+#   cache.index = TRUE)
+# idx.file <-  tempfile()
+# SaveAnnoyIndex(object = n.annoy.ob, file = idx.file)
+# nao2 <- LoadAnnoyIndex(object = n.annoy.ob, file = idx.file)
+#
+# test_that("Saving/Loading annoy index", {
+#   expect_error(SaveAnnoyIndex(object = n.rann.ob, file = idx.file))
+#   expect_equal(head(Indices(n.annoy.ob)), head(Indices(nao2)))
+#   expect_equal(head(Distances(n.annoy.ob)), head(Distances(nao2)))
+#   expect_false(is.null(x = Index(nao2)))
+# })
+
+# Tests for FetchData
+# ------------------------------------------------------------------------------
+context("FetchData")
+
+# Features to test:
+# able to pull cell embeddings, data, metadata
+# subset of cells
+
+test_that("Fetching a subset of cells works", {
+  x <- FetchData(object = pbmc_small, cells = colnames(x = pbmc_small)[1:10], vars = rownames(x = pbmc_small)[1])
+  expect_equal(rownames(x = x), colnames(x = pbmc_small)[1:10])
+  random.cells <- sample(x = colnames(x = pbmc_small), size = 10)
+  x <- FetchData(object = pbmc_small, cells = random.cells, vars = rownames(x = pbmc_small)[1])
+  expect_equal(rownames(x = x), random.cells)
+  x <- FetchData(object = pbmc_small, cells = 1:10, vars = rownames(x = pbmc_small)[1])
+  expect_equal(rownames(x = x), colnames(x = pbmc_small)[1:10])
+})
+
+suppressWarnings(pbmc_small[["RNA2"]] <- pbmc_small[["RNA"]])
+Key(pbmc_small[["RNA2"]]) <- "rna2_"
+
+test_that("Fetching keyed variables works", {
+  x <- FetchData(object = pbmc_small, vars = c(paste0("rna_", rownames(x = pbmc_small)[1:5]), paste0("rna2_", rownames(x = pbmc_small)[1:5])))
+  expect_equal(colnames(x = x), c(paste0("rna_", rownames(x = pbmc_small)[1:5]), paste0("rna2_", rownames(x = pbmc_small)[1:5])))
+  x <- FetchData(object = pbmc_small, vars = c(paste0("rna_", rownames(x = pbmc_small)[1:5]), paste0("PC_", 1:5)))
+  expect_equal(colnames(x = x), c(paste0("rna_", rownames(x = pbmc_small)[1:5]), paste0("PC_", 1:5)))
+})
+
+test_that("Fetching embeddings/loadings not present returns warning or errors", {
+  expect_warning(FetchData(object = pbmc_small, vars = c("PC_1", "PC_100")))
+  expect_error(FetchData(object = pbmc_small, vars = "PC_100"))
+})
+
+bad.gene <- GetAssayData(object = pbmc_small[["RNA"]], slot = "data")
+rownames(x = bad.gene)[1] <- paste0("rna_", rownames(x = bad.gene)[1])
+pbmc_small[["RNA"]]@data <- bad.gene
+
+# Tests for WhichCells
+# ------------------------------------------------------------------------------
+
+test_that("Specifying cells works", {
+  test.cells <- Cells(x = pbmc_small)[1:10]
+  expect_equal(WhichCells(object = pbmc_small, cells = test.cells), test.cells)
+  expect_equal(WhichCells(object = pbmc_small, cells = test.cells, invert = TRUE), setdiff(Cells(x = pbmc_small), test.cells))
+})
+
+test_that("Specifying idents works", {
+  c12 <- WhichCells(object = pbmc_small, idents = c(1, 2))
+  expect_equal(length(x = c12), 44)
+  expect_equal(c12[44], "CTTGATTGATCTTC")
+  expect_equal(c12, WhichCells(object = pbmc_small, idents = 0, invert = TRUE))
+})
+
+test_that("downsample works", {
+  expect_equal(length(x = WhichCells(object = pbmc_small, downsample = 5)), 15)
+  expect_equal(length(x = WhichCells(object = pbmc_small, downsample = 100)), 80)
+})
+
+test_that("passing an expression works", {
+  lyz.pos <- WhichCells(object = pbmc_small, expression = LYZ > 1)
+  expect_true(all(GetAssayData(object = pbmc_small, slot = "data")["LYZ", lyz.pos] > 1))
+  # multiple values in expression
+  lyz.pos <- WhichCells(object = pbmc_small, expression = LYZ > 1 & groups == "g1")
+  expect_equal(length(x = lyz.pos), 30)
+  expect_equal(lyz.pos[30], "CTTGATTGATCTTC")
+})
+
+# Tests for small other functions
+# ------------------------------------------------------------------------------
+test_that("Top works", {
+  dat <- Embeddings(object = pbmc_small[['pca']])[, 1, drop = FALSE]
+  expect_warning(Top(data = dat, num = 1000, balanced = FALSE))
+  tpc1 <- Top(data = dat, num = 20, balanced = FALSE)
+  expect_equal(length(x = tpc1), 20)
+  expect_equal(tpc1[1], "ACGTGATGCCATGA")
+  expect_equal(tpc1[20], "GTCATACTTCGCCT")
+  tpc1b <- Top(data = dat, num = 20, balanced = TRUE)
+  expect_equal(length(x = tpc1b), 2)
+  expect_equal(names(tpc1b), c("positive", "negative"))
+  expect_equal(length(tpc1b[[1]]), 10)
+  expect_equal(length(tpc1b[[2]]), 10)
+  expect_equal(tpc1b[[1]][1], "GTCATACTTCGCCT")
+  expect_equal(tpc1b[[1]][10], "CTTGATTGATCTTC")
+  expect_equal(tpc1b[[2]][1], "ACGTGATGCCATGA")
+  expect_equal(tpc1b[[2]][10], "ATTGTAGATTCCCG")
+  tpc1.sub <- Top(data = dat[1:79, , drop = FALSE], num = 79, balanced = TRUE)
+  expect_equal(length(tpc1.sub[[1]]), 40)
+  expect_equal(length(tpc1.sub[[2]]), 39)
+})