fastMNN: Fast mutual nearest neighbors correction

In LTLA/batchelor: Single-Cell Batch Correction Methods

Fast mutual nearest neighbors correction

Description

Correct for batch effects in single-cell expression data using a fast version of the mutual nearest neighbors (MNN) method.

Usage

fastMNN(
  ...,
  batch = NULL,
  k = 20,
  prop.k = NULL,
  restrict = NULL,
  cos.norm = TRUE,
  ndist = 3,
  d = 50,
  deferred = TRUE,
  weights = NULL,
  get.variance = FALSE,
  merge.order = NULL,
  auto.merge = FALSE,
  min.batch.skip = 0,
  subset.row = NULL,
  correct.all = FALSE,
  assay.type = "logcounts",
  BSPARAM = IrlbaParam(),
  BNPARAM = KmknnParam(),
  BPPARAM = SerialParam()
)

Arguments

...	One or more log-expression matrices where genes correspond to rows and cells correspond to columns. Alternatively, one or more SingleCellExperiment objects can be supplied containing a log-expression matrix in the assay.type assay. Each object should contain the same number of rows, corresponding to the same genes in the same order. Objects of different types can be mixed together. If multiple objects are supplied, each object is assumed to contain all and only cells from a single batch. If a single object is supplied, it is assumed to contain cells from all batches, so batch should also be specified. Alternatively, one or more lists of matrices or SingleCellExperiments can be provided; this is flattened as if the objects inside each list were passed directly to ....
batch	A vector or factor specifying the batch of origin for all cells when only a single object is supplied in .... This is ignored if multiple objects are present.
k	An integer scalar specifying the number of nearest neighbors to consider when identifying MNNs.
prop.k	A numeric scalar in (0, 1) specifying the proportion of cells in each dataset to use for mutual nearest neighbor searching. If set, the number of nearest neighbors used for the MNN search in each batch is redefined as max(k, prop.k*N) where N is the number of cells in that batch.
restrict	A list of length equal to the number of objects in .... Each entry of the list corresponds to one batch and specifies the cells to use when computing the correction.
cos.norm	A logical scalar indicating whether cosine normalization should be performed on the input data prior to PCA.
ndist	A numeric scalar specifying the threshold beyond which neighbours are to be ignored when computing correction vectors. Each threshold is defined as a multiple of the number of median distances.
d	Numeric scalar specifying the number of dimensions to use for dimensionality reduction in multiBatchPCA. If NA, no dimensionality reduction is performed and any matrices in ... are used as-is.
deferred	Logical scalar indicating whether to defer centering/scaling, see multiBatchPCA for details.
weights, get.variance	Further arguments to pass to multiBatchPCA.
merge.order	An integer vector containing the linear merge order of batches in .... Alternatively, a list of lists representing a tree structure specifying a hierarchical merge order.
auto.merge	Logical scalar indicating whether to automatically identify the “best” merge order.
min.batch.skip	Numeric scalar specifying the minimum relative magnitude of the batch effect, below which no correction will be performed at a given merge step.
subset.row	A vector specifying which features to use for correction.
correct.all	Logical scalar indicating whether corrected expression values should be computed for genes not in subset.row. Only relevant if subset.row is not NULL.
assay.type	A string or integer scalar specifying the assay containing the log-expression values. Only used for SingleCellExperiment inputs.
BSPARAM	A BiocSingularParam object specifying the algorithm to use for PCA in multiBatchPCA.
BNPARAM	A BiocNeighborParam object specifying the nearest neighbor algorithm.
BPPARAM	A BiocParallelParam object specifying whether the PCA and nearest-neighbor searches should be parallelized.

Details

This function provides a variant of the mnnCorrect function, modified for speed and more robust performance. In particular:

• It performs a multi-sample PCA via multiBatchPCA and subsequently performs all calculations in the PC space. This reduces computational work and provides some denoising for improved neighbour detection. As a result, though, the corrected output cannot be interpreted on a gene level and is useful only for cell-level comparisons, e.g., clustering and visualization.

• The correction vector for each cell is directly computed from its k nearest neighbours in the same batch. Specifically, only the k nearest neighbouring cells that also participate in MNN pairs are used. Each MNN-participating neighbour is weighted by distance from the current cell, using a tricube scheme with bandwidth equal to the median distance multiplied by ndist. This ensures that the correction vector only uses information from the closest cells, improving the fidelity of local correction.

