# Extending dispatch to more batch correction methods

Aaron Lun1

1Cancer Research UK Cambridge Institute, Cambridge, United Kingdom

#### Revised: 11 February 2019

#### Package

batchelor 1.2.4

# 1 Overview

The batchelor package provides the `batchCorrect()` generic that dispatches on its `PARAM` argument. Users writing code using `batchCorrect()` can easily change one method for another by simply modifying the class of object supplied as `PARAM`. For example:

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

    # Switching easily between batch correction methods.
    m.out <- batchCorrect(B1, B2, PARAM=ClassicMnnParam())
    f.out <- batchCorrect(B1, B2, PARAM=FastMnnParam(d=20))
    r.out <- batchCorrect(B1, B2, PARAM=RescaleParam(pseudo.count=0))

Developers of other packages can extend this further by adding their batch correction methods to this dispatch system. This improves interoperability across packages by allowing users to easily experiment with different methods.

# 2 Setting up

You will need to `Imports: batchelor, methods` in your `DESCRIPTION` file. You will also need to add `import(methods)`, `importFrom(batchelor, "batchCorrect")` and `importClassesFrom(batchelor, "BatchelorParam")` in your `NAMESPACE` file.

Obviously, you will also need to have a function that implements your batch correction method. For demonstration purposes, we will use an identity function that simply returns the input values1 1 Not a very good correction, but that’s not the point right now.. This is implemented like so:

    noCorrect <- function(...)
    # Takes a set of batches and returns them without modification.
    {
       do.call(cbind, list(...))
    }

# 3 Deriving a `BatchelorParam` subclass

We need to define a new `BatchelorParam` subclass that instructs the `batchCorrect()` generic to dispatch to our new method. This is most easily done like so:

    NothingParam <- setClass("NothingParam", contains="BatchelorParam")

Note that `BatchelorParam` itself is derived from a `SimpleList` and can be modified with standard list operators like `$`.

    nothing <- NothingParam()
    nothing

    ## NothingParam of length 0

    nothing$some_value <- 1
    nothing

    ## NothingParam of length 1
    ## names(1): some_value

If no parameters are set, the default values in the function will be used2 2 Here there are none in `noCorrect()`, but presumably your function is more complex than that.. Additional slots can be specified in the class definition if there are important parameters that need to be manually specified by the user.

# 4 Defining a `batchCorrect` method

## 4.1 Input

The `batchCorrect()` generic looks like this:

    batchCorrect

    ## standardGeneric for "batchCorrect" defined from package "batchelor"
    ##
    ## function (..., batch = NULL, restrict = NULL, subset.row = NULL,
    ##     correct.all = FALSE, assay.type = NULL, get.spikes = FALSE,
    ##     PARAM)
    ## standardGeneric("batchCorrect")
    ## <bytecode: 0x0000000024c847b8>
    ## <environment: 0x0000000024c57dd0>
    ## Methods may be defined for arguments: PARAM
    ## Use  showMethods("batchCorrect")  for currently available ones.

Any implemented method must accept one or more matrix-like objects containing single-cell gene expression matrices in `...`. Rows are assumed to be genes and columns are assumed to be cells. If only one object is supplied, `batch` must be specified to indicate the batches to which each cell belongs.

Alternatively, one or more `SingleCellExperiment` objects can be supplied, containing the gene expression matrix in the `assay.type` assay. These should not be mixed with matrix-like objects, i.e., if one object is a `SingleCellExperiment`, all objects should be `SingleCellExperiment`s.

The `subset.row=` argument specifies a subset of genes on which to perform the correction. The `get.spikes=` argument specifies whether spike-in transcripts should be used for `SingleCellExperiment` inputs. The `correct.all=` argument specifies whether corrected values should be returned for all genes, by “extrapolating” from the results for the genes that were used3 3 If your method cannot support this option, setting it to `TRUE` should raise an error.. See the Output section below for the expected output from each combination of these settings.

The `restrict=` argument allows users to compute the correction using only a subset of cells in each batch (e.g., from cell controls). The correction is then “extrapolated” to all cells in the batch4 4 Again, if your method cannot support this, any non-`NULL` value of `restrict` should raise an error., such that corrected values are returned for all cells.

## 4.2 Output

