--- output: html_document bibliography: ref.bib --- # Correction diagnostics ## Motivation Ideally, batch correction would remove the differences between batches while preserving the heterogeneity within batches. In the corrected data, cells of the same type should be intermingled and indistinguishable even if they come from different batches, while cells of different types should remain well-separated. Unfortunately, we rarely have prior knowledge of the underlying types of the cells, making it difficult to unambiguously determine whether differences between batches represent geniune biology or incomplete correction. Indeed, it could be said that all correction methods are at least somewhat incorrect (Section \@ref(sacrificing-differences)), though that not preclude them from being useful. In this chapter, we will describe a few diagnostics that, when combined with biological context, can be used to identify potential problems with the correction. We will recycle the `mnn.out`, `tab.mnn` and `clusters.mnn` objects that were produced in Section \@ref(mnn-correction). For the sake of brevity, we will not reproduce the relevant code - see Chapter \@ref(integrating-datasets) for more details. ```r mnn.out ``` ``` ## class: SingleCellExperiment ## dim: 13431 6791 ## metadata(2): merge.info pca.info ## assays(1): reconstructed ## rownames(13431): ENSG00000239945 ENSG00000228463 ... ENSG00000198695 ## ENSG00000198727 ## rowData names(1): rotation ## colnames: NULL ## colData names(1): batch ## reducedDimNames(1): corrected ## mainExpName: NULL ## altExpNames(0): ``` ```r tab.mnn ``` ``` ## Batch ## Cluster 1 2 ## 1 337 606 ## 2 289 542 ## 3 152 181 ## 4 12 4 ## 5 517 467 ## 6 17 19 ## 7 313 661 ## 8 162 118 ## 9 11 56 ## 10 547 1083 ## 11 17 59 ## 12 16 58 ## 13 144 93 ## 14 67 191 ## 15 4 36 ## 16 4 8 ``` ## Mixing between batches The simplest way to quantify the degree of mixing across batches is to test each cluster for imbalances in the contribution from each batch [@buttner2019test]. This is done by applying Pearson's chi-squared test to each row of `tab.mnn` where the expected proportions under the null hypothesis proportional to the total number of cells per batch. Low $p$-values indicate that there are significant imbalances In practice, this strategy is most suited to experiments where the batches are technical replicates with identical population composition; it is usually too stringent for batches with more biological variation, where proportions can genuinely vary even in the absence of any batch effect. ```r library(batchelor) p.values <- clusterAbundanceTest(tab.mnn) p.values ``` ``` ## 1 2 3 4 5 6 7 8 ## 9.047e-02 3.093e-02 6.700e-03 2.627e-03 8.424e-20 2.775e-01 5.546e-05 2.274e-11 ## 9 10 11 12 13 14 15 16 ## 2.136e-04 5.480e-05 4.019e-03 2.972e-03 1.538e-12 3.936e-05 2.197e-04 7.172e-01 ``` We favor a more qualitative approach where we compute the variance in the log-normalized abundances across batches for each cluster. A highly variable cluster has large relative differences in cell abundance across batches; this may be an indicator for incomplete batch correction, e.g., if the same cell type in two batches was not combined into a single cluster in the corrected data. We can then focus our attention on these clusters to determine whether they might pose a problem for downstream interpretation. Of course, a large variance can also be caused by genuinely batch-specific populations, so some prior knowledge about the biological context is necessary to distinguish between these two possibilities. For the PBMC dataset, none of the most variable clusters are overtly batch-specific, consistent with the fact that our batches are effectively replicates. ```r rv <- clusterAbundanceVar(tab.mnn) # Also printing the percentage of cells in each cluster in each batch: percent <- t(t(tab.mnn)/colSums(tab.mnn)) * 100 df <- DataFrame(Batch=unclass(percent), var=rv) df[order(df$var, decreasing=TRUE),] ``` ``` ## DataFrame with 16 rows and 3 columns ## Batch.1 Batch.2 var ## ## 15 0.153315 0.860832 0.934778 ## 13 5.519356 2.223816 0.728465 ## 9 0.421617 1.339072 0.707757 ## 8 6.209276 2.821616 0.563419 ## 4 0.459946 0.095648 0.452565 ## ... ... ... ... ## 6 0.651591 0.454328 0.05689945 ## 10 20.965887 25.896700 0.04527468 ## 2 11.077041 12.960306 0.02443988 ## 1 12.916826 14.490674 0.01318296 ## 16 0.153315 0.191296 0.00689661 ``` ## Preserving biological heterogeneity Another useful diagnostic check is to compare the pre-correction clustering of each batch to the clustering of the same cells in the corrected data. Accurate data integration should preserve population structure within each batch as there is no batch effect to remove between cells in the same batch. This check complements the previously mentioned diagnostics that only focus on the removal of differences between batches. Specifically, it protects us against scenarios where the correction method simply aggregates all cells together, which would achieve perfect mixing but also discard the biological heterogeneity of interest. To illustrate, we will use clustering results from the analysis of each batch of the PBMC dataset:
```r #--- loading ---# library(TENxPBMCData) all.sce <- list( pbmc3k=TENxPBMCData('pbmc3k'), pbmc4k=TENxPBMCData('pbmc4k'), pbmc8k=TENxPBMCData('pbmc8k') ) #--- quality-control ---# library(scater) stats <- high.mito <- list() for (n in names(all.sce)) { current <- all.sce[[n]] is.mito <- grep("MT", rowData(current)$Symbol_TENx) stats[[n]] <- perCellQCMetrics(current, subsets=list(Mito=is.mito)) high.mito[[n]] <- isOutlier(stats[[n]]$subsets_Mito_percent, type="higher") all.sce[[n]] <- current[,!high.mito[[n]]] } #--- normalization ---# all.sce <- lapply(all.sce, logNormCounts) #--- variance-modelling ---# library(scran) all.dec <- lapply(all.sce, modelGeneVar) all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1) #--- dimensionality-reduction ---# library(BiocSingular) set.seed(10000) all.sce <- mapply(FUN=runPCA, x=all.sce, subset_row=all.hvgs, MoreArgs=list(ncomponents=25, BSPARAM=RandomParam()), SIMPLIFY=FALSE) set.seed(100000) all.sce <- lapply(all.sce, runTSNE, dimred="PCA") set.seed(1000000) all.sce <- lapply(all.sce, runUMAP, dimred="PCA") #--- clustering ---# for (n in names(all.sce)) { g <- buildSNNGraph(all.sce[[n]], k=10, use.dimred='PCA') clust <- igraph::cluster_walktrap(g)$membership colLabels(all.sce[[n]]) <- factor(clust) } ```
```r pbmc3k <- all.sce$pbmc3k table(colLabels(pbmc3k)) ``` ``` ## ## 1 2 3 4 5 6 7 8 9 10 ## 487 154 603 514 31 150 179 333 147 11 ``` ```r pbmc4k <- all.sce$pbmc4k table(colLabels(pbmc4k)) ``` ``` ## ## 1 2 3 4 5 6 7 8 9 10 11 12 13 ## 497 185 569 786 373 232 44 1023 77 218 88 54 36 ``` Ideally, we should see a many-to-1 mapping where the post-correction clustering is nested inside the pre-correction clustering. This indicates that any within-batch structure was preserved after correction while acknowledging that greater resolution is possible with more cells. We quantify this mapping using the `nestedClusters()` function from the *[bluster](https://bioconductor.org/packages/3.18/bluster)* package, which identifies the nesting of post-correction clusters within the pre-correction clusters. Well-nested clusters have high `max` values, indicating that most of their cells are derived from a single pre-correction cluster. ```r library(bluster) tab3k <- nestedClusters(ref=paste("before", colLabels(pbmc3k)), alt=paste("after", clusters.mnn[mnn.out$batch==1])) tab3k$alt.mapping ``` ``` ## DataFrame with 16 rows and 2 columns ## max which ## ## after 1 0.985163 before 8 ## after 10 0.919561 before 3 ## after 11 1.000000 before 5 ## after 12 0.812500 before 5 ## after 13 0.993056 before 9 ## ... ... ... ## after 5 0.874275 before 4 ## after 6 0.882353 before 4 ## after 7 1.000000 before 1 ## after 8 0.981481 before 1 ## after 9 1.000000 before 1 ``` We can visualize this mapping for the PBMC dataset in Figure \@ref(fig:heat-after-mnn). Ideally, each row should have a single dominant entry close to unity. Horizontal stripes are more concerning as these indicate that multiple pre-correction clusters were merged together, though the exact level of concern will depend on whether specific clusters of interest are gained or lost. In practice, more discrepancies can be expected even when the correction is perfect, due to the existence of closely related clusters that were arbitrarily separated in the within-batch clustering. ```r library(pheatmap) # For the first batch: heat3k <- pheatmap(tab3k$proportions, cluster_row=FALSE, cluster_col=FALSE, main="PBMC 3K comparison", silent=TRUE) # For the second batch: tab4k <- nestedClusters(ref=paste("before", colLabels(pbmc4k)), alt=paste("after", clusters.mnn[mnn.out$batch==2])) heat4k <- pheatmap(tab4k$proportions, cluster_row=FALSE, cluster_col=FALSE, main="PBMC 4K comparison", silent=TRUE) gridExtra::grid.arrange(heat3k[[4]], heat4k[[4]]) ```
Comparison between the clusterings obtained before (columns) and after MNN correction (rows). One heatmap is generated for each of the PBMC 3K and 4K datasets, where each entry is colored according to the proportion of cells distributed along each row (i.e., the row sums equal unity).

(\#fig:heat-after-mnn)Comparison between the clusterings obtained before (columns) and after MNN correction (rows). One heatmap is generated for each of the PBMC 3K and 4K datasets, where each entry is colored according to the proportion of cells distributed along each row (i.e., the row sums equal unity).

We use the adjusted Rand index ([Advanced Section 5.3](http://bioconductor.org/books/3.18/OSCA.advanced/clustering-redux.html#comparing-different-clusterings)) to quantify the agreement between the clusterings before and after batch correction. Recall that larger indices are more desirable as this indicates that within-batch heterogeneity is preserved, though this must be balanced against the ability of each method to actually perform batch correction. ```r library(bluster) ri3k <- pairwiseRand(clusters.mnn[mnn.out$batch==1], colLabels(pbmc3k), mode="index") ri3k ``` ``` ## [1] 0.7361 ``` ```r ri4k <- pairwiseRand(clusters.mnn[mnn.out$batch==2], colLabels(pbmc4k), mode="index") ri4k ``` ``` ## [1] 0.8301 ``` We can also break down the ARI into per-cluster ratios for more detailed diagnostics (Figure \@ref(fig:rand-after-mnn)). For example, we could see low ratios off the diagonal if distinct clusters in the within-batch clustering were incorrectly aggregated in the merged clustering. Conversely, we might see low ratios on the diagonal if the correction inflated or introduced spurious heterogeneity inside a within-batch cluster. ```r # For the first batch. tab <- pairwiseRand(colLabels(pbmc3k), clusters.mnn[mnn.out$batch==1]) heat3k <- pheatmap(tab, cluster_row=FALSE, cluster_col=FALSE, col=rev(viridis::magma(100)), main="PBMC 3K probabilities", silent=TRUE) # For the second batch. tab <- pairwiseRand(colLabels(pbmc4k), clusters.mnn[mnn.out$batch==2]) heat4k <- pheatmap(tab, cluster_row=FALSE, cluster_col=FALSE, col=rev(viridis::magma(100)), main="PBMC 4K probabilities", silent=TRUE) gridExtra::grid.arrange(heat3k[[4]], heat4k[[4]]) ```
ARI-derived ratios for the within-batch clusters after comparison to the merged clusters obtained after MNN correction. One heatmap is generated for each of the PBMC 3K and 4K datasets.

(\#fig:rand-after-mnn)ARI-derived ratios for the within-batch clusters after comparison to the merged clusters obtained after MNN correction. One heatmap is generated for each of the PBMC 3K and 4K datasets.

## MNN-specific diagnostics For `fastMNN()`, one useful diagnostic is the proportion of variance within each batch that is lost during MNN correction. Specifically, this refers to the within-batch variance that is removed during orthogonalization with respect to the average correction vector at each merge step. This is returned via the `lost.var` field in the metadata of `mnn.out`, which contains a matrix of the variance lost in each batch (column) at each merge step (row). ```r metadata(mnn.out)$merge.info$lost.var ``` ``` ## [,1] [,2] ## [1,] 0.006617 0.003315 ``` Large proportions of lost variance (>10%) suggest that correction is removing genuine biological heterogeneity. This would occur due to violations of the assumption of orthogonality between the batch effect and the biological subspace [@haghverdi2018batch]. In this case, the proportion of lost variance is small, indicating that non-orthogonality is not a major concern. Another MNN-related diagnostic involves examining the variance in the differences in expression between MNN pairs. A small variance indicates that the correction had little effect - either there was no batch effect, or any batch effect was simply a constant shift across all cells. On the other hand, a large variance indicates that the correction was highly non-linear, most likely involving subpopulation-specific batch effects. This computation is achieved using the `mnnDeltaVariance()` function on the MNN pairings produced by `fastMNN()`. ```r library(batchelor) common <- rownames(mnn.out) vars <- mnnDeltaVariance(pbmc3k[common,], pbmc4k[common,], pairs=metadata(mnn.out)$merge.info$pairs) vars[order(vars$adjusted, decreasing=TRUE),] ``` ``` ## DataFrame with 13431 rows and 4 columns ## mean total trend adjusted ## ## ENSG00000111796 0.404086 1.487536 0.657367 0.830169 ## ENSG00000257764 0.320638 1.044040 0.548209 0.495831 ## ENSG00000170345 2.543318 1.342864 0.948955 0.393909 ## ENSG00000187118 0.313542 0.916383 0.538257 0.378126 ## ENSG00000177606 1.799190 1.400792 1.033322 0.367470 ## ... ... ... ... ... ## ENSG00000204482 0.786304 0.537895 0.97769 -0.439796 ## ENSG00000011600 1.176075 0.628421 1.07374 -0.445319 ## ENSG00000101439 1.135830 0.626148 1.07300 -0.446850 ## ENSG00000158869 0.927395 0.575360 1.03902 -0.463658 ## ENSG00000105374 1.219511 0.536416 1.07336 -0.536943 ``` Such genes with large variances are particularly interesting as they exhibit complex differences between batches that may reflect real biology. For example, in Figure \@ref(fig:mnn-delta-var-pbmc), the _KLRB1_-positive clusters in the second batch lack any counterpart in the first batch, despite the two batches being replicates. This may represent some kind of batch-specific state in two otherwise identical populations, though whether this is biological or technical in nature is open for interpretation. ```r library(scater) top <- rownames(vars)[order(vars$adjusted, decreasing=TRUE)[1]] gridExtra::grid.arrange( plotExpression(pbmc3k, x="label", features=top) + ggtitle("3k"), plotExpression(pbmc4k, x="label", features=top) + ggtitle("4k"), ncol=2 ) ```
Distribution of the expression of the gene with the largest variance of MNN pair differences in each batch of the the PBMC dataset.

(\#fig:mnn-delta-var-pbmc)Distribution of the expression of the gene with the largest variance of MNN pair differences in each batch of the the PBMC dataset.

## Session Info {-}
``` R version 4.3.2 Patched (2023-11-13 r85521) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 22.04.3 LTS Matrix products: default BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0 locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_GB LC_COLLATE=C [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C time zone: America/New_York tzcode source: system (glibc) attached base packages: [1] stats4 stats graphics grDevices utils datasets methods [8] base other attached packages: [1] scater_1.30.1 ggplot2_3.4.4 [3] scuttle_1.12.0 HDF5Array_1.30.1 [5] rhdf5_2.46.1 DelayedArray_0.28.0 [7] SparseArray_1.2.4 S4Arrays_1.2.0 [9] abind_1.4-5 Matrix_1.6-5 [11] pheatmap_1.0.12 bluster_1.12.0 [13] batchelor_1.18.1 BiocSingular_1.18.0 [15] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0 [17] Biobase_2.62.0 GenomicRanges_1.54.1 [19] GenomeInfoDb_1.38.6 IRanges_2.36.0 [21] S4Vectors_0.40.2 BiocGenerics_0.48.1 [23] MatrixGenerics_1.14.0 matrixStats_1.2.0 [25] BiocStyle_2.30.0 rebook_1.12.0 loaded via a namespace (and not attached): [1] bitops_1.0-7 gridExtra_2.3 [3] CodeDepends_0.6.5 rlang_1.1.3 [5] magrittr_2.0.3 compiler_4.3.2 [7] dir.expiry_1.10.0 DelayedMatrixStats_1.24.0 [9] vctrs_0.6.5 pkgconfig_2.0.3 [11] crayon_1.5.2 fastmap_1.1.1 [13] XVector_0.42.0 labeling_0.4.3 [15] utf8_1.2.4 rmarkdown_2.25 [17] ggbeeswarm_0.7.2 graph_1.80.0 [19] xfun_0.42 zlibbioc_1.48.0 [21] cachem_1.0.8 beachmat_2.18.1 [23] jsonlite_1.8.8 highr_0.10 [25] rhdf5filters_1.14.1 Rhdf5lib_1.24.2 [27] BiocParallel_1.36.0 irlba_2.3.5.1 [29] parallel_4.3.2 cluster_2.1.6 [31] R6_2.5.1 bslib_0.6.1 [33] RColorBrewer_1.1-3 limma_3.58.1 [35] jquerylib_0.1.4 Rcpp_1.0.12 [37] bookdown_0.37 knitr_1.45 [39] igraph_2.0.1.1 tidyselect_1.2.0 [41] yaml_2.3.8 viridis_0.6.5 [43] codetools_0.2-19 lattice_0.22-5 [45] tibble_3.2.1 withr_3.0.0 [47] evaluate_0.23 pillar_1.9.0 [49] BiocManager_1.30.22 filelock_1.0.3 [51] generics_0.1.3 RCurl_1.98-1.14 [53] sparseMatrixStats_1.14.0 munsell_0.5.0 [55] scales_1.3.0 glue_1.7.0 [57] metapod_1.10.1 tools_4.3.2 [59] BiocNeighbors_1.20.2 ScaledMatrix_1.10.0 [61] locfit_1.5-9.8 scran_1.30.2 [63] XML_3.99-0.16.1 cowplot_1.1.3 [65] grid_4.3.2 edgeR_4.0.15 [67] colorspace_2.1-0 GenomeInfoDbData_1.2.11 [69] beeswarm_0.4.0 vipor_0.4.7 [71] cli_3.6.2 rsvd_1.0.5 [73] fansi_1.0.6 viridisLite_0.4.2 [75] dplyr_1.1.4 ResidualMatrix_1.12.0 [77] gtable_0.3.4 sass_0.4.8 [79] digest_0.6.34 ggrepel_0.9.5 [81] dqrng_0.3.2 farver_2.1.1 [83] htmltools_0.5.7 lifecycle_1.0.4 [85] statmod_1.5.0 ```