diff --git a/DESCRIPTION b/DESCRIPTION
index 956e2ffab..2c56f1980 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: Seurat
-Version: 4.0.0
-Date: 2021-01-27
+Version: 4.0.1
+Date: 2021-03-17
 Title: Tools for Single Cell Genomics
 Description: A toolkit for quality control, analysis, and exploration of single cell RNA sequencing data. 'Seurat' aims to enable users to identify and interpret sources of heterogeneity from single cell transcriptomic measurements, and to integrate diverse types of single cell data. 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>, Stuart T, Butler A, et al (2019) <doi:10.1016/j.cell.2019.05.031>, and Hao, Hao, et al (2020) <doi:10.1101/2020.10.12.335331> for more details.
 Authors@R: c(
@@ -68,7 +68,8 @@ Imports:
     sctransform (>= 0.3.2),
     SeuratObject,
     shiny,
-    spatstat,
+    spatstat.core,
+    spatstat.geom,
     stats,
     tibble,
     tools,
diff --git a/NAMESPACE b/NAMESPACE
index 21e7c6faa..4e32d9ac6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -40,6 +40,7 @@ S3method(GetImage,VisiumV1)
 S3method(GetTissueCoordinates,STARmap)
 S3method(GetTissueCoordinates,SlideSeq)
 S3method(GetTissueCoordinates,VisiumV1)
+S3method(HVFInfo,SCTAssay)
 S3method(IntegrateEmbeddings,IntegrationAnchorSet)
 S3method(IntegrateEmbeddings,TransferAnchorSet)
 S3method(MappingScore,AnchorSet)
@@ -261,6 +262,7 @@ export(Radius)
 export(Read10X)
 export(Read10X_Image)
 export(Read10X_h5)
+export(ReadMtx)
 export(ReadSlideSeq)
 export(Reductions)
 export(RegroupIdents)
@@ -555,7 +557,9 @@ importFrom(grid,unit)
 importFrom(grid,viewport)
 importFrom(httr,GET)
 importFrom(httr,accept_json)
+importFrom(httr,build_url)
 importFrom(httr,content)
+importFrom(httr,parse_url)
 importFrom(httr,status_code)
 importFrom(httr,timeout)
 importFrom(ica,icafast)
@@ -636,8 +640,8 @@ importFrom(shiny,sliderInput)
 importFrom(shiny,stopApp)
 importFrom(shiny,updateSelectInput)
 importFrom(shiny,verbatimTextOutput)
-importFrom(spatstat,markvario)
-importFrom(spatstat,ppp)
+importFrom(spatstat.core,markvario)
+importFrom(spatstat.geom,ppp)
 importFrom(stats,aggregate)
 importFrom(stats,anova)
 importFrom(stats,approxfun)
diff --git a/NEWS.md b/NEWS.md
index 28229078b..382a2eb64 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,23 @@
-# Seurat 4.0.0 (2020-01-27) 
+# Seurat 4.0.1 (2020-03-17)
 ## Added
-- Expose `FoldChange()` component in `FindMarkers()`. 
+- Add direction option to `PlotClusterTree()`
+- Add `cols` parameter to `JackStrawPlot()`
+- Add `ReadMtx()` to read local and remote mtx files with associated cell and feature name files
+
+## Changes
+- Equality added to differential expression thresholds in `FindMarkers` (e.g, >= logfc.threshold rather than >)
+- `Read10X()` now prepends dataset number for first dataset when reading multiple datasets
+- Bug fix for `subset.AnchorSet()`
+- Bug fix for fold change values in `FindMarkers()` when setting a different pseudocount ([#4111](https://github.com/satijalab/seurat/pull/4111))
+- Bug fix for `RunLDA()` related to proper passing of assay parameter.
+- When using `order=TRUE` in `SingleDimPlot()`, print NA points under all others.
+- Remove default parameter value for `data.dir` in `Read10X()`
+- Import spatstat fxns from subpackages (spatstat.core, spatstat.geom)
+- `RunUMAP` now checks for graph/neighbor consistency
+
+# Seurat 4.0.0 (2020-01-27)
+## Added
+- Expose `FoldChange()` component in `FindMarkers()`.
 - Add the `merge.DimReduc` method
 - Add `IntegrateEmbeddings()` to correct embeddings of `DimReduc`s
 - Add `ProjectUMAP()` to project query cells into a reference UMAP space
@@ -11,20 +28,21 @@
 - Add `FindSubCluster()` to further cluster existing clusters
 - Add supervised PCA functionality via `RunSPCA()`
 - Add functionality to enable weighted nearest neighbor analyses via `FindMultiModalNeighbors()`
-- Add neighbor visualization plot via `NNPlot()`. 
+- Add neighbor visualization plot via `NNPlot()`.
 - Add `PredictAssay()` to impute expression or embeddings from nearest neighbors
 - Add `Graphs()` function to access the names of the stored Graph objects or pull a specific one
 - Add checks for NA, NaN, logical, non-integer, and infinite values during CreateAssayObject and NormalizeData.default
 - Add `AnnotateAnchors()` to aid in AnchorSet interpretation as well as `subset.AnchorSet()`
 - Add flexibility of choice for cell column in `Read10X()`
 - Add rasterization option to `FeatureScatter()` and `VariableFeaturePlot()`
+- Add step1 feature parameters in the SCTModel via `PrepVSTResults()`
 
 ## Changes
 - Default neighbor finding algorithm changed from "rann" to "annoy"
 - Default `ncells` parameter in `SCTransform()` changed to 5000
 - Default fold change in `FindMarkers()` changed from ln to log2
 - Implementation improvements to `AverageExpression()`
-- `AnchorSet` class re-implemented as a virtual class from which `IntegrationAnchorSet` and `TransferAnchorSet` now inherit. 
+- `AnchorSet` class re-implemented as a virtual class from which `IntegrationAnchorSet` and `TransferAnchorSet` now inherit.
 - Point size in `VlnPlot()` now set automatically if not specified
 - Return the sample.tree properly when integrating with a single reference dataset
 - Replace `as.character.quosure` usage with `as_label` due to deprecation
@@ -32,6 +50,7 @@
 - Default rasterization limit in `DimPlot()` and `FeaturePlot()` changed from 50,000 to 100,000
 - `SCTransform()` now returns a formalized `Assay` subclass `SCTAssay()`
 - When using `normalization.method='SCT'` in `FindTransferAnchors()`, normalize query using reference SCT model when possible.
+- Change default Neighbor name in `FindNeighbors` to `Assay.nn`
 
 ## Removed
 - `CreateGeneActivityMatrix` replaced by `Signac::GeneActivity()`
diff --git a/R/clustering.R b/R/clustering.R
index 690b1dfea..de0cc96c5 100644
--- a/R/clustering.R
+++ b/R/clustering.R
@@ -559,7 +559,11 @@ FindNeighbors.default <- function(
   # find the k-nearest neighbors for each single cell
   if (!distance.matrix) {
     if (verbose) {
-      message("Computing nearest neighbor graph")
+      if (return.neighbor) {
+        message("Computing nearest neighbors")
+      } else {
+        message("Computing nearest neighbor graph")
+      }
     }
     nn.ranked <- NNHelper(
       data = object,
@@ -709,6 +713,11 @@ FindNeighbors.dist <- function(
 #' @param do.plot Plot SNN graph on tSNE coordinates
 #' @param graph.name Optional naming parameter for stored (S)NN graph
 #' (or Neighbor object, if return.neighbor = TRUE). Default is assay.name_(s)nn.
+#' To store both the neighbor graph and the shared nearest neighbor (SNN) graph,
+#' you must supply a vector containing two names to the \code{graph.name}
+#' parameter. The first element in the vector will be used to store the nearest
+#' neighbor (NN) graph, and the second element used to store the SNN graph. If
+#' only one name is supplied, only the NN graph is stored.
 #'
 #' @importFrom igraph graph.adjacency plot.igraph E
 #'
@@ -787,7 +796,15 @@ FindNeighbors.Seurat <- function(
   if (length(x = neighbor.graphs) == 1) {
     neighbor.graphs <- list(nn = neighbor.graphs)
   }
-  graph.name <- graph.name %||% paste0(assay, "_", names(x = neighbor.graphs))
+  graph.name <- graph.name %||%
+    if (return.neighbor) {
+      paste0(assay, ".", names(x = neighbor.graphs))
+    } else {
+      paste0(assay, "_", names(x = neighbor.graphs))
+    }
+  if (length(x = graph.name) == 1) {
+    message("Only one graph name supplied, storing nearest-neighbor graph only")
+  }
   for (ii in 1:length(x = graph.name)) {
     if (inherits(x = neighbor.graphs[[ii]], what = "Graph")) {
       DefaultAssay(object = neighbor.graphs[[ii]]) <- assay
diff --git a/R/convenience.R b/R/convenience.R
index dac371fc5..69e8a7ddf 100644
--- a/R/convenience.R
+++ b/R/convenience.R
@@ -32,6 +32,7 @@ PCAPlot <- function(object, ...) {
 
 #' @rdname SpatialPlot
 #' @concept convenience
+#' @concept spatial
 #' @export
 #'
 SpatialDimPlot <- function(
@@ -82,6 +83,7 @@ SpatialDimPlot <- function(
 
 #' @rdname SpatialPlot
 #' @concept convenience
+#' @concept spatial
 #' @export
 #'
 SpatialFeaturePlot <- function(
diff --git a/R/differential_expression.R b/R/differential_expression.R
index 3abf832e9..25a91e562 100644
--- a/R/differential_expression.R
+++ b/R/differential_expression.R
@@ -513,14 +513,14 @@ FindMarkers.default <- function(
   # feature selection (based on percentages)
   alpha.min <- pmax(fc.results$pct.1, fc.results$pct.2)
   names(x = alpha.min) <- rownames(x = fc.results)
-  features <- names(x = which(x = alpha.min > min.pct))
+  features <- names(x = which(x = alpha.min >= min.pct))
   if (length(x = features) == 0) {
     warning("No features pass min.pct threshold; returning empty data.frame")
     return(fc.results[features, ])
   }
   alpha.diff <- alpha.min - pmin(fc.results$pct.1, fc.results$pct.2)
   features <- names(
-    x = which(x = alpha.min > min.pct & alpha.diff > min.diff.pct)
+    x = which(x = alpha.min >= min.pct & alpha.diff >= min.diff.pct)
   )
   if (length(x = features) == 0) {
     warning("No features pass min.diff.pct threshold; returning empty data.frame")
@@ -531,9 +531,9 @@ FindMarkers.default <- function(
     total.diff <- fc.results[, 1] #first column is logFC
     names(total.diff) <- rownames(fc.results)
     features.diff <- if (only.pos) {
-      names(x = which(x = total.diff > logfc.threshold))
+      names(x = which(x = total.diff >= logfc.threshold))
     } else {
-      names(x = which(x = abs(x = total.diff) > logfc.threshold))
+      names(x = which(x = abs(x = total.diff) >= logfc.threshold))
     }
     features <- intersect(x = features, y = features.diff)
     if (length(x = features) == 0) {
@@ -627,6 +627,7 @@ FindMarkers.Assay <- function(
     cells.1 = cells.1,
     cells.2 = cells.2,
     features = features,
+    pseudocount.use = pseudocount.use,
     mean.fxn = mean.fxn,
     fc.name = fc.name,
     base = base
@@ -680,6 +681,7 @@ FindMarkers.DimReduc <- function(
   min.cells.group = 3,
   pseudocount.use = 1,
   mean.fxn = rowMeans,
+  fc.name = NULL,
   ...
 
 ) {
@@ -705,7 +707,9 @@ FindMarkers.DimReduc <- function(
     object = object,
     cells.1 = cells.1,
     cells.2 = cells.2,
-    features = features
+    features = features,
+    mean.fxn = mean.fxn,
+    fc.name = fc.name
   )
   # subsample cell groups if they are too large
   if (max.cells.per.ident < Inf) {
diff --git a/R/dimensional_reduction.R b/R/dimensional_reduction.R
index 1ce3b6dd5..a91ce1923 100644
--- a/R/dimensional_reduction.R
+++ b/R/dimensional_reduction.R
@@ -1654,9 +1654,27 @@ RunUMAP.Seurat <- function(
         call. = FALSE
       )
     }
-  }  else if (!is.null( x = nn.name)){
+  }  else if (!is.null(x = nn.name)) {
+    if (!inherits(x = object[[nn.name]], what = "Neighbor")) {
+      stop(
+        "Please specify a Neighbor object name, ",
+        "instead of the name of a ", 
+        class(object[[nn.name]]), 
+        " object",   
+        call. = FALSE
+      )
+    }
     data.use <- object[[nn.name]]
   } else if (!is.null(x = graph)) {
+    if (!inherits(x = object[[graph]], what = "Graph")) {
+      stop(
+        "Please specify a Graph object name, ",
+        "instead of the name of a ", 
+        class(object[[graph]]), 
+        " object",   
+        call. = FALSE
+      )
+    }
     data.use <- object[[graph]]
   } else {
     stop("Please specify one of dims, features, or graph")
diff --git a/R/integration.R b/R/integration.R
index 13afef5b1..19195e145 100644
--- a/R/integration.R
+++ b/R/integration.R
@@ -492,7 +492,8 @@ FindIntegrationAnchors <- function(
 #' }
 #' @param reference.reduction Name of dimensional reduction to use from the
 #' reference if running the pcaproject workflow. Optionally enables reuse of
-#' precomputed reference dimensional reduction.
+#' precomputed reference dimensional reduction. If NULL (default), use a PCA
+#' computed on the reference object.
 #' @param project.query Project the PCA from the query dataset onto the
 #' reference. Use only in rare cases where the query dataset has a much larger
 #' cell number, but the reference dataset has a unique assay for transfer. In
@@ -627,25 +628,25 @@ FindTransferAnchors <- function(
   if (normalization.method == "SCT") {
       # ensure all residuals required are computed
       query <- suppressWarnings(expr = GetResidual(object = query, assay = query.assay, features = features, verbose = FALSE))
-        if (is.null(x = reference.reduction)) {
-      reference <- suppressWarnings(expr = GetResidual(object = reference, assay = reference.assay, features = features, verbose = FALSE))
-      features <- intersect(
-        x = features,
-        y = intersect(
-          x = rownames(x = GetAssayData(object = query[[query.assay]], slot = "scale.data")),
-          y = rownames(x = GetAssayData(object = reference[[reference.assay]], slot = "scale.data"))
+      if (is.null(x = reference.reduction)) {
+        reference <- suppressWarnings(expr = GetResidual(object = reference, assay = reference.assay, features = features, verbose = FALSE))
+        features <- intersect(
+          x = features,
+          y = intersect(
+            x = rownames(x = GetAssayData(object = query[[query.assay]], slot = "scale.data")),
+            y = rownames(x = GetAssayData(object = reference[[reference.assay]], slot = "scale.data"))
+          )
         )
-      )
-      reference[[reference.assay]] <- as(
-        object = CreateAssayObject(
-          data = GetAssayData(object = reference[[reference.assay]], slot = "scale.data")[features, ]),
-        Class = "SCTAssay"
-      )
-      reference <- SetAssayData(
-        object = reference,
-        slot = "scale.data",
-        assay = reference.assay,
-        new.data =  as.matrix(x = GetAssayData(object = reference[[reference.assay]], slot = "data"))
+        reference[[reference.assay]] <- as(
+          object = CreateAssayObject(
+            data = GetAssayData(object = reference[[reference.assay]], slot = "scale.data")[features, ]),
+          Class = "SCTAssay"
+        )
+        reference <- SetAssayData(
+          object = reference,
+          slot = "scale.data",
+          assay = reference.assay,
+          new.data =  as.matrix(x = GetAssayData(object = reference[[reference.assay]], slot = "data"))
       )
     }
     query[[query.assay]] <- as(
@@ -2675,6 +2676,7 @@ TransferData <- function(
 #' @rdname AnnotateAnchors
 #' @export
 #' @method AnnotateAnchors default
+#' @concept integration
 #'
 AnnotateAnchors.default <- function(
   anchors,
diff --git a/R/mixscape.R b/R/mixscape.R
index 289e99a81..f0a68eabc 100644
--- a/R/mixscape.R
+++ b/R/mixscape.R
@@ -532,7 +532,7 @@ RunLDA.default <- function(
 #'
 RunLDA.Assay <- function(
   object,
-  assay,
+  assay = NULL,
   labels,
   features = NULL,
   verbose = TRUE,
@@ -595,6 +595,7 @@ RunLDA.Seurat <- function(
   assay.data <- GetAssay(object = object, assay = assay)
   reduction.data <- RunLDA(
     object = assay.data,
+    assay = assay,
     labels = labels,
     features = features,
     verbose = verbose,
diff --git a/R/objects.R b/R/objects.R
index 6ba81b047..f79dcb24f 100644
--- a/R/objects.R
+++ b/R/objects.R
@@ -164,6 +164,7 @@ IntegrationData <- setClass(
 #'
 #' @name SCTAssay-class
 #' @rdname SCTAssay-class
+#' @concept objects
 #'
 #' @examples
 #' \dontrun{
@@ -194,6 +195,7 @@ SCTModel <- setClass(
 #'
 #' @name SCTAssay-class
 #' @rdname SCTAssay-class
+#' @concept objects
 #'
 #' @examples
 #' # SCTAssay objects are generated from SCTransform
@@ -218,6 +220,7 @@ SCTAssay <- setClass(
 #'
 #' @rdname ScaleFactors
 #' @concept objects
+#' @concept spatial
 #' @export
 #'
 scalefactors <- function(spot, fiducial, hires, lowres) {
@@ -239,7 +242,7 @@ setOldClass(Classes = c('scalefactors'))
 #'
 #' @inheritSection SeuratObject::SpatialImage Slots
 #' @slot coordinates ...
-#' @slot ...
+#' @concept spatial
 #'
 SlideSeq <- setClass(
   Class = 'SlideSeq',
@@ -253,7 +256,8 @@ SlideSeq <- setClass(
 #'
 #'
 #' @inheritSection SeuratObject::SpatialImage Slots
-#' @slot ...
+#' @concept objects
+#' @concept spatial
 #'
 STARmap <- setClass(
   Class = 'STARmap',
@@ -279,6 +283,7 @@ STARmap <- setClass(
 #' @name VisiumV1-class
 #' @rdname VisiumV1-class
 #' @concept objects
+#' @concept spatial
 #' @exportClass VisiumV1
 #'
 VisiumV1 <- setClass(
@@ -306,7 +311,6 @@ setClass(Class = 'SliceImage', contains = 'VisiumV1')
 #' named by image name.
 #'
 #' @return A vector of cell names
-#' @concept objects
 #'
 #' @examples
 #' \dontrun{
@@ -345,7 +349,7 @@ CellsByImage <- function(object, images = NULL, unlist = FALSE) {
 #' @importFrom Matrix colSums rowSums
 #'
 #' @export
-#'
+#' @concept objects
 #'
 CreateSCTAssayObject <- function(
   counts,
@@ -533,6 +537,7 @@ DietSeurat <- function(
 #' circular filter
 #'
 #' @concept objects
+#' @concept spatial
 #' @examples
 #' \dontrun{
 #' # This example uses the ssHippo dataset which you can download
@@ -1189,6 +1194,7 @@ Cells.SCTModel <- function(x) {
 
 #' @rdname Cells
 #' @concept objects
+#' @concept spatial
 #' @method Cells SlideSeq
 #' @export
 #'
@@ -1200,6 +1206,7 @@ Cells.SlideSeq <- function(x) {
 
 #' @rdname Cells
 #' @concept objects
+#' @concept spatial
 #' @method Cells STARmap
 #' @export
 #'
@@ -1250,6 +1257,7 @@ GetAssay.Seurat <- function(object, assay = NULL, ...) {
 #' @rdname GetImage
 #' @method GetImage SlideSeq
 #' @concept objects
+#' @concept spatial
 #' @export
 #'
 #' @seealso \code{\link[SeuratObject:GetImage]{SeuratObject::GetImage}}
@@ -1266,6 +1274,7 @@ GetImage.SlideSeq <- function(
 #' @rdname GetImage
 #' @method GetImage STARmap
 #' @concept objects
+#' @concept spatial
 #' @export
 #'
 GetImage.STARmap <- function(
@@ -1283,6 +1292,7 @@ GetImage.STARmap <- function(
 #'
 #' @rdname GetImage
 #' @concept objects
+#' @concept spatial
 #' @method GetImage VisiumV1
 #' @export
 #'
@@ -1326,6 +1336,7 @@ GetImage.VisiumV1 <- function(
 #' @rdname GetTissueCoordinates
 #' @method GetTissueCoordinates SlideSeq
 #' @concept objects
+#' @concept spatial
 #' @export
 #'
 #' @seealso \code{\link[SeuratObject:GetTissueCoordinates]{SeuratObject::GetTissueCoordinates}}
@@ -1344,6 +1355,7 @@ GetTissueCoordinates.SlideSeq <- function(object, ...) {
 #' @rdname GetTissueCoordinates
 #' @method GetTissueCoordinates STARmap
 #' @concept objects
+#' @concept spatial
 #' @export
 #'
 GetTissueCoordinates.STARmap <- function(object, qhulls = FALSE, ...) {
@@ -1360,6 +1372,7 @@ GetTissueCoordinates.STARmap <- function(object, qhulls = FALSE, ...) {
 #' @rdname GetTissueCoordinates
 #' @method GetTissueCoordinates VisiumV1
 #' @concept objects
+#' @concept spatial
 #' @export
 #'
 GetTissueCoordinates.VisiumV1 <- function(
@@ -1380,12 +1393,49 @@ GetTissueCoordinates.VisiumV1 <- function(
   return(coordinates)
 }
 
+#' Get Variable Feature Information
+#'
+#' Get variable feature information from \code{\link{SCTAssay}} objects
+#'
+#' @inheritParams SeuratObject::HVFInfo
+#'
+#' @export
+#' @method HVFInfo SCTAssay
+#'
+#' @seealso \code{\link[SeuratObject]{HVFInfo}}
+#'
+#' @examples
+#' # Get the HVF info directly from an SCTAssay object
+#' pbmc_small <- SCTransform(pbmc_small)
+#' HVFInfo(pbmc_small[["SCT"]], selection.method = 'sct')[1:5, ]
+#'
+HVFInfo.SCTAssay <- function(object, selection.method, status = FALSE, ...) {
+  CheckDots(...)
+  disp.methods <- c('mean.var.plot', 'dispersion', 'disp')
+  if (tolower(x = selection.method) %in% disp.methods) {
+    selection.method <- 'mvp'
+  }
+  selection.method <- switch(
+    EXPR = tolower(x = selection.method),
+    'sctransform' = 'sct',
+    selection.method
+  )
+  vars <- c('gmean', 'variance', 'residual_variance')
+  hvf.info <- SCTResults(object = object, slot = "feature.attributes")[,vars]
+  if (status) {
+    hvf.info$variable <- FALSE
+    hvf.info[VariableFeatures(object = object), "variable"] <-  TRUE
+  }
+  return(hvf.info)
+}
+
 #' Get Spot Radius
 #'
 #' @inheritParams SeuratObject::Radius
 #'
 #' @rdname Radius
 #' @concept objects
+#' @concept spatial
 #' @method Radius SlideSeq
 #' @export
 #'
@@ -1397,6 +1447,7 @@ Radius.SlideSeq <- function(object) {
 
 #' @rdname Radius
 #' @concept objects
+#' @concept spatial
 #' @method Radius STARmap
 #' @export
 #'
@@ -1406,6 +1457,7 @@ Radius.STARmap <- function(object) {
 
 #' @rdname Radius
 #' @concept objects
+#' @concept spatial
 #' @method Radius VisiumV1
 #' @export
 #'
@@ -1415,6 +1467,7 @@ Radius.VisiumV1 <- function(object) {
 
 #' @rdname RenameCells
 #' @export
+#' @concept objects
 #' @method RenameCells SCTAssay
 #'
 RenameCells.SCTAssay <- function(object, new.names = NULL, ...) {
@@ -1448,6 +1501,7 @@ RenameCells.SCTAssay <- function(object, new.names = NULL, ...) {
 #' @inheritParams SeuratObject::RenameCells
 #'
 #' @rdname RenameCells
+#' @concept objects
 #' @method RenameCells SlideSeq
 #' @export
 #'
@@ -1458,6 +1512,7 @@ RenameCells.SlideSeq <- function(object, new.names = NULL, ...) {
 }
 
 #' @rdname RenameCells
+#' @concept objects
 #' @method RenameCells STARmap
 #' @export
 #'
@@ -1471,6 +1526,7 @@ RenameCells.STARmap <- function(object, new.names = NULL, ...) {
 }
 
 #' @rdname RenameCells
+#' @concept objects
 #' @method RenameCells VisiumV1
 #' @export
 #'
@@ -1505,6 +1561,7 @@ SCTResults.SCTModel <- function(object, slot, ...) {
 }
 
 #' @rdname SCTResults
+#' @concept objects
 #' @export
 #' @method SCTResults<- SCTModel
 #'
@@ -1531,6 +1588,7 @@ SCTResults.SCTModel <- function(object, slot, ...) {
 #' the slot directly).
 #'
 #' @rdname SCTResults
+#' @concept objects
 #' @export
 #' @method SCTResults SCTAssay
 #'
@@ -1554,6 +1612,7 @@ SCTResults.SCTAssay <- function(object, slot, model = NULL, ...) {
 }
 
 #' @rdname SCTResults
+#' @concept objects
 #' @export
 #' @method SCTResults<- SCTAssay
 #'
@@ -1588,6 +1647,7 @@ SCTResults.SCTAssay <- function(object, slot, model = NULL, ...) {
 #'
 #' @rdname SCTResults
 #' @export
+#' @concept objects
 #' @method SCTResults Seurat
 #'
 SCTResults.Seurat <- function(object, assay = "SCT", slot, model = NULL, ...) {
@@ -1598,6 +1658,7 @@ SCTResults.Seurat <- function(object, assay = "SCT", slot, model = NULL, ...) {
 #' @rdname ScaleFactors
 #' @method ScaleFactors VisiumV1
 #' @export
+#' @concept spatial
 #'
 ScaleFactors.VisiumV1 <- function(object, ...) {
   return(slot(object = object, name = 'scale.factors'))
@@ -1606,6 +1667,7 @@ ScaleFactors.VisiumV1 <- function(object, ...) {
 #' @rdname ScaleFactors
 #' @method ScaleFactors VisiumV1
 #' @export
+#' @concept spatial
 #'
 ScaleFactors.VisiumV1 <- function(object, ...) {
   return(slot(object = object, name = 'scale.factors'))
@@ -1677,6 +1739,7 @@ dim.VisiumV1 <- function(x) {
 #' @return \code{levels}: SCT model names
 #'
 #' @export
+#' @concept objects
 #' @method levels SCTAssay
 #'
 #' @examples
@@ -1699,6 +1762,7 @@ levels.SCTAssay <- function(x) {
 #' @return \code{levels<-}: \code{x} with updated SCT model names
 #'
 #' @export
+#' @concept objects
 #' @method levels<- SCTAssay
 #'
 "levels<-.SCTAssay" <- function(x, value) {
@@ -1726,6 +1790,7 @@ levels.SCTAssay <- function(x) {
 #' populated with NAs.
 #' @export
 #' @method merge SCTAssay
+#' @concept objects
 #'
 merge.SCTAssay <- function(
   x = NULL,
@@ -1861,6 +1926,7 @@ merge.SCTAssay <- function(
 #'
 #' @export
 #' @method subset AnchorSet
+#' @concept objects
 #'
 subset.AnchorSet <- function(
   x,
@@ -1922,11 +1988,14 @@ subset.AnchorSet <- function(
   }
   # Filter based on ident pairings
   if (!is.null(x = group.by)) {
-    anchors <- AnnotateAnchors(object = x, vars = group.by)
+    anchors <- AnnotateAnchors(anchors = x, vars = group.by)
     if (!is.null(x = disallowed.ident.pairs) && !is.null(x = ident.matrix)) {
       stop("Please use either disallowed.ident.pairs OR ident.matrix, not both.")
     }
-    unique.ids <- unique(x = c(anchors[, paste0("cell1.", group.by)], anchors[, paste0("cell2.", group.by)]))
+    unique.ids <- unique(x = c(
+      as.character(x = anchors[, paste0("cell1.", group.by)]),
+      as.character(x = anchors[, paste0("cell2.", group.by)]))
+    )
     unique.ids <- unique.ids[!is.na(x = unique.ids)]
     num.ids <- length(x = unique.ids)
     if (!is.null(x = disallowed.ident.pairs)) {
@@ -1963,6 +2032,7 @@ subset.AnchorSet <- function(
 
 #' @export
 #' @method subset SCTAssay
+#' @concept objects
 #'
 subset.SCTAssay <- function(x, cells = NULL, features = NULL, ...) {
   x <- NextMethod()
@@ -2008,13 +2078,12 @@ subset.VisiumV1 <- function(x, cells, ...) {
   return(x)
 }
 
-
-
 #' Update pre-V4 Assays generated with SCTransform in the Seurat to the new
 #' SCTAssay class
 #
 #' @param object A Seurat object
 #' @export
+#' @concept objects
 #' @return A Seurat object with updated SCTAssays
 #'
 UpdateSCTAssays <- function(object) {
@@ -2373,6 +2442,11 @@ PrepVSTResults <- function(vst.res, cell.names) {
       y = rownames(feature.attrs)
     )
     feature.attrs[genes_step1,"genes_log_gmean_step1"] <- TRUE
+
+    # add parameters from step1
+    feature.attrs[, paste0("step1_", colnames(vst.res$model_pars))] <- NA
+    feature.attrs[genes_step1, paste0("step1_", colnames(vst.res$model_pars))] <- vst.res$model_pars[genes_step1,]
+
   }
   # Prepare clipping information
   clips <- list(
diff --git a/R/preprocessing.R b/R/preprocessing.R
index 301cb3d3d..74574c850 100644
--- a/R/preprocessing.R
+++ b/R/preprocessing.R
@@ -771,7 +771,7 @@ MULTIseqDemux <- function(
 #' }
 #'
 Read10X <- function(
-  data.dir = NULL,
+  data.dir,
   gene.column = 2,
   cell.column = 1,
   unique.features = TRUE,
@@ -821,7 +821,7 @@ Read10X <- function(
       )))
     }
     if (is.null(x = names(x = data.dir))) {
-      if (i < 2) {
+      if (length(x = data.dir) < 2) {
         colnames(x = data) <- cell.names
       } else {
         colnames(x = data) <- paste0(i, "_", cell.names)
@@ -1028,6 +1028,201 @@ Read10X_Image <- function(image.dir, filter.matrix = TRUE, ...) {
   ))
 }
 
+#' Load in data from remote or local mtx files
+#'
+#' Enables easy loading of sparse data matrices
+#'
+#' @param mtx Name or remote URL of the mtx file
+#' @param cells Name or remote URL of the cells/barcodes file
+#' @param features Name or remote URL of the features/genes file
+#' @param cell.column Specify which column of cells file to use for cell names; default is 1
+#' @param feature.column Specify which column of features files to use for feature/gene names; default is 2
+#' @param skip.cell Number of lines to skip in the cells file before beginning to read cell names
+#' @param skip.feature Number of lines to skip in the features file before beginning to gene names
+#' @param unique.features Make feature names unique (default TRUE)
+#' @param strip.suffix Remove trailing "-1" if present in all cell barcodes.
+#'
+#' @return A sparse matrix containing the expression data.
+#'
+#' @importFrom Matrix readMM
+#' @importFrom utils read.delim
+#' @importFrom httr build_url parse_url
+#' @importFrom tools file_ext
+#'
+#'
+#' @export
+#' @concept preprocessing
+#'
+#' @examples
+#' \dontrun{
+#' # For local files:
+#'
+#' expression_matrix <- ReadMtx(
+#'   mtx = "count_matrix.mtx.gz", features = "features.tsv.gz",
+#'   cells = "barcodes.tsv.gz"
+#' )
+#' seurat_object <- CreateSeuratObject(counts = expression_matrix)
+#'
+#' # For remote files:
+#'
+#' expression_matrix <- ReadMtx(mtx = "http://localhost/matrix.mtx",
+#' cells = "http://localhost/barcodes.tsv",
+#' features = "http://localhost/genes.tsv")
+#' seurat_object <- CreateSeuratObject(counts = data)
+#' }
+#'
+ReadMtx <- function(
+  mtx,
+  cells,
+  features,
+  cell.column = 1,
+  feature.column = 2,
+  skip.cell = 0,
+  skip.feature = 0,
+  unique.features = TRUE,
+  strip.suffix = FALSE
+) {
+  all.files <- list(
+    "expression matrix" = mtx,
+    "barcode list" = cells,
+    "feature list" = features
+  )
+  for (i in seq_along(along.with = all.files)) {
+    uri <- all.files[[i]]
+    err <- paste("Cannot find", names(x = all.files)[i], "at", uri)
+    uri <- build_url(url = parse_url(url = uri))
+    if (grepl(pattern = '^:///', x = uri)) {
+      uri <- gsub(pattern = '^:///', replacement = '', x = uri)
+      if (!file.exists(uri)) {
+        stop(err, call. = FALSE)
+      }
+    } else {
+      if (!Online(url = uri, seconds = 2L)) {
+        stop(err, call. = FALSE)
+      }
+      if (file_ext(uri) == 'gz') {
+        con <- url(description = uri)
+        uri <- gzcon(con = con, text = TRUE)
+      }
+    }
+    all.files[[i]] <- uri
+  }
+  cell.barcodes <- read.table(
+    file = all.files[['barcode list']],
+    header = FALSE,
+    sep = '\t',
+    row.names = NULL,
+    skip = skip.cell
+  )
+  feature.names <- read.table(
+    file = all.files[['feature list']],
+    header = FALSE,
+    sep = '\t',
+    row.names = NULL,
+    skip = skip.feature
+  )
+  # read barcodes
+  bcols <- ncol(x = cell.barcodes)
+  if (bcols < cell.column) {
+    stop(
+      "cell.column was set to ",
+      cell.column,
+      " but ",
+      cells,
+      " only has ",
+      bcols,
+      " columns.",
+      " Try setting the cell.column argument to a value <= to ",
+      bcols,
+      "."
+    )
+  }
+  cell.names <- cell.barcodes[, cell.column]
+  if (all(grepl(pattern = "\\-1$", x = cell.names)) & strip.suffix) {
+    cell.names <- as.vector(x = as.character(x = sapply(
+      X = cell.names,
+      FUN = ExtractField,
+      field = 1,
+      delim = "-"
+    )))
+  }
+  # read features
+  fcols <- ncol(x = feature.names)
+  if (fcols < feature.column) {
+    stop(
+      "feature.column was set to ",
+      feature.column,
+      " but ",
+      features,
+      " only has ",
+      fcols, " column(s).",
+      " Try setting the feature.column argument to a value <= to ",
+      fcols,
+      "."
+    )
+  }
+  if (any(is.na(x = feature.names[, feature.column]))) {
+    na.features <- which(x = is.na(x = feature.names[, feature.column]))
+    replacement.column <- ifelse(test = feature.column == 2, yes = 1, no = 2)
+    if (replacement.column > fcols) {
+      stop(
+        "Some features names are NA in column ",
+        feature.column,
+        ". Try specifiying a different column.",
+        call. = FALSE
+        )
+    } else {
+      warning(
+        "Some features names are NA in column ",
+        feature.column,
+        ". Replacing NA names with ID from column ",
+        replacement.column,
+        ".",
+        call. = FALSE
+        )
+    }
+    feature.names[na.features, feature.column] <- feature.names[na.features, replacement.column]
+  }
+  feature.names <- feature.names[, feature.column]
+  if (unique.features) {
+    feature.names <- make.unique(names = feature.names)
+  }
+  data <- readMM(file = all.files[['expression matrix']])
+  if (length(x = cell.names) != ncol(x = data)) {
+    stop(
+      "Matrix has ",
+      ncol(data),
+      " columns but found ", length(cell.names),
+      " barcodes. ",
+      ifelse(
+        test = length(x = cell.names) > ncol(x = data),
+        yes = "Try increasing `skip.cell`. ",
+        no = ""
+      ),
+      call. = FALSE
+      )
+  }
+  if (length(x = feature.names) != nrow(x = data)) {
+    stop(
+      "Matrix has ",
+      ncol(data),
+      " rows but found ", length(feature.names),
+      " features. ",
+      ifelse(
+        test = length(x = feature.names) > nrow(x = data),
+        yes = "Try increasing `skip.feature`. ",
+        no = ""
+      ),
+      call. = FALSE
+      )
+  }
+
+  colnames(x = data) <- cell.names
+  rownames(x = data) <- feature.names
+  data <- as(data, Class = "dgCMatrix")
+  return(data)
+}
+
 #' Load Slide-seq spatial data
 #'
 #' @param coord.file Path to csv file containing bead coordinate positions
@@ -1107,7 +1302,8 @@ RelativeCounts <- function(data, scale.factor = 1, verbose = TRUE) {
 #' @param data Matrix containing the data used as "marks" (e.g. gene expression)
 #' @param ... Arguments passed to markvario
 #'
-#' @importFrom spatstat markvario ppp
+#' @importFrom spatstat.core markvario
+#' @importFrom spatstat.geom ppp
 #'
 #' @export
 #' @concept preprocessing
@@ -1962,6 +2158,7 @@ FindVariableFeatures.Seurat <- function(
 #' @method FindSpatiallyVariableFeatures default
 #' @rdname FindSpatiallyVariableFeatures
 #' @concept preprocessing
+#' @concept spatial
 #' @export
 #'
 #'
@@ -2014,6 +2211,7 @@ FindSpatiallyVariableFeatures.default <- function(
 #' @method FindSpatiallyVariableFeatures Assay
 #' @rdname FindSpatiallyVariableFeatures
 #' @concept preprocessing
+#' @concept spatial
 #' @export
 #'
 FindSpatiallyVariableFeatures.Assay <- function(
@@ -2078,6 +2276,7 @@ FindSpatiallyVariableFeatures.Assay <- function(
 #' @method FindSpatiallyVariableFeatures Seurat
 #' @rdname FindSpatiallyVariableFeatures
 #' @concept preprocessing
+#' @concept spatial
 #' @export
 #'
 FindSpatiallyVariableFeatures.Seurat <- function(
diff --git a/R/utilities.R b/R/utilities.R
index fe227d1f7..287f07019 100644
--- a/R/utilities.R
+++ b/R/utilities.R
@@ -1984,14 +1984,14 @@ ModifyParam <- function(param, value) {
   env2[[param]] <- value
 }
 
-# Give hints for old paramters and their newer counterparts
+# Give hints for old parameters and their newer counterparts
 #
 # This is a non-exhaustive list. If your function isn't working properly based
 # on the parameters you give it, please read the documentation for your function
 #
-# @param param A vector of paramters to get hints for
+# @param param A vector of parameters to get hints for
 #
-# @return Parameter hints for the specified paramters
+# @return Parameter hints for the specified parameters
 #
 OldParamHints <- function(param) {
   param.conversion <- c(
@@ -2017,6 +2017,40 @@ OldParamHints <- function(param) {
   return(param.conversion[param])
 }
 
+# Check if a web resource is available
+#
+# @param url A URL
+# @param strict Perform a strict web availability test
+# @param seconds Timeout in seconds
+#
+# @return \code{TRUE} if \url{is available} otherwise \code{FALSE}
+#
+#' @importFrom httr GET status_code timeout
+#
+# @keywords internal
+#
+Online <- function(url, strict = FALSE, seconds = 5L) {
+  if (isTRUE(x = strict)) {
+    code <- 200L
+    comp <- identical
+  } else {
+    code <- 404L
+    comp <- Negate(f = identical)
+  }
+  request <- tryCatch(
+    expr = GET(url = url, timeout(seconds = seconds)),
+    error = function(err) {
+      code <- if (grepl(pattern = '^Timeout was reached', x = err$message)) {
+        408L
+      } else {
+        404L
+      }
+      return(code)
+    }
+  )
+  return(comp(x = status_code(x = request), y = code))
+}
+
 # Check the existence of a package
 #
 # @param ... Package names
@@ -2044,7 +2078,7 @@ PackageCheck <- function(..., error = TRUE) {
 
 # Parenting parameters from one environment to the next
 #
-# This function allows one to modifiy a parameter in a parent environement
+# This function allows one to modify a parameter in a parent environment
 # The primary use of this is to ensure logging functions store correct parameters
 # if they've been modified by a child function or method
 #
diff --git a/R/visualization.R b/R/visualization.R
index fb4cd9fae..25d74c015 100644
--- a/R/visualization.R
+++ b/R/visualization.R
@@ -1210,7 +1210,7 @@ FeaturePlot <- function(
   }
   # Make the plots
   for (i in 1:length(x = levels(x = data$split))) {
-    # Figre out which split we're working with
+    # Figure out which split we're working with
     ident <- levels(x = data$split)[i]
     data.plot <- data[as.character(x = data$split) == ident, , drop = FALSE]
     # Blend expression values
@@ -2010,10 +2010,25 @@ VariableFeaturePlot <- function(
     pt.size = pt.size,
     raster = raster
   )
+  if (length(x = unique(x = var.status)) == 1) {
+    switch(
+      EXPR = var.status[1],
+      'yes' = {
+        cols <- cols[2]
+        labels.legend <- 'Variable'
+      },
+      'no' = {
+        cols <- cols[1]
+        labels.legend <- 'Non-variable'
+      }
+    )
+  } else {
+    labels.legend <- c('Non-variable', 'Variable')
+  }
   plot <- plot +
     labs(title = NULL, x = axis.labels[1], y = axis.labels[2]) +
     scale_color_manual(
-      labels = paste(c('Non-variable', 'Variable'), 'count:', table(var.status)),
+      labels = paste(labels.legend, 'count:', table(var.status)),
       values = cols
     )
   if (log) {
@@ -2084,6 +2099,7 @@ PolyDimPlot <- function(
 #'
 #' @export
 #' @concept visualization
+#' @concept spatial
 #'
 PolyFeaturePlot <- function(
   object,
@@ -2207,6 +2223,7 @@ PolyFeaturePlot <- function(
 #'
 #' @export
 #' @concept visualization
+#' @concept spatial
 #'
 #' @examples
 #' \dontrun{
@@ -2421,7 +2438,7 @@ LinkedDimPlot <- function(
 #'
 #' @export
 #' @concept visualization
-#'
+#' @concept spatial
 LinkedFeaturePlot <- function(
   object,
   feature,
@@ -2554,6 +2571,7 @@ LinkedFeaturePlot <- function(
 #'
 #' @export
 #' @concept visualization
+#' @concept spatial
 #'
 ISpatialDimPlot <- function(
   object,
@@ -2679,7 +2697,7 @@ ISpatialDimPlot <- function(
 #'
 #' @export
 #' @concept visualization
-#'
+#' @concept spatial
 ISpatialFeaturePlot <- function(
   object,
   feature,
@@ -2926,6 +2944,7 @@ ISpatialFeaturePlot <- function(
 #' @importFrom patchwork wrap_plots
 #' @export
 #' @concept visualization
+#' @concept spatial
 #'
 #' @examples
 #' \dontrun{
@@ -3688,6 +3707,10 @@ GroupCorrelationPlot <- function(
 #'
 #' @param object Seurat object
 #' @param dims Dims to plot
+#' @param cols Vector of colors, each color corresponds to an individual PC. This may also be a single character
+#' or numeric value corresponding to a palette as specified by \code{\link[RColorBrewer]{brewer.pal.info}}.
+#' By default, ggplot2 assigns colors. We also include a number of palettes from the pals package.
+#' See \code{\link{DiscretePalette}} for details.
 #' @param reduction reduction to pull jackstraw info from
 #' @param xmax X-axis maximum on each QQ plot.
 #' @param ymax Y-axis maximum on each QQ plot.
@@ -3698,6 +3721,7 @@ GroupCorrelationPlot <- function(
 #' @seealso \code{\link{ScoreJackStraw}}
 #'
 #' @importFrom stats qunif
+#' @importFrom scales hue_pal
 #' @importFrom ggplot2 ggplot aes_string stat_qq labs xlim ylim
 #' coord_flip geom_abline guides guide_legend
 #' @importFrom cowplot theme_cowplot
@@ -3712,6 +3736,7 @@ GroupCorrelationPlot <- function(
 JackStrawPlot <- function(
   object,
   dims = 1:5,
+  cols = NULL,
   reduction = 'pca',
   xmax = 0.1,
   ymax = 0.3
@@ -3740,9 +3765,16 @@ JackStrawPlot <- function(
     x = data.plot$PC.Score,
     levels = paste0("PC ", score.df[, "PC"], ": ", sprintf("%1.3g", score.df[, "Score"]))
   )
+  if (is.null(x = cols)) {
+    cols <- hue_pal()(length(x = dims))
+  }
+  if (length(x = cols) < length(x = dims)) {
+    stop("Not enough colors for the number of dims selected")
+  }
   gp <- ggplot(data = data.plot, mapping = aes_string(sample = 'Value', color = 'PC.Score')) +
     stat_qq(distribution = qunif) +
     labs(x = "Theoretical [runif(1000)]", y = "Empirical") +
+    scale_color_manual(values = cols) +
     xlim(0, ymax) +
     ylim(0, xmax) +
     coord_flip() +
@@ -3757,6 +3789,8 @@ JackStrawPlot <- function(
 #' Plots previously computed tree (from BuildClusterTree)
 #'
 #' @param object Seurat object
+#' @param direction A character string specifying the direction of the tree (default is downwards)
+#' Possible options: "rightwards", "leftwards", "upwards", and "downwards".
 #' @param \dots Additional arguments to
 #' \code{\link[ape:plot.phylo]{ape::plot.phylo}}
 #'
@@ -3770,7 +3804,7 @@ JackStrawPlot <- function(
 #' pbmc_small <- BuildClusterTree(object = pbmc_small)
 #' PlotClusterTree(object = pbmc_small)
 #'
-PlotClusterTree <- function(object, ...) {
+PlotClusterTree <- function(object, direction = "downwards", ...) {
   if (!PackageCheck('ape', error = FALSE)) {
     stop(cluster.ape, call. = FALSE)
   }
@@ -3778,7 +3812,7 @@ PlotClusterTree <- function(object, ...) {
     stop("Phylogenetic tree does not exist, build using BuildClusterTree")
   }
   data.tree <- Tool(object = object, slot = "BuildClusterTree")
-  ape::plot.phylo(x = data.tree, direction = "downwards", ...)
+  ape::plot.phylo(x = data.tree, direction = direction, ...)
   ape::nodelabels()
 }
 
@@ -7042,7 +7076,7 @@ SingleDimPlot <- function(
   if (!is.null(x = order) && !is.null(x = col.by)) {
     if (typeof(x = order) == "logical") {
       if (order) {
-        data <- data[order(data[, col.by]), ]
+        data <- data[order(!is.na(x = data[, col.by]), data[, col.by]), ]
       }
     } else {
       order <- rev(x = c(
diff --git a/README.md b/README.md
index 619d7af2c..ae591d755 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-[![Build Status](https://travis-ci.com/satijalab/seurat.svg?branch=master)](https://travis-ci.com/satijalab/seurat)
+[![Build Status](https://travis-ci.com/satijalab/seurat.svg?branch=master)](https://travis-ci.com/github/satijalab/seurat)
 [![AppVeyor build status](https://ci.appveyor.com/api/projects/status/github/satijalab/seurat?branch=master&svg=true)](https://ci.appveyor.com/project/satijalab/seurat)
 [![CRAN Version](https://www.r-pkg.org/badges/version/Seurat)](https://cran.r-project.org/package=Seurat)
 [![CRAN Downloads](https://cranlogs.r-pkg.org/badges/Seurat)](https://cran.r-project.org/package=Seurat)
diff --git a/_pkgdown.yaml b/_pkgdown.yaml
index 7f753f206..4f2139554 100644
--- a/_pkgdown.yaml
+++ b/_pkgdown.yaml
@@ -113,6 +113,12 @@ reference:
 - contents:
   - has_concept("integration")
 
+- title: "Spatial"
+  desc: >
+    Functions related to the analysis of spatially-resolved single-cell data
+- contents:
+  - has_concept("spatial")
+
 - title: "Mixscape"
   desc: >
       Functions related to the mixscape algorithm
diff --git a/cran-comments.md b/cran-comments.md
index e7b0941f7..eed9624bd 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,7 +1,6 @@
-# Seurat v4.0.0
+# Seurat v4.0.1
 
 ## Test environments
-* local Ubuntu 16.04.6 install, R 4.0.3
 * local Ubuntu 18.04.4 install, R 4.0.3
 * local Ubuntu 20.04 install, R 4.0.3
 * local Windows 10 install, R 4.0.3
@@ -15,8 +14,8 @@ There were no ERRORs, WARNINGs, or NOTEs
 
 ## Downstream dependencies
 
-There is one package that depends on Seurat: tidyseurat; this update does not impact its functionality
+There are no packages that depend on Seurat.
 
-There are three packages that imports Seurat: scMappR, Signac, and SoupX; this update does not impact their functionality
+There are six packages that imports Seurat: CDSeq, DUBStepR, scMappR, Signac, SignacX, and SoupX; this update does not impact their functionality
 
-There are ten packages that suggest Seurat: BisqueRNA, clustree, DIscBIO, nanny, rliger, Rmagic, scSorter, singleCellHaystack, treefit, VAM; this update does not impact their functionality.
+There are eleven packages that suggest Seurat: BisqueRNA, clustree, conos, DIscBIO, dyngen, nanny, rliger, Rmagic, singleCellHaystack, treefit, VAM; this update does not impact their functionality.
diff --git a/man/AnnotateAnchors.Rd b/man/AnnotateAnchors.Rd
index 6fb16c299..6a6a5a375 100644
--- a/man/AnnotateAnchors.Rd
+++ b/man/AnnotateAnchors.Rd
@@ -61,3 +61,4 @@ metadata
 \description{
 Add info to anchor matrix
 }
+\concept{integration}
diff --git a/man/Cells.Rd b/man/Cells.Rd
index aa8598f44..8da74d202 100644
--- a/man/Cells.Rd
+++ b/man/Cells.Rd
@@ -25,3 +25,4 @@ Get Cell Names
 \code{\link[SeuratObject:Cells]{SeuratObject::Cells}}
 }
 \concept{objects}
+\concept{spatial}
diff --git a/man/CellsByImage.Rd b/man/CellsByImage.Rd
index 3ad9f224a..ca0acaa3c 100644
--- a/man/CellsByImage.Rd
+++ b/man/CellsByImage.Rd
@@ -26,4 +26,3 @@ CellsByImage(object = object, images = "slice1")
 }
 
 }
-\concept{objects}
diff --git a/man/CreateSCTAssayObject.Rd b/man/CreateSCTAssayObject.Rd
index 2d9eeeb55..c3c1406c8 100644
--- a/man/CreateSCTAssayObject.Rd
+++ b/man/CreateSCTAssayObject.Rd
@@ -40,3 +40,4 @@ The expected format of the input matrix is features x cells.
 Non-unique cell or feature names are not allowed. Please make unique before
 calling this function.
 }
+\concept{objects}
diff --git a/man/FilterSlideSeq.Rd b/man/FilterSlideSeq.Rd
index 7187e8ede..a9b873db1 100644
--- a/man/FilterSlideSeq.Rd
+++ b/man/FilterSlideSeq.Rd
@@ -50,3 +50,4 @@ ssHippo.filtered <- FilterSlideSeq(ssHippo, radius = 2300)
 }
 }
 \concept{objects}
+\concept{spatial}
diff --git a/man/FindMarkers.Rd b/man/FindMarkers.Rd
index 19b23a781..099691bb4 100644
--- a/man/FindMarkers.Rd
+++ b/man/FindMarkers.Rd
@@ -76,6 +76,7 @@ FindMarkers(object, ...)
   min.cells.group = 3,
   pseudocount.use = 1,
   mean.fxn = rowMeans,
+  fc.name = NULL,
   ...
 )
 
diff --git a/man/FindNeighbors.Rd b/man/FindNeighbors.Rd
index e4c080b63..4771bf893 100644
--- a/man/FindNeighbors.Rd
+++ b/man/FindNeighbors.Rd
@@ -150,7 +150,12 @@ is \code{NULL}}
 \item{do.plot}{Plot SNN graph on tSNE coordinates}
 
 \item{graph.name}{Optional naming parameter for stored (S)NN graph
-(or Neighbor object, if return.neighbor = TRUE). Default is assay.name_(s)nn.}
+(or Neighbor object, if return.neighbor = TRUE). Default is assay.name_(s)nn.
+To store both the neighbor graph and the shared nearest neighbor (SNN) graph,
+you must supply a vector containing two names to the \code{graph.name}
+parameter. The first element in the vector will be used to store the nearest
+neighbor (NN) graph, and the second element used to store the SNN graph. If
+only one name is supplied, only the NN graph is stored.}
 }
 \value{
 This function can either return a \code{\link{Neighbor}} object
diff --git a/man/FindSpatiallyVariableFeatures.Rd b/man/FindSpatiallyVariableFeatures.Rd
index 7a710c248..a94625155 100644
--- a/man/FindSpatiallyVariableFeatures.Rd
+++ b/man/FindSpatiallyVariableFeatures.Rd
@@ -89,3 +89,4 @@ Identify features whose variability in expression can be explained to some
 degree by spatial location.
 }
 \concept{preprocessing}
+\concept{spatial}
diff --git a/man/FindTransferAnchors.Rd b/man/FindTransferAnchors.Rd
index 8d42ea0f8..43e7ae9f2 100644
--- a/man/FindTransferAnchors.Rd
+++ b/man/FindTransferAnchors.Rd
@@ -59,7 +59,8 @@ Options are:
 
 \item{reference.reduction}{Name of dimensional reduction to use from the
 reference if running the pcaproject workflow. Optionally enables reuse of
-precomputed reference dimensional reduction.}
+precomputed reference dimensional reduction. If NULL (default), use a PCA
+computed on the reference object.}
 
 \item{project.query}{Project the PCA from the query dataset onto the
 reference. Use only in rare cases where the query dataset has a much larger
diff --git a/man/GetImage.Rd b/man/GetImage.Rd
index c2daba055..a0d134863 100644
--- a/man/GetImage.Rd
+++ b/man/GetImage.Rd
@@ -27,3 +27,4 @@ Get Image Data
 \code{\link[SeuratObject:GetImage]{SeuratObject::GetImage}}
 }
 \concept{objects}
+\concept{spatial}
diff --git a/man/GetTissueCoordinates.Rd b/man/GetTissueCoordinates.Rd
index e7bd6a516..5ccb69d68 100644
--- a/man/GetTissueCoordinates.Rd
+++ b/man/GetTissueCoordinates.Rd
@@ -36,3 +36,4 @@ Get Tissue Coordinates
 \code{\link[SeuratObject:GetTissueCoordinates]{SeuratObject::GetTissueCoordinates}}
 }
 \concept{objects}
+\concept{spatial}
diff --git a/man/HVFInfo.SCTAssay.Rd b/man/HVFInfo.SCTAssay.Rd
new file mode 100644
index 000000000..6e26995d8
--- /dev/null
+++ b/man/HVFInfo.SCTAssay.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/objects.R
+\name{HVFInfo.SCTAssay}
+\alias{HVFInfo.SCTAssay}
+\title{Get Variable Feature Information}
+\usage{
+\method{HVFInfo}{SCTAssay}(object, selection.method, status = FALSE, ...)
+}
+\arguments{
+\item{object}{An object}
+
+\item{selection.method}{Which method to pull. For \code{HVFInfo} and
+\code{VariableFeatures}, choose one from one of the
+following:
+\itemize{
+ \item \dQuote{vst}
+ \item \dQuote{sctransform} or \dQuote{sct}
+ \item \dQuote{mean.var.plot}, \dQuote{dispersion}, \dQuote{mvp}, or
+  \dQuote{disp}
+}
+For \code{SVFInfo} and \code{SpatiallyVariableFeatures}, choose from:
+\itemize{
+ \item \dQuote{markvariogram}
+ \item \dQuote{moransi}
+}}
+
+\item{status}{Add variable status to the resulting data frame}
+
+\item{...}{Arguments passed to other methods}
+}
+\description{
+Get variable feature information from \code{\link{SCTAssay}} objects
+}
+\examples{
+# Get the HVF info directly from an SCTAssay object
+pbmc_small <- SCTransform(pbmc_small)
+HVFInfo(pbmc_small[["SCT"]], selection.method = 'sct')[1:5, ]
+
+}
+\seealso{
+\code{\link[SeuratObject]{HVFInfo}}
+}
diff --git a/man/ISpatialDimPlot.Rd b/man/ISpatialDimPlot.Rd
index ff6c5360c..9f4f65b90 100644
--- a/man/ISpatialDimPlot.Rd
+++ b/man/ISpatialDimPlot.Rd
@@ -23,4 +23,5 @@ Returns final plot as a ggplot object
 \description{
 Visualize clusters spatially and interactively
 }
+\concept{spatial}
 \concept{visualization}
diff --git a/man/ISpatialFeaturePlot.Rd b/man/ISpatialFeaturePlot.Rd
index 9c175eb82..3e4c022c1 100644
--- a/man/ISpatialFeaturePlot.Rd
+++ b/man/ISpatialFeaturePlot.Rd
@@ -30,4 +30,5 @@ Returns final plot as a ggplot object
 \description{
 Visualize features spatially and interactively
 }
+\concept{spatial}
 \concept{visualization}
diff --git a/man/JackStrawPlot.Rd b/man/JackStrawPlot.Rd
index c1176f60b..4b3ce4c52 100644
--- a/man/JackStrawPlot.Rd
+++ b/man/JackStrawPlot.Rd
@@ -4,13 +4,25 @@
 \alias{JackStrawPlot}
 \title{JackStraw Plot}
 \usage{
-JackStrawPlot(object, dims = 1:5, reduction = "pca", xmax = 0.1, ymax = 0.3)
+JackStrawPlot(
+  object,
+  dims = 1:5,
+  cols = NULL,
+  reduction = "pca",
+  xmax = 0.1,
+  ymax = 0.3
+)
 }
 \arguments{
 \item{object}{Seurat object}
 
 \item{dims}{Dims to plot}
 
+\item{cols}{Vector of colors, each color corresponds to an individual PC. This may also be a single character
+or numeric value corresponding to a palette as specified by \code{\link[RColorBrewer]{brewer.pal.info}}.
+By default, ggplot2 assigns colors. We also include a number of palettes from the pals package.
+See \code{\link{DiscretePalette}} for details.}
+
 \item{reduction}{reduction to pull jackstraw info from}
 
 \item{xmax}{X-axis maximum on each QQ plot.}
diff --git a/man/LinkedPlots.Rd b/man/LinkedPlots.Rd
index f67d8f340..ff3bd2aa0 100644
--- a/man/LinkedPlots.Rd
+++ b/man/LinkedPlots.Rd
@@ -66,4 +66,5 @@ LinkedFeaturePlot(seurat.object, feature = 'Hpca')
 }
 
 }
+\concept{spatial}
 \concept{visualization}
diff --git a/man/PlotClusterTree.Rd b/man/PlotClusterTree.Rd
index ffcf41865..c88db7b3c 100644
--- a/man/PlotClusterTree.Rd
+++ b/man/PlotClusterTree.Rd
@@ -4,11 +4,14 @@
 \alias{PlotClusterTree}
 \title{Plot clusters as a tree}
 \usage{
-PlotClusterTree(object, ...)
+PlotClusterTree(object, direction = "downwards", ...)
 }
 \arguments{
 \item{object}{Seurat object}
 
+\item{direction}{A character string specifying the direction of the tree (default is downwards)
+Possible options: "rightwards", "leftwards", "upwards", and "downwards".}
+
 \item{\dots}{Additional arguments to
 \code{\link[ape:plot.phylo]{ape::plot.phylo}}}
 }
diff --git a/man/PolyFeaturePlot.Rd b/man/PolyFeaturePlot.Rd
index 56b441b52..1eacd0ecd 100644
--- a/man/PolyFeaturePlot.Rd
+++ b/man/PolyFeaturePlot.Rd
@@ -49,4 +49,5 @@ Returns a ggplot object
 \description{
 Plot cells as polygons, rather than single points. Color cells by any value accessible by \code{\link{FetchData}}.
 }
+\concept{spatial}
 \concept{visualization}
diff --git a/man/Radius.Rd b/man/Radius.Rd
index 6c0263d99..a0b1a74db 100644
--- a/man/Radius.Rd
+++ b/man/Radius.Rd
@@ -22,3 +22,4 @@ Get Spot Radius
 \code{\link[SeuratObject:Radius]{SeuratObject::Radius}}
 }
 \concept{objects}
+\concept{spatial}
diff --git a/man/Read10X.Rd b/man/Read10X.Rd
index a98d18338..c6691c460 100644
--- a/man/Read10X.Rd
+++ b/man/Read10X.Rd
@@ -5,7 +5,7 @@
 \title{Load in data from 10X}
 \usage{
 Read10X(
-  data.dir = NULL,
+  data.dir,
   gene.column = 2,
   cell.column = 1,
   unique.features = TRUE,
diff --git a/man/ReadMtx.Rd b/man/ReadMtx.Rd
new file mode 100644
index 000000000..e24620851
--- /dev/null
+++ b/man/ReadMtx.Rd
@@ -0,0 +1,63 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/preprocessing.R
+\name{ReadMtx}
+\alias{ReadMtx}
+\title{Load in data from remote or local mtx files}
+\usage{
+ReadMtx(
+  mtx,
+  cells,
+  features,
+  cell.column = 1,
+  feature.column = 2,
+  skip.cell = 0,
+  skip.feature = 0,
+  unique.features = TRUE,
+  strip.suffix = FALSE
+)
+}
+\arguments{
+\item{mtx}{Name or remote URL of the mtx file}
+
+\item{cells}{Name or remote URL of the cells/barcodes file}
+
+\item{features}{Name or remote URL of the features/genes file}
+
+\item{cell.column}{Specify which column of cells file to use for cell names; default is 1}
+
+\item{feature.column}{Specify which column of features files to use for feature/gene names; default is 2}
+
+\item{skip.cell}{Number of lines to skip in the cells file before beginning to read cell names}
+
+\item{skip.feature}{Number of lines to skip in the features file before beginning to gene names}
+
+\item{unique.features}{Make feature names unique (default TRUE)}
+
+\item{strip.suffix}{Remove trailing "-1" if present in all cell barcodes.}
+}
+\value{
+A sparse matrix containing the expression data.
+}
+\description{
+Enables easy loading of sparse data matrices
+}
+\examples{
+\dontrun{
+# For local files:
+
+expression_matrix <- ReadMtx(
+  mtx = "count_matrix.mtx.gz", features = "features.tsv.gz",
+  cells = "barcodes.tsv.gz"
+)
+seurat_object <- CreateSeuratObject(counts = expression_matrix)
+
+# For remote files:
+
+expression_matrix <- ReadMtx(mtx = "http://localhost/matrix.mtx",
+cells = "http://localhost/barcodes.tsv",
+features = "http://localhost/genes.tsv")
+seurat_object <- CreateSeuratObject(counts = data)
+}
+
+}
+\concept{preprocessing}
diff --git a/man/RenameCells.Rd b/man/RenameCells.Rd
index f95dc7b57..105b7d1e9 100644
--- a/man/RenameCells.Rd
+++ b/man/RenameCells.Rd
@@ -28,3 +28,4 @@ Rename Cells in an Object
 \seealso{
 \code{\link[SeuratObject:RenameCells]{SeuratObject::RenameCells}}
 }
+\concept{objects}
diff --git a/man/RunLDA.Rd b/man/RunLDA.Rd
index a7fd0dccd..22ef6635f 100644
--- a/man/RunLDA.Rd
+++ b/man/RunLDA.Rd
@@ -23,7 +23,7 @@ RunLDA(object, ...)
 
 \method{RunLDA}{Assay}(
   object,
-  assay,
+  assay = NULL,
   labels,
   features = NULL,
   verbose = TRUE,
diff --git a/man/SCTAssay-class.Rd b/man/SCTAssay-class.Rd
index 1cdc600b8..347894fbb 100644
--- a/man/SCTAssay-class.Rd
+++ b/man/SCTAssay-class.Rd
@@ -93,3 +93,4 @@ levels(pbmc_small[['SCT']])
 
 \code{\link{Assay}}
 }
+\concept{objects}
diff --git a/man/SCTResults.Rd b/man/SCTResults.Rd
index 181843eab..6092ca259 100644
--- a/man/SCTResults.Rd
+++ b/man/SCTResults.Rd
@@ -48,3 +48,4 @@ the slot directly).
 Pull the \code{\link{SCTResults}} information from an \code{\link{SCTAssay}}
 object.
 }
+\concept{objects}
diff --git a/man/STARmap-class.Rd b/man/STARmap-class.Rd
index 66b744ef9..7984f21c3 100644
--- a/man/STARmap-class.Rd
+++ b/man/STARmap-class.Rd
@@ -10,12 +10,6 @@ The STARmap class
 }
 \section{Slots}{
 
-\describe{
-\item{\code{...}}{}
-}}
-
-\section{Slots}{
-
 
 \describe{
 \item{\code{assay}}{Name of assay to associate image data with; will give this image
@@ -26,3 +20,5 @@ in a \code{Seurat} object}
 }
 }
 
+\concept{objects}
+\concept{spatial}
diff --git a/man/ScaleFactors.Rd b/man/ScaleFactors.Rd
index 93744a61a..bac4a5fee 100644
--- a/man/ScaleFactors.Rd
+++ b/man/ScaleFactors.Rd
@@ -37,3 +37,4 @@ Get image scale factors
 \code{scalefactors} objects can be created with \code{scalefactors()}
 }
 \concept{objects}
+\concept{spatial}
diff --git a/man/SlideSeq-class.Rd b/man/SlideSeq-class.Rd
index 1c3231f31..a81a02e3c 100644
--- a/man/SlideSeq-class.Rd
+++ b/man/SlideSeq-class.Rd
@@ -12,8 +12,6 @@ The SlideSeq class represents spatial information from the Slide-seq platform
 
 \describe{
 \item{\code{coordinates}}{...}
-
-\item{\code{...}}{}
 }}
 
 \section{Slots}{
@@ -28,3 +26,4 @@ in a \code{Seurat} object}
 }
 }
 
+\concept{spatial}
diff --git a/man/SpatialPlot.Rd b/man/SpatialPlot.Rd
index dd13d6ea9..d131cd820 100644
--- a/man/SpatialPlot.Rd
+++ b/man/SpatialPlot.Rd
@@ -175,4 +175,5 @@ SpatialDimPlot(seurat.object, group.by = "clusters")
 
 }
 \concept{convenience}
+\concept{spatial}
 \concept{visualization}
diff --git a/man/UpdateSCTAssays.Rd b/man/UpdateSCTAssays.Rd
index 611120589..92ae5edd0 100644
--- a/man/UpdateSCTAssays.Rd
+++ b/man/UpdateSCTAssays.Rd
@@ -17,3 +17,4 @@ A Seurat object with updated SCTAssays
 Update pre-V4 Assays generated with SCTransform in the Seurat to the new
 SCTAssay class
 }
+\concept{objects}
diff --git a/man/VisiumV1-class.Rd b/man/VisiumV1-class.Rd
index e42bd19be..bc2bd3a87 100644
--- a/man/VisiumV1-class.Rd
+++ b/man/VisiumV1-class.Rd
@@ -24,3 +24,4 @@ platform
 }}
 
 \concept{objects}
+\concept{spatial}
diff --git a/man/as.sparse.Rd b/man/as.sparse.Rd
index b0cdaae48..1872e4fe8 100644
--- a/man/as.sparse.Rd
+++ b/man/as.sparse.Rd
@@ -20,8 +20,8 @@
 
 \item{...}{Arguments passed to other methods}
 
-\item{row.names}{\code{NULL} or a character vector giving the row
-    names for the data frame.  Missing values are not allowed.}
+\item{row.names}{\code{NULL} or a character vector giving the row names for
+the data; missing values are not allowed}
 
 \item{optional}{logical. If \code{TRUE}, setting row names and
     converting column names (to syntactic names: see
diff --git a/man/merge.SCTAssay.Rd b/man/merge.SCTAssay.Rd
index 6c66cc8af..f976cbc1d 100644
--- a/man/merge.SCTAssay.Rd
+++ b/man/merge.SCTAssay.Rd
@@ -34,3 +34,4 @@ populated with NAs.}
 \description{
 Merge SCTAssay objects
 }
+\concept{objects}
diff --git a/man/subset.AnchorSet.Rd b/man/subset.AnchorSet.Rd
index 0bc058e19..6f5adf634 100644
--- a/man/subset.AnchorSet.Rd
+++ b/man/subset.AnchorSet.Rd
@@ -45,3 +45,4 @@ filtered out
 \description{
 Subset an AnchorSet object
 }
+\concept{objects}
diff --git a/tests/testthat/test_differential_expression.R b/tests/testthat/test_differential_expression.R
index d8ed8ce57..81e231f76 100644
--- a/tests/testthat/test_differential_expression.R
+++ b/tests/testthat/test_differential_expression.R
@@ -56,7 +56,7 @@ test_that("features parameter behaves correctly ", {
 
 results <- suppressWarnings(FindMarkers(object = pbmc_small, ident.1 = Cells(x = pbmc_small)[1:40], ident.2 = Cells(x = pbmc_small)[41:80], verbose = FALSE, base = exp(1)))
 test_that("passing cell names works", {
-  expect_equal(nrow(x = results), 176)
+  expect_equal(nrow(x = results), 190)
   expect_equal(results[1, "p_val"], 0.0001690882)
   expect_equal(results[1, "avg_logFC"], -1.790824, tolerance = 1e-6)
   expect_equal(results[1, "pct.1"], 0.075)
@@ -78,18 +78,18 @@ test_that("logfc.threshold warns when none met", {
 
 results <- suppressWarnings(FindMarkers(object = pbmc_small, ident.1 = 0, ident.2 = 1, min.pct = 0.5, verbose = FALSE, base = exp(1)))
 test_that("min.pct works", {
-  expect_equal(nrow(x = results), 63)
+  expect_equal(nrow(x = results), 65)
   expect_gte(min(apply(X = results, MARGIN = 1, FUN = function(x) max(x[3], x[4]))), 0.5)
 })
 
-results <- expect_warning(FindMarkers(object = pbmc_small, ident.1 = 0, ident.2 = 1, min.pct = 1.0, verbose = FALSE, base = exp(1)))
+results <- expect_warning(FindMarkers(object = pbmc_small, ident.1 = 0, ident.2 = 1, min.pct = 2.0, verbose = FALSE, base = exp(1)))
 test_that("min.pct warns when none met", {
   expect_equal(nrow(x = results), 0)
 })
 
 results <- suppressWarnings(FindMarkers(object = pbmc_small, ident.1 = 0, ident.2 = 1, min.diff.pct = 0.5, verbose = FALSE, base = exp(1)))
 test_that("min.diff.pct works", {
-  expect_equal(nrow(x = results), 43)
+  expect_equal(nrow(x = results), 44)
   expect_gte(min(apply(X = results, MARGIN = 1, FUN = function(x) abs(x[4] - x[3]))), 0.5)
 })
 
@@ -150,7 +150,7 @@ Idents(object = t2) <- "groups"
 results2 <- suppressWarnings(FindMarkers(object = t2, ident.1 = "g1", ident.2 = "g2", verbose = FALSE, base = exp(1)))
 
 test_that("subset.ident works", {
-  expect_equal(nrow(x = results), 114)
+  expect_equal(nrow(x = results), 127)
   expect_equal(results, results2)
   expect_equal(results[1, "p_val"], 0.01293720)
   expect_equal(results[1, "avg_logFC"], 1.799280, tolerance = 1e-6)
@@ -261,7 +261,7 @@ if (requireNamespace('metap', quietly = TRUE)) {
     expect_equal(markers[1, "g1_p_val_adj"], 9.077279e-06)
     expect_equal(markers[1, "max_pval"], 4.983576e-05)
     expect_equal(markers[1, "minimump_p_val"], 7.893286e-08)
-    expect_equal(nrow(markers), 162)
+    expect_equal(nrow(markers), 179)
     expect_equal(rownames(markers)[1], "HLA-DRB1")
     expect_equal(markers[, "max_pval"], unname(obj = apply(X = markers, MARGIN = 1, FUN = function(x) max(x[c("g1_p_val", "g2_p_val")]))))
   })
@@ -286,7 +286,7 @@ if (requireNamespace('metap', quietly = TRUE)) {
     # expect_equal(markers.missing[1, "g2_pct.1"], 0.062)
     expect_equal(markers.missing[1, "g2_pct.2"], 0.95)
     expect_equal(markers.missing[1, "g2_p_val_adj"], 3.847695e-11)
-    expect_equal(nrow(markers.missing), 190)
+    expect_equal(nrow(markers.missing), 205)
     expect_equal(rownames(markers.missing)[1], "HLA-DPB1")
   })
 }
diff --git a/tests/testthat/test_integration.R b/tests/testthat/test_integration.R
index f3f65c310..ff5f57e01 100644
--- a/tests/testthat/test_integration.R
+++ b/tests/testthat/test_integration.R
@@ -194,10 +194,10 @@ test_that("FindTransferAnchors with project.query and reference.reduction works"
 ref <- FindNeighbors(object = ref, reduction = "pca", dims = 1:30, return.neighbor = TRUE, k.param = 31, verbose = FALSE, l2.norm = TRUE, nn.method = "annoy")
 test_that("FindTransferAnchors with reference.neighbors precomputed works", {
   skip_on_cran()
-  anchors <- FindTransferAnchors(reference = ref, query = query, reference.neighbors = "RNA_nn", k.filter = 50)
+  anchors <- FindTransferAnchors(reference = ref, query = query, reference.neighbors = "RNA.nn", k.filter = 50)
   expect_error(FindTransferAnchors(reference = ref, query = query, reference.neighbors = "BAD", k.filter = 50))
-  expect_error(FindTransferAnchors(reference = ref, query = query, reference.neighbors = "RNA_nn", k.filter = 50, k.score = 31))
-  expect_error(FindTransferAnchors(reference = ref, query = query, reference.neighbors = "RNA_nn", k.filter = 50, k.anchor = 31))
+  expect_error(FindTransferAnchors(reference = ref, query = query, reference.neighbors = "RNA.nn", k.filter = 50, k.score = 31))
+  expect_error(FindTransferAnchors(reference = ref, query = query, reference.neighbors = "RNA.nn", k.filter = 50, k.anchor = 31))
   co <- anchors@object.list[[1]]
   expect_equal(dim(co), c(100, 160))
   expect_equal(Reductions(co), c("pcaproject", "pcaproject.l2"))
diff --git a/tests/testthat/test_read_mtx.R b/tests/testthat/test_read_mtx.R
new file mode 100644
index 000000000..75707decc
--- /dev/null
+++ b/tests/testthat/test_read_mtx.R
@@ -0,0 +1,30 @@
+context("ReadMtx")
+
+test_that("skip.cell and skip.feature work", {
+  skip_on_cran()
+  mtx <- "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE126nnn/GSE126836/suppl/GSE126836_SN_MD5828_matrix.mtx.gz"
+  features <- "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE126nnn/GSE126836/suppl/GSE126836_SN_MD5828_genes.csv.gz"
+  cells <- "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE126nnn/GSE126836/suppl/GSE126836_SN_MD5828_barcodes.csv.gz"
+  counts1 <- ReadMtx(mtx = mtx, cells = cells, features = features, feature.column = 1, skip.cell = 1, skip.feature = 1)
+  expect_is(counts1, "dgCMatrix")
+  expect_equal(ncol(counts1), 1436)
+  expect_equal(nrow(counts1), 29445)
+  expect_equal(colnames(counts1)[5], "MD5828a_GGGCATCCAATGAAAC-1")
+  expect_equal(rownames(counts1)[2], "A1BG-AS1")
+})
+
+
+test_that("ReadMtx works", {
+  skip_on_cran()
+  mtx <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE127774&format=file&file=GSE127774%5FACC%5FB%5Fmatrix%2Emtx%2Egz"
+  cells <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE127774&format=file&file=GSE127774%5FACC%5FB%5Fbarcodes%2Etsv%2Egz"
+  features <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE127774&format=file&file=GSE127774%5FACC%5FB%5Fgenes%2Etsv%2Egz"
+  counts2 <- ReadMtx(mtx = mtx, cells = cells, features = features, feature.column = 1)
+  expect_is(counts2, "dgCMatrix")
+  expect_equal(ncol(counts2), 22063)
+  expect_equal(nrow(counts2), 22530)
+  expect_equal(colnames(counts2)[1], "AAACCTGAGCAATCTC-1")
+  expect_equal(rownames(counts2)[2], "ENSPPAG00000040697")
+})
+
+