• Issues with “kissing” are avoided with a two-step procedure that removes variation along the batch effect vector. First, the average correction vector across all MNN pairs is computed. Cell coordinates are adjusted such that all cells in a single batch have the same position along this vector. The correction vectors are then recalculated with the adjusted coordinates (but the same MNN pairs).

The default setting of cos.norm=TRUE provides some protection against differences in scaling between log-expression matrices from batches that are normalized separately (see cosineNorm for details). However, if possible, we recommend using the output of multiBatchNorm as input to fastMNN. This will equalize coverage on the count level before the log-transformation, which is a more accurate rescaling than cosine normalization on the log-values.

The batch argument allows users to easily perform batch correction when all cells have already been combined into a single object. This avoids the need to manually split the matrix or SingleCellExperiment object into separate objects for input into fastMNN. In this situation, the order of input batches is defined by the order of levels in batch.

Value

A SingleCellExperiment is returned where each row is a gene and each column is a cell. This contains:

• A corrected matrix in the reducedDims slot, containing corrected low-dimensional coordinates for each cell. This has number of columns equal to d and number of rows equal to the total number of cells in ....

• A batch column in the colData slot, containing the batch of origin for each row (i.e., cell) in corrected.

• A rotation column the rowData slot, containing the rotation matrix used for the PCA. This has d columns and number of rows equal to the number of genes to report (see the “Choice of genes” section).

• A reconstructed matrix in the assays slot, containing the low-rank reconstruction of the expression matrix. This can be interpreted as per-gene corrected log-expression values (after cosine normalization, if cos.norm=TRUE) but should not be used for quantitative analyses. This has number of rows equal to the number of input genes if subset.row=NULL or correct.all=TRUE, otherwise each row corresponds to a gene in subset.row.

Cells in the output object are always ordered in the same manner as supplied in .... For a single input object, cells will be reported in the same order as they are arranged in that object. In cases with multiple input objects, the cell identities are simply concatenated from successive objects, i.e., all cells from the first object (in their provided order), then all cells from the second object, and so on. This is true regardless of the value of merge.order and auto.merge, which only affects the internal merge order of the batches.

Additionally, the metadata of the output object contains:

• merge.info, a DataFrame of diagnostic information about each merge step. See the “Merge diagnostics” section for more details.

• pca.info, a list of metadata produced by multiBatchPCA, such as the variance explained when get.variance=TRUE.

Controlling the merge order

By default, batches are merged in the user-supplied order in ..., i.e., the first batch is merged with the second batch, the third batch is merged with the combined first-second batch, the fourth batch is merged with the combined first-second-third batch and so on. We refer to this approach as a progressive merge. When batch is supplied for a single object in ..., the ordering of batches in a progressive merge is determined by the ordering of factor levels in batch.

If merge.order is an integer vector, it is treated as an ordering permutation with which to perform a progressive merge. For example, if merge.order=c(4,1,3,2), batches 4 and 1 in ... are merged first; batch 3 is merged with the combined 4+1 batch; and then batch 2 is merged with the combined 4+1+3 batch. This is often more convenient than changing the order manually in ..., which would alter the order of batches in the output corrected matrix.

If merge.order is a character vector, it is treated as an ordering permutation for named batches. If merge.order is a factor, it is coerced into a character vector and also treated as a permutation of names.

If merge.order is a nested list, it is treated as a tree that specifies a hierarchical merge. Each element of the list should either be a string or integer scalar, corresponding to a leaf node that specifies a batch; or another list, corresponding to an internal node that should contain at least two children; or an integer or character vector of length 2 or more, again corresponding to an internal node.

• For example, list(list(1,2), list(3,4)) indicates that batch 1 should be merged with batch 2; batch 3 should be merged with batch 4; and that, finally, the combined batches 1+2 and 3+4 should be merged.

• More than two children per node are supported and will result in a progressive merge within that node. For example, list(list(1,2,3), list(4,5,6)) will merge batch 1 with 2, then 1+2 with 3; batch 4 with 5, and then 4+5 with 6; and finally, 1+2+3 with 4+5+6.

