Fast integration using reciprocal PCA (RPCA)
Compiled: January 11, 2022
Source:vignettes/integration_rpca.Rmd
integration_rpca.RmdIn this vignette, we present a slightly modified workflow for the integration of scRNA-seq datasets. Instead of utilizing canonical correlation analysis (‘CCA’) to identify anchors, we instead utilize reciprocal PCA (‘RPCA’). When determining anchors between any two datasets using RPCA, we project each dataset into the others PCA space and constrain the anchors by the same mutual neighborhood requirement. The commands for both workflows are largely identical, but the two methods may be applied in different context.
By identifying shared sources of variation between datasets, CCA is well-suited for identifying anchors when cell types are conserved, but there are very substantial differences in gene expression across experiments. CCA-based integration therefore enables integrative analysis when experimental conditions or disease states introduce very strong expression shifts, or when integrating datasets across modalities and species. However, CCA-based integration may also lead to overcorrection, especially when a large proportion of cells are non-overlapping across datasets.
RPCA-based integration runs significantly faster, and also represents a more conservative approach where cells in different biological states are less likely to ‘align’ after integration. We therefore,recommend RPCA during integrative analysis where: * A substantial fraction of cells in one dataset have no matching type in the other * Datasets originate from the same platform (i.e. multiple lanes of 10x genomics) * There are a large number of datasets or cells to integrate (see INSERT LINK for more tips on integrating large datasets)
Below, we demonstrate the use of reciprocal PCA to align the same stimulated and resting datasets first analyzed in our introduction to scRNA-seq integration vignette. While the list of commands is nearly identical, this workflow requires users to run principal components analysis (PCA) individually on each dataset prior to integration. Users should also set the ‘reduction’ argument to ‘rpca’, when running FindIntegrationAnchors().
library(SeuratData)
# install dataset
InstallData("ifnb")
# load dataset
LoadData("ifnb")
# split the dataset into a list of two seurat objects (stim and CTRL)
ifnb.list <- SplitObject(ifnb, split.by = "stim")
# normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
# select features that are repeatedly variable across datasets for integration run PCA on each
# dataset using these features
features <- SelectIntegrationFeatures(object.list = ifnb.list)
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
x <- ScaleData(x, features = features, verbose = FALSE)
x <- RunPCA(x, features = features, verbose = FALSE)
})Perform integration
We then identify anchors using the FindIntegrationAnchors() function, which takes a list of Seurat objects as input, and use these anchors to integrate the two datasets together with IntegrateData().
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features, reduction = "rpca")
# this command creates an 'integrated' data assay
immune.combined <- IntegrateData(anchorset = immune.anchors)Now we can run a single integrated analysis on all cells!
# specify that we will perform downstream analysis on the corrected data note that the
# original unmodified data still resides in the 'RNA' assay
DefaultAssay(immune.combined) <- "integrated"
# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)
# Visualization
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", group.by = "seurat_annotations", label = TRUE,
repel = TRUE)
p1 + p2
Modifying the strength of integration
The results show that rpca-based integration is more conservative, and in this case, do not perfectly align a subset of cells (which are naive and memory T cells) across experiments. You can increase the strength of alignment by increasing the k.anchor parameter, which is set to 5 by default. Increasing this parameter to 20 will assist in aligning these populations.
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features, reduction = "rpca",
k.anchor = 20)
immune.combined <- IntegrateData(anchorset = immune.anchors)
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 13999
## Number of edges: 589767
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9094
## Number of communities: 15
## Elapsed time: 4 seconds
# Visualization
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2
Now that the datasets have been integrated, you can follow the previous steps in the introduction to scRNA-seq integration vignette to identify cell types and cell type-specific responses.
Performing integration on datasets normalized with SCTransform
As an additional example, we repeat the analyses performed above, but normalize the datasets using SCTransform. We may choose to set the method parameter to glmGamPoi (install here) in order to enable faster estimation of regression parameters in SCTransform().
LoadData("ifnb")
ifnb.list <- SplitObject(ifnb, split.by = "stim")
ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform, method = "glmGamPoi")
features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)
ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)
ifnb.list <- lapply(X = ifnb.list, FUN = RunPCA, features = features)
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",
anchor.features = features, dims = 1:30, reduction = "rpca", k.anchor = 20)
immune.combined.sct <- IntegrateData(anchorset = immune.anchors, normalization.method = "SCT", dims = 1:30)
immune.combined.sct <- RunPCA(immune.combined.sct, verbose = FALSE)
immune.combined.sct <- RunUMAP(immune.combined.sct, reduction = "pca", dims = 1:30)
# Visualization
p1 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "seurat_annotations", label = TRUE,
repel = TRUE)
p1 + p2
Session Info
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [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
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.3.5 SeuratObject_4.0.4
## [3] Seurat_4.0.6 thp1.eccite.SeuratData_3.1.5
## [5] stxBrain.SeuratData_0.1.1 ssHippo.SeuratData_3.1.4
## [7] pbmcsca.SeuratData_3.0.0 pbmcMultiome.SeuratData_0.1.1
## [9] pbmc3k.SeuratData_3.1.4 panc8.SeuratData_3.0.2
## [11] ifnb.SeuratData_3.1.0 hcabm40k.SeuratData_3.0.0
## [13] bmcite.SeuratData_0.3.0 SeuratData_0.2.1
##
## loaded via a namespace (and not attached):
## [1] systemfonts_1.0.2 plyr_1.8.6
## [3] igraph_1.2.11 lazyeval_0.2.2
## [5] splines_4.1.0 listenv_0.8.0
## [7] scattermore_0.7 GenomeInfoDb_1.28.1
## [9] digest_0.6.29 htmltools_0.5.2
## [11] fansi_1.0.0 magrittr_2.0.1
## [13] memoise_2.0.0 tensor_1.5
## [15] cluster_2.1.2 ROCR_1.0-11
## [17] globals_0.14.0 matrixStats_0.61.0
## [19] pkgdown_1.6.1 spatstat.sparse_2.1-0
## [21] colorspace_2.0-2 rappdirs_0.3.3
## [23] ggrepel_0.9.1 textshaping_0.3.5
## [25] xfun_0.25 dplyr_1.0.7
## [27] RCurl_1.98-1.3 crayon_1.4.2
## [29] jsonlite_1.7.2 spatstat.data_2.1-2
## [31] survival_3.2-11 zoo_1.8-9
## [33] glue_1.6.0 polyclip_1.10-0
## [35] gtable_0.3.0 zlibbioc_1.38.0
## [37] XVector_0.32.0 leiden_0.3.9
## [39] DelayedArray_0.18.0 future.apply_1.8.1
## [41] BiocGenerics_0.38.0 abind_1.4-5
## [43] scales_1.1.1 DBI_1.1.1
## [45] miniUI_0.1.1.1 Rcpp_1.0.7
## [47] viridisLite_0.4.0 xtable_1.8-4
## [49] reticulate_1.22 spatstat.core_2.3-2
## [51] stats4_4.1.0 htmlwidgets_1.5.4
## [53] httr_1.4.2 RColorBrewer_1.1-2
## [55] ellipsis_0.3.2 ica_1.0-2
## [57] pkgconfig_2.0.3 farver_2.1.0
## [59] sass_0.4.0 uwot_0.1.11
## [61] deldir_1.0-6 utf8_1.2.2
## [63] tidyselect_1.1.1 labeling_0.4.2
## [65] rlang_0.4.12 reshape2_1.4.4
## [67] later_1.3.0 munsell_0.5.0
## [69] tools_4.1.0 cachem_1.0.6
## [71] cli_3.1.0 generics_0.1.1
## [73] ggridges_0.5.3 evaluate_0.14
## [75] stringr_1.4.0 fastmap_1.1.0
## [77] yaml_2.2.1 ragg_1.1.3
## [79] goftest_1.2-3 knitr_1.33
## [81] fs_1.5.2 fitdistrplus_1.1-6
## [83] purrr_0.3.4 RANN_2.6.1
## [85] sparseMatrixStats_1.4.2 pbapply_1.5-0
## [87] future_1.23.0 nlme_3.1-152
## [89] mime_0.12 formatR_1.11
## [91] compiler_4.1.0 plotly_4.10.0
## [93] png_0.1-7 spatstat.utils_2.3-0
## [95] tibble_3.1.6 bslib_0.3.1
## [97] glmGamPoi_1.4.0 stringi_1.7.6
## [99] highr_0.9 desc_1.4.0
## [101] RSpectra_0.16-0 lattice_0.20-44
## [103] Matrix_1.3-3 vctrs_0.3.8
## [105] pillar_1.6.4 lifecycle_1.0.1
## [107] spatstat.geom_2.3-1 lmtest_0.9-39
## [109] jquerylib_0.1.4 RcppAnnoy_0.0.19
## [111] bitops_1.0-7 data.table_1.14.2
## [113] cowplot_1.1.1 irlba_2.3.5
## [115] GenomicRanges_1.44.0 httpuv_1.6.5
## [117] patchwork_1.1.1 R6_2.5.1
## [119] promises_1.2.0.1 KernSmooth_2.23-20
## [121] gridExtra_2.3 IRanges_2.26.0
## [123] parallelly_1.30.0 codetools_0.2-18
## [125] MASS_7.3-54 assertthat_0.2.1
## [127] SummarizedExperiment_1.22.0 rprojroot_2.0.2
## [129] withr_2.4.3 sctransform_0.3.2
## [131] GenomeInfoDbData_1.2.6 S4Vectors_0.30.0
## [133] mgcv_1.8-35 parallel_4.1.0
## [135] grid_4.1.0 rpart_4.1-15
## [137] tidyr_1.1.4 DelayedMatrixStats_1.14.2
## [139] rmarkdown_2.10 MatrixGenerics_1.4.2
## [141] Rtsne_0.15 Biobase_2.52.0
## [143] shiny_1.7.1