Any implemented method must return a `SingleCellExperiment` where the first assay contains corrected gene expression values for all genes. Corrected values should be returned:

- for all genes, if `correct.all=TRUE` for any input.
- for all genes, if `subset.row=NULL` for matrix-like inputs.
- for selected genes, if `subset.row` is not `NULL` for matrix-like inputs.
- for all genes, if `subset.row=NULL` and `get.spikes=TRUE` for `SingleCellExperiment` inputs.
- for all non-spike-in transcripts, if `subset.row=NULL` and `get.spikes=FALSE` for `SingleCellExperiment` inputs.
- for selected genes, if `subset.row` is not `NULL` and `get.spikes=FALSE` for matrix-like inputs.
- for selected genes that are not spike-in transcripts, if `subset.row` is not `NULL` and `get.spikes=FALSE` for `SingleCellExperiment` inputs.

This looks more complicated than it really is, as developers can just modify `subset.row` to exclude spike-ins when `get.spikes=FALSE`.

Cells should be reported in the same order that they are supplied in `...`. In cases with multiple batches, the cell identities are simply concatenated from successive objects in their specified order, i.e., all cells from the first object (in their provided order), then all cells from the second object, and so on. If there is only a single batch, the order of cells in that batch should be preserved.

The output object should have row names equal to the row names of the input objects. Column names should be equivalent to the concatenated column names of the input objects, unless all are `NULL`, in which case the column names in the output can be `NULL`. In situations where some input objects have column names, and others are `NULL`, those missing column names should be filled in with empty strings. This represents the expected behaviour when `cbind`ing multiple matrices together.

Finally, the `colData` slot should contain ‘batch’, a vector specifying the batch of origin for each cell.

## 4.3 Demonstration

Finally, we define a method that calls our `noCorrect` function while satisfying all of the above input/output requirements. To be clear, it is not mandatory to lay out the code as shown below; this is simply one way that all the requirements can be satisfied. We have used some internal batchelor functions for brevity - please contact us if you find these useful and want them to be exported.

    setMethod("batchCorrect", "NothingParam", function(..., batch = NULL,
        restrict=NULL, subset.row = NULL, correct.all = FALSE,
        assay.type = "logcounts", get.spikes = FALSE, PARAM)
    {
        batches <- list(...)
        checkBatchConsistency(batches)

        # Pulling out information from the SCE objects.
        is.sce <- checkIfSCE(batches)
        if (any(is.sce)) {
            sce.batches <- batches[is.sce]
            checkSpikeConsistency(sce.batches)
            subset.row <- batchelor:::.SCE_subset_genes(subset.row, sce.batches[[1]], get.spikes)
            batches[is.sce] <- lapply(sce.batches, assay, i=assay.type)
        }

        # Subsetting by 'batch', if only one object is supplied.
        do.split <- length(batches)==1L
        if (do.split) {
            divided <- divideIntoBatches(batches[[1]], batch=batch, restrict=restrict)
            batches <- divided$batches
            restrict <- divided$restricted
        }

        # Subsetting by row.
        # This is a per-gene "method", so correct.all=TRUE will ignore subset.row.
        # More complex methods will need to handle this differently.
        if (correct.all) {
            subset.row <- NULL
        } else if (!is.null(subset.row)) {
            subset.row <- normalizeSingleBracketSubscript(originals[[1]], subset.row)
            batches <- lapply(batches, "[", i=subset.row, , drop=FALSE)
        }

        # Don't really need to consider restrict!=NULL here, as this function
        # doesn't do anything with the cells anyway.
        output <- do.call(noCorrect, batches)

        # Reordering the output for correctness if it was previously split.
        if (do.split) {
            d.reo <- divided$reorder
            output <- output[,d.reo,drop=FALSE]
        }

        ncells.per.batch <- vapply(batches, FUN=ncol, FUN.VALUE=0L)
        batch.names <- names(batches)
        if (is.null(batch.names)) {
            batch.names <- seq_along(batches)
        }

        SingleCellExperiment(list(corrected=output),
            colData=DataFrame(batch=rep(batch.names, ncells.per.batch)))
    })