• The same approach can be used for integer or character vectors, e.g., list(1:3, 4:6) has the same effect as above.

Note that, while batches can be specified by name (character) or index (integer), users cannot use both in the same tree. Factors can be provided and are coerced to character vectors.

The merge order may occasionally be important as it determines the number of MNN pairs available at each merge step. MNN pairs results in greater stability of the batch vectors and increased likelihood of identifying shared subpopulations, which are important to the precision and accuracy of the MNN-based correction, respectively.

• In a progressive merge, the reference increases in size at each step, ensuring that more cells are available to identify MNN pairs in later merges. We suggest setting the largest, most heterogeneous batch as the first reference, which favors detection of sufficient MNN pairs between the first and other batches. Conversely, if two small batches without shared populations are supplied first, the wrong MNN pairs will be detected and the result of the merge will be incorrect.

• A merge tree is useful for merging together batches that are known to be more closely related (e.g., replicates) before attempting difficult merges involving more dissimilar batches. The idea is to increase the number of cells and thus MNN pairs prior to merging batches with few shared subpopulations. By comparison, performing the more difficult merges first is more likely to introduce errors whereby distinct subpopulations are incorrectly placed together, which is propagated to later steps as the initial merge is used as a reference for subsequent merges.

• If auto.merge=TRUE, merge steps are chosen to maximize the number of MNN pairs at each step. The aim is to improve the stability of the correction by first merging more similar batches with more MNN pairs. This can be somewhat time-consuming as MNN pairs need to be iteratively recomputed for all possible batch pairings.

The order of cells in the output is never affected by the setting of merge.order. It depends only on the order of objects in ... and the order of cells within each object.

Choice of genes

All genes are used with the default setting of subset.row=NULL. Users can set subset.row to subset the inputs to highly variable genes or marker genes. This improves the quality of the PCA and identification of MNN pairs by reducing the noise from irrelevant genes. Users should not be too restrictive with subsetting, as high dimensionality is required to satisfy the orthogonality assumption in MNN detection.

By default, only the selected genes are used to compute rotation vectors and a low-rank corrected expression matrix. However, setting correct.all=TRUE will return rotation vectors spanning all genes in the supplied input data. This ensures that corrected values are returned in reconstructed for all input genes, e.g., in correctExperiments. Note that this setting will not alter the corrected low-dimension coordinates, nor the rotation values for the selected genes.

If d=NA, any gene not in subset.row will have reconstructed values of zero when correct.all=TRUE. Without a PCA, we cannot easily extrapolate the correction to other genes.

Using restriction

See ?"batchelor-restrict" for a description of the restrict argument. Specifically, fastMNN will only use the restricted subset of cells in each batch to identify MNN pairs and the center of the orthogonalization. It will then extrapolate the correction to all cells in each batch.

Note that all cells are used to perform the PCA, regardless of whether restrict is set. This is generally desirable in applications where restrict is useful. For example, constructing the projection vectors with only control cells will not guarantee resolution of unique non-control populations in each batch.

However, this also means that the corrected values for the restricted cells will differ from the output when the inputs are directly subsetted to only contain the restricted cells. If this is not desirable, users can perform the PCA manually and apply reducedMNN instead.

Merge diagnostics

We can consider fastMNN's operation in terms of pairwise merge steps. Each merge step involves two mutually exclusive sets of cells, a “left” set and “right” set. Each set may consist of cells from different batches if those batches were merged in a previous step. The merge will then create a new set of cells that combines the left and right sets. Iteratively repeating this process with the newly formed sets will eventually merge all batches together.

The output metadata contains merge.info, a DataFrame where each row corresponds to a merge step. It contains the following fields:

• left, a List of integer or character vectors. Each vector specifies the batches in the left set at a given merge step.

• right, a similar List of integer or character vectors. Each vector specifies the batches in the right set at a given merge step.

• pairs, a List of DataFrames specifying which pairs of cells were identified as MNNs at each step. In each DataFrame, each row corresponds to a single MNN pair and specifies the paired cells that were in the left and right sets, respectively. Note that the indices refer to those paired cells in the output ordering of cells, i.e., users can identify the paired cells at each step by column-indexing the output of the fastMNN function.

• batch.size, a numeric vector specifying the relative magnitude of the batch effect at each merge, see “Orthogonalization details”.

