In this tutorial we will look at different ways of integrating multiple single cell RNA-seq datasets. We will explore two different methods to correct for batch effects across datasets. We will also look at a quantitative measure to assess the quality of the integrated data. Seurat uses the data integration method presented in Comprehensive Integration of Single Cell Data, while Scran and Scanpy use a mutual Nearest neighbour method (MNN). Below you can find a list of the most recent methods for single data integration:
| Markdown | Language | Library | Ref |
|---|---|---|---|
| CCA | R | Seurat | Cell |
| MNN | R/Python | Scater/Scanpy | Nat. Biotech. |
| Conos | R | conos | Nat. Methods |
| Scanorama | Python | scanorama | Nat. Biotech. |
Let’s first load necessary libraries and the data saved in the previous lab.
suppressPackageStartupMessages({
library(scater)
library(scran)
library(cowplot)
library(ggplot2)
library(rafalib)
# library(venn)
})
sce <- readRDS("data/results/covid_qc_dm.rds")
print(reducedDims(sce))## List of length 8
## names(8): PCA UMAP tSNE_on_PCA ... UMAP_on_ScaleData KNN UMAP_on_GraphWe split the combined object into a list, with each dataset as an element. We perform standard preprocessing (log-normalization), and identify variable features individually for each dataset based on a variance stabilizing transformation (“vst”).
sce.list <- lapply(unique(sce$sample), function(x) {
x <- sce[, sce$sample == x]
})
mypar(1, 3)
hvgs_per_dataset <- lapply(sce.list, function(x) {
x <- computeSumFactors(x, sizes = c(20, 40, 60, 80))
x <- logNormCounts(x)
var.out <- modelGeneVar(x, method = "loess")
hvg.out <- var.out[which(var.out$FDR <= 0.05 & var.out$bio >= 0.2), ]
hvg.out <- hvg.out[order(hvg.out$bio, decreasing = TRUE), ]
return(rownames(hvg.out))
})## Warning in (function (x, sizes, min.mean = NULL, positive = FALSE, scaling =
## NULL) : encountered non-positive size factor estimates
## Warning in (function (x, sizes, min.mean = NULL, positive = FALSE, scaling =
## NULL) : encountered non-positive size factor estimates
## Warning in (function (x, sizes, min.mean = NULL, positive = FALSE, scaling =
## NULL) : encountered non-positive size factor estimates
## Warning in (function (x, sizes, min.mean = NULL, positive = FALSE, scaling =
## NULL) : encountered non-positive size factor estimatesnames(hvgs_per_dataset) <- unique(sce$sample)
# venn::venn(hvgs_per_dataset,opacity = .4,zcolor = scales::hue_pal()(3),cexsn
# = 1,cexil = 1,lwd=1,col='white',borders = NA)
temp <- unique(unlist(hvgs_per_dataset))
overlap <- sapply(hvgs_per_dataset, function(x) {
temp %in% x
})
pheatmap::pheatmap(t(overlap * 1), cluster_rows = F, color = c("grey90", "grey20"))
The mutual nearest neighbors (MNN) approach within the scran package
utilizes a novel approach to adjust for batch effects. The
fastMNN() function returns a representation of the data
with reduced dimensionality, which can be used in a similar fashion to
other lower-dimensional representations such as PCA. In particular, this
representation can be used for downstream methods such as clustering.
The BNPARAM can be used to specify the specific nearest neighbors method
to use from the BiocNeighbors package. Here we make use of the Annoy library via the
BiocNeighbors::AnnoyParam() argument. We save the
reduced-dimension MNN representation into the reducedDims slot of our
sce object.
mnn_out <- batchelor::fastMNN(sce, subset.row = unique(unlist(hvgs_per_dataset)),
batch = factor(sce$sample), k = 20, d = 50)NOTE: fastMNN() does not produce a
batch-corrected expression matrix.
mnn_out <- t(reducedDim(mnn_out, "corrected"))
colnames(mnn_out) <- unlist(lapply(sce.list, function(x) {
colnames(x)
}))
mnn_out <- mnn_out[, colnames(sce)]
rownames(mnn_out) <- paste0("dim", 1:50)
reducedDim(sce, "MNN") <- t(mnn_out)We can observe that a new assay slot is now created under the name
MNN.
reducedDims(sce)## List of length 9
## names(9): PCA UMAP tSNE_on_PCA UMAP_on_PCA ... KNN UMAP_on_Graph MNNThus, the result from fastMNN() should solely be treated
as a reduced dimensionality representation, suitable for direct
plotting, TSNE/UMAP, clustering, and trajectory analysis that relies on
such results.
set.seed(42)
sce <- runTSNE(sce, dimred = "MNN", n_dimred = 50, perplexity = 30, name = "tSNE_on_MNN")
sce <- runUMAP(sce, dimred = "MNN", n_dimred = 50, ncomponents = 2, name = "UMAP_on_MNN")We can now plot the un-integrated and the integrated space reduced dimensions.
plot_grid(ncol = 3,
plotReducedDim(sce,dimred = "PCA",colour_by = "sample", point_size = 0.6)+ ggplot2::ggtitle(label ="PCA"),
plotReducedDim(sce,dimred = "tSNE_on_PCA",colour_by = "sample", point_size = 0.6)+ ggplot2::ggtitle(label ="tSNE_on_PCA"),
plotReducedDim(sce,dimred = "UMAP_on_PCA",colour_by = "sample",point_size = 0.6)+ ggplot2::ggtitle(label ="UMAP_on_PCA"),
plotReducedDim(sce,dimred = "MNN",colour_by = "sample", point_size = 0.6)+ ggplot2::ggtitle(label ="MNN"),
plotReducedDim(sce,dimred = "tSNE_on_MNN",colour_by = "sample", point_size = 0.6)+ ggplot2::ggtitle(label ="tSNE_on_MNN"),
plotReducedDim(sce,dimred = "UMAP_on_MNN",colour_by = "sample", point_size = 0.6)+ ggplot2::ggtitle(label ="UMAP_on_MNN")
)
Let’s plot some marker genes for different celltypes onto the embedding. Some genes are:
| Markers | Cell Type |
|---|---|
| CD3E | T cells |
| CD3E CD4 | CD4+ T cells |
| CD3E CD8A | CD8+ T cells |
| GNLY, NKG7 | NK cells |
| MS4A1 | B cells |
| CD14, LYZ, CST3, MS4A7 | CD14+ Monocytes |
| FCGR3A, LYZ, CST3, MS4A7 | FCGR3A+ Monocytes |
| FCER1A, CST3 | DCs |
plotlist <- list()
for (i in c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7",
"FCGR3A", "CST3", "FCER1A")) {
plotlist[[i]] <- plotReducedDim(sce, dimred = "UMAP_on_MNN", colour_by = i, by_exprs_values = "logcounts",
point_size = 0.6) + scale_fill_gradientn(colours = colorRampPalette(c("grey90",
"orange3", "firebrick", "firebrick", "red", "red"))(10)) + ggtitle(label = i) +
theme(plot.title = element_text(size = 20))
}
plot_grid(ncol = 3, plotlist = plotlist)
#INTEG_R1:
#INTEG_R2:
library(harmony)## Loading required package: RcppreducedDimNames(sce)
sce <- RunHarmony(sce, group.by.vars = "sample", reduction.save = "harmony", reduction = "PCA",
dims.use = 1:50)## Harmony 1/10## Harmony 2/10## Harmony 3/10## Harmony 4/10## Harmony 5/10## Harmony converged after 5 iterations# Here we use all PCs computed from Harmony for UMAP calculation
sce <- runUMAP(sce, dimred = "harmony", n_dimred = 50, ncomponents = 2, name = "UMAP_on_Harmony")## [1] "PCA" "UMAP" "tSNE_on_PCA"
## [4] "UMAP_on_PCA" "UMAP10_on_PCA" "UMAP_on_ScaleData"
## [7] "KNN" "UMAP_on_Graph" "MNN"
## [10] "tSNE_on_MNN" "UMAP_on_MNN"#INTEG_R3:
#INTEG_R4:
hvgs <- unique(unlist(hvgs_per_dataset))
scelist <- list()
genelist <- list()
for (i in 1:length(sce.list)) {
scelist[[i]] <- t(as.matrix(logcounts(sce.list[[i]])[hvgs, ]))
genelist[[i]] <- hvgs
}
lapply(scelist, dim)## [[1]]
## [1] 901 444
##
## [[2]]
## [1] 598 444
##
## [[3]]
## [1] 1052 444
##
## [[4]]
## [1] 1062 444
##
## [[5]]
## [1] 1175 444
##
## [[6]]
## [1] 1108 444#INTEG_R5:
library(reticulate)
scanorama <- import("scanorama")
integrated.data <- scanorama$integrate(datasets_full = scelist, genes_list = genelist)
intdimred <- do.call(rbind, integrated.data[[1]])
colnames(intdimred) <- paste0("PC_", 1:100)
rownames(intdimred) <- colnames(logcounts(sce))
# Add standard deviations in order to draw Elbow Plots in Seurat
stdevs <- apply(intdimred, MARGIN = 2, FUN = sd)
attr(intdimred, "varExplained") <- stdevs
reducedDim(sce, "Scanorama_PCA") <- intdimred
# Here we use all PCs computed from Scanorama for UMAP calculation
sce <- runUMAP(sce, dimred = "Scanorama_PCA", n_dimred = 50, ncomponents = 2, name = "UMAP_on_Scanorama")#INTEG_R6:
p1 <- plotReducedDim(sce, dimred = "UMAP_on_PCA", colour_by = "sample", point_size = 0.6) +
ggplot2::ggtitle(label = "UMAP_on_PCA")
p2 <- plotReducedDim(sce, dimred = "UMAP_on_MNN", colour_by = "sample", point_size = 0.6) +
ggplot2::ggtitle(label = "UMAP_on_MNN")
p3 <- plotReducedDim(sce, dimred = "UMAP_on_Harmony", colour_by = "sample", point_size = 0.6) +
ggplot2::ggtitle(label = "UMAP_on_Harmony")
p4 <- plotReducedDim(sce, dimred = "UMAP_on_Scanorama", colour_by = "sample", point_size = 0.6) +
ggplot2::ggtitle(label = "UMAP_on_Scanorama")
leg <- get_legend(p1)
gridExtra::grid.arrange(gridExtra::arrangeGrob(p1 + Seurat::NoLegend() + Seurat::NoAxes(),
p2 + Seurat::NoLegend() + Seurat::NoAxes(), p3 + Seurat::NoLegend() + Seurat::NoAxes(),
p4 + Seurat::NoLegend() + Seurat::NoAxes(), nrow = 2), leg, ncol = 2, widths = c(8,
2))
#INTEG_R7:
Finally, lets save the integrated data for further analysis.
saveRDS(sce, "data/results/covid_qc_dr_int.rds")Session Info
sessionInfo()## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS/LAPACK: /Users/asabjor/miniconda3/envs/scRNAseq2023/lib/libopenblasp-r0.3.21.dylib
##
## locale:
## [1] C/UTF-8/C/C/C/C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] reticulate_1.27 harmony_1.0
## [3] Rcpp_1.0.10 rafalib_1.0.0
## [5] cowplot_1.1.1 scran_1.22.1
## [7] scater_1.22.0 ggplot2_3.4.0
## [9] scuttle_1.4.0 SingleCellExperiment_1.16.0
## [11] SummarizedExperiment_1.24.0 Biobase_2.54.0
## [13] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1
## [15] IRanges_2.28.0 S4Vectors_0.32.4
## [17] BiocGenerics_0.40.0 MatrixGenerics_1.6.0
## [19] matrixStats_0.63.0 RJSONIO_1.3-1.7
## [21] optparse_1.7.3
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 spatstat.explore_3.0-5
## [3] tidyselect_1.2.0 htmlwidgets_1.6.1
## [5] grid_4.1.3 BiocParallel_1.28.3
## [7] Rtsne_0.16 munsell_0.5.0
## [9] ScaledMatrix_1.2.0 codetools_0.2-18
## [11] ica_1.0-3 statmod_1.5.0
## [13] future_1.30.0 miniUI_0.1.1.1
## [15] withr_2.5.0 batchelor_1.10.0
## [17] spatstat.random_3.0-1 colorspace_2.1-0
## [19] progressr_0.13.0 highr_0.10
## [21] knitr_1.41 Seurat_4.3.0
## [23] ROCR_1.0-11 tensor_1.5
## [25] listenv_0.9.0 labeling_0.4.2
## [27] GenomeInfoDbData_1.2.7 polyclip_1.10-4
## [29] farver_2.1.1 pheatmap_1.0.12
## [31] rprojroot_2.0.3 parallelly_1.34.0
## [33] vctrs_0.5.2 generics_0.1.3
## [35] xfun_0.36 R6_2.5.1
## [37] ggbeeswarm_0.7.1 rsvd_1.0.5
## [39] locfit_1.5-9.7 bitops_1.0-7
## [41] spatstat.utils_3.0-1 cachem_1.0.6
## [43] DelayedArray_0.20.0 assertthat_0.2.1
## [45] promises_1.2.0.1 scales_1.2.1
## [47] beeswarm_0.4.0 gtable_0.3.1
## [49] beachmat_2.10.0 globals_0.16.2
## [51] goftest_1.2-3 rlang_1.0.6
## [53] splines_4.1.3 lazyeval_0.2.2
## [55] spatstat.geom_3.0-5 yaml_2.3.7
## [57] reshape2_1.4.4 abind_1.4-5
## [59] httpuv_1.6.8 tools_4.1.3
## [61] ellipsis_0.3.2 jquerylib_0.1.4
## [63] RColorBrewer_1.1-3 ggridges_0.5.4
## [65] plyr_1.8.8 sparseMatrixStats_1.6.0
## [67] zlibbioc_1.40.0 purrr_1.0.1
## [69] RCurl_1.98-1.9 deldir_1.0-6
## [71] pbapply_1.7-0 viridis_0.6.2
## [73] zoo_1.8-11 SeuratObject_4.1.3
## [75] ggrepel_0.9.2 cluster_2.1.4
## [77] here_1.0.1 magrittr_2.0.3
## [79] data.table_1.14.6 scattermore_0.8
## [81] ResidualMatrix_1.4.0 lmtest_0.9-40
## [83] RANN_2.6.1 fitdistrplus_1.1-8
## [85] patchwork_1.1.2 mime_0.12
## [87] evaluate_0.20 xtable_1.8-4
## [89] gridExtra_2.3 compiler_4.1.3
## [91] tibble_3.1.8 KernSmooth_2.23-20
## [93] htmltools_0.5.4 later_1.3.0
## [95] tidyr_1.2.1 DBI_1.1.3
## [97] formatR_1.14 MASS_7.3-58.2
## [99] Matrix_1.5-3 getopt_1.20.3
## [101] cli_3.6.0 parallel_4.1.3
## [103] metapod_1.2.0 igraph_1.3.5
## [105] pkgconfig_2.0.3 sp_1.6-0
## [107] plotly_4.10.1 spatstat.sparse_3.0-0
## [109] vipor_0.4.5 bslib_0.4.2
## [111] dqrng_0.3.0 XVector_0.34.0
## [113] stringr_1.5.0 digest_0.6.31
## [115] sctransform_0.3.5 RcppAnnoy_0.0.20
## [117] spatstat.data_3.0-0 rmarkdown_2.20
## [119] leiden_0.4.3 uwot_0.1.14
## [121] edgeR_3.36.0 DelayedMatrixStats_1.16.0
## [123] shiny_1.7.4 nlme_3.1-161
## [125] lifecycle_1.0.3 jsonlite_1.8.4
## [127] BiocNeighbors_1.12.0 viridisLite_0.4.1
## [129] limma_3.50.3 fansi_1.0.4
## [131] pillar_1.8.1 lattice_0.20-45
## [133] fastmap_1.1.0 httr_1.4.4
## [135] survival_3.5-0 glue_1.6.2
## [137] png_0.1-8 bluster_1.4.0
## [139] stringi_1.7.12 sass_0.4.5
## [141] BiocSingular_1.10.0 dplyr_1.0.10
## [143] irlba_2.3.5.1 future.apply_1.10.0