And it works5 5 In a strictly programming sense, as the method itself does no correction at all.:

    n.out <- batchCorrect(B1, B2, PARAM=NothingParam())
    n.out

    ## class: SingleCellExperiment
    ## dim: 200 100
    ## metadata(0):
    ## assays(1): corrected
    ## rownames: NULL
    ## rowData names(0):
    ## colnames: NULL
    ## colData names(1): batch
    ## reducedDimNames(0):
    ## spikeNames(0):
    ## altExpNames(0):

Remember to export both the new method and the `NothingParam` class and constructor.

# 5 Session information

    sessionInfo()

    ## R version 3.6.2 (2019-12-12)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] scater_1.14.6               ggplot2_3.2.1
    ##  [3] scran_1.14.5                scRNAseq_2.0.2
    ##  [5] batchelor_1.2.4             SingleCellExperiment_1.8.0
    ##  [7] SummarizedExperiment_1.16.1 DelayedArray_0.12.1
    ##  [9] BiocParallel_1.20.1         matrixStats_0.55.0
    ## [11] Biobase_2.46.0              GenomicRanges_1.38.0
    ## [13] GenomeInfoDb_1.22.0         IRanges_2.20.1
    ## [15] S4Vectors_0.24.1            BiocGenerics_0.32.0
    ## [17] knitr_1.26                  BiocStyle_2.14.2
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] bitops_1.0-6                  bit64_0.9-7
    ##  [3] httr_1.4.1                    tools_3.6.2
    ##  [5] backports_1.1.5               R6_2.4.1
    ##  [7] irlba_2.3.3                   vipor_0.4.5
    ##  [9] DBI_1.1.0                     lazyeval_0.2.2
    ## [11] colorspace_1.4-1              withr_2.1.2
    ## [13] tidyselect_0.2.5              gridExtra_2.3
    ## [15] bit_1.1-14                    curl_4.3
    ## [17] compiler_3.6.2                BiocNeighbors_1.4.1
    ## [19] labeling_0.3                  bookdown_0.16
    ## [21] scales_1.1.0                  rappdirs_0.3.1
    ## [23] stringr_1.4.0                 digest_0.6.23
    ## [25] rmarkdown_2.0                 XVector_0.26.0
    ## [27] pkgconfig_2.0.3               htmltools_0.4.0
    ## [29] limma_3.42.0                  fastmap_1.0.1
    ## [31] dbplyr_1.4.2                  rlang_0.4.2
    ## [33] RSQLite_2.1.5                 shiny_1.4.0
    ## [35] DelayedMatrixStats_1.8.0      farver_2.0.1
    ## [37] dplyr_0.8.3                   RCurl_1.95-4.12
    ## [39] magrittr_1.5                  BiocSingular_1.2.1
    ## [41] GenomeInfoDbData_1.2.2        Matrix_1.2-18
    ## [43] Rcpp_1.0.3                    ggbeeswarm_0.6.0
    ## [45] munsell_0.5.0                 viridis_0.5.1
    ## [47] lifecycle_0.1.0               edgeR_3.28.0
    ## [49] stringi_1.4.3                 yaml_2.2.0
    ## [51] zlibbioc_1.32.0               BiocFileCache_1.10.2
    ## [53] AnnotationHub_2.18.0          grid_3.6.2
    ## [55] blob_1.2.0                    dqrng_0.2.1
    ## [57] promises_1.1.0                ExperimentHub_1.12.0
    ## [59] crayon_1.3.4                  lattice_0.20-38
    ## [61] cowplot_1.0.0                 locfit_1.5-9.1
    ## [63] zeallot_0.1.0                 pillar_1.4.3
    ## [65] igraph_1.2.4.2                glue_1.3.1
    ## [67] BiocVersion_3.10.1            evaluate_0.14
    ## [69] BiocManager_1.30.10           vctrs_0.2.1
    ## [71] httpuv_1.5.2                  gtable_0.3.0
    ## [73] purrr_0.3.3                   assertthat_0.2.1
    ## [75] xfun_0.11                     mime_0.8
    ## [77] rsvd_1.0.2                    xtable_1.8-4
    ## [79] later_1.0.0                   viridisLite_0.3.0
    ## [81] tibble_2.1.3                  AnnotationDbi_1.48.0
    ## [83] beeswarm_0.2.3                memoise_1.1.0
    ## [85] statmod_1.4.32                interactiveDisplayBase_1.24.0