• skipped, a logical vector indicating whether the correction was skipped if the magnitude of the batch effect was below min.batch.skip.

• lost.var, a numeric matrix specifying the percentage of variance lost due to orthogonalization at each merge step. This is reported separately for each batch (columns, ordered according to the input order, not the merge order).

Specifying the number of neighbors

The threshold to define nearest neighbors is defined by k, which is passed to findMutualNN to identify MNN pairs. The size of k can be roughly interpreted as the anticipated minimum size of a shared subpopulation in each batch. If a batch has fewer than k cells of a shared subpopulation, there is an increased risk that its counterparts in other batches will form incorrect MNN pairs.

From the perspective of the algorithm, larger values allow for more MNN pairs to be obtained, which improves the stability of the correction vectors. Larger values also increase robustness against non-orthogonality, by ignoring a certain level of biological variation when identifying pairs. This can be used to avoid the kissing problem where MNN pairs are only detected on the “surface” of the distribution. However, values of k should not be too large, as this would result in MNN pairs being inappropriately identified between biologically distinct populations.

In practice, increasing k will generally result in more aggressive merging as the algorithm is more generous in matching subpopulations across batches. We suggest starting with the default k and increasing it if one is confident that the same cell types are not adequately merged across batches. This is better than starting with a large k as incorrect merging is much harder to diagnose than insufficient merging.

An additional consideration is that the effect of any given k will vary with the number of cells in each batch. With more cells, a larger k may be preferable to achieve better merging in the presence of non-orthogonality. We can achieve this by setting prop.k, e.g., prop.k=0.05 will set k to 5% of the number of cells in each batch. This allows the choice of k to adapt to the size of each batch at each merge step and handles asymmetry in batch sizes (via the k1 and k2 arguments in findMutualNN).

Orthogonalization details

fastMNN will compute the percentage of variance that is lost from each batch during orthogonalization at each merge step. This represents the variance in each batch that is parallel to the average correction vectors (and hence removed during orthogonalization) at each merge step. Large proportions suggest that there is biological structure that is parallel to the batch effect, corresponding to violations of the assumption that the batch effect is orthogonal to the biological subspace. A rule of thumb is that more than 10% of lost variance is cause for closer examination, though this is highly dependent on the context, e.g., a large lost proportion may be fine if the population structure is still approximately preserved while a small lost proportion may be unacceptable if an important rare subpopulation can no longer be distinguished.

In rare cases, orthogonalization may cause problems in the absence of a batch effect, resulting in large losses of variance. To avoid this, fastMNN will not perform any correction if the relative magnitude of the batch effect is less than min.batch.skip. The relative magnitude is defined as the L2 norm of the average correction vector divided by the root-mean-square of the L2 norms of the per-MNN pair correction vectors. This will be large when the per-pair vectors are all pointing in the same direction, and small when the per-pair vectors point in random directions due to the absence of a consistent batch effect. If a large loss of variance is observed along with a small batch effect in a given merge step, users can set min.batch.skip to simply skip correction in that step.

Author(s)

Aaron Lun

References

Haghverdi L, Lun ATL, Morgan MD, Marioni JC (2018). Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nat. Biotechnol. 36(5):421

Lun ATL (2018). Further MNN algorithm development. https://MarioniLab.github.io/FurtherMNN2018/theory/description.html

See Also

cosineNorm and multiBatchPCA, to obtain the values to be corrected.

reducedMNN, for a version of the function that operates in low-dimensional space.

mnnCorrect, for the “classic” version of the MNN correction algorithm.

clusterMNN, for the cluster-based version of this approach.

multiModalMNN from the mumosa package, which extends this to data in the altExps.

mnnDeltaVariance, to obtain further per-gene diagnostics on the behavior of the correction.

Examples

B1 <- matrix(rnorm(10000, -1), ncol=50) # Batch 1 
B2 <- matrix(rnorm(10000, 1), ncol=50) # Batch 2
out <- fastMNN(B1, B2)

# Corrected values for use in clustering, etc.
str(reducedDim(out)) 

# Extracting corrected expression values for gene 10.
summary(assay(out)[10,])

LTLA/batchelor documentation built on July 10, 2024, 9:09 p.m.