Comparing integration tools on heterogeneous T cell datasets

Introduction

Here we will compare integration results of STACAS with other R packages on a collection of scRNA-seq datasets of mouse tumor-infiltrating T cells from multiple studies from Andreatta et al. (2021) Nat Comms. It consists of seven datasets from six different studies covering tumor and lymph node samples, comprising studies with both CD4+ and CD8+ T cells (MC38_dLN, Ekiz and Xiong), only CD8+ T cells (Carmona, Singer) or only CD4+ T cells (Magen_dLN and Magen_TILs). The collection of datasets is therefore imbalanced in terms of cell types, posing a challenging task to integration tools.

The data are available on figshare at: figshare/12478571

R environment

Get and load some useful packages

renv::restore()

if (!require("remotes", quietly = TRUE))
    install.packages("remotes")
library(remotes)

if (!requireNamespace("STACAS", quietly = TRUE))
  remotes::install_github("carmonalab/STACAS")

if (!requireNamespace("harmony", quietly = TRUE))
  install.packages("Harmony")

if (!requireNamespace("SeuratWrappers", quietly = TRUE))
remotes::install_github('satijalab/seurat-wrappers')

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

if (!require("batchelor", quietly = TRUE))
BiocManager::install("batchelor")

if (!require("scIntegrationMetrics", quietly = TRUE))
install_github("carmonalab/scIntegrationMetrics") #calculates LISI and Silhouette
library(Seurat)
library(dplyr)
library(ggplot2)
library(STACAS)
library(scGate)
library(harmony)
library(SeuratWrappers)
library(batchelor)
library(tidyr)
library(scIntegrationMetrics)
library(patchwork)
library(tidytext)

seed = 1234
set.seed(seed)

Download and load datasets

The commands below will download the collection of datasets as a Seurat object.

download <- T
where <- 'aux'
dir.create(where, showWarnings = FALSE)

rds.path <- sprintf("%s/ref_TILAtlas_mouse_v1.rds", where)

if(download){
  options(timeout=500)
  url <- "https://figshare.com/ndownloader/files/41398167"
  download.file(url = url, destfile = rds.path)
}

object <- readRDS(rds.path)
DefaultAssay(object) <- "RNA"

Annotated subtypes and dataset/study are stored in functional.cluster and Study metadata column, respectively.

meta.batch <- "Study"
meta.label <- "functional.cluster"

For integration metrics, we will not consider labels that are contributed by a single dataset/batch for more than thrMaxPerBatch (90%).

batchLabels.table <- table(object@meta.data[[meta.label]], object@meta.data[[meta.batch]])

batchLabels <- round( batchLabels.table / rowSums(batchLabels.table) * 100) 
batchLabels
##                     
##                      Carmona Ekiz Magen_dLN Magen_TILs MC38_dLN Singer Xiong
##   CD8_Tex                 18   10         0          0        1      8    63
##   CD8_Tpex                27    3         0          4        3      7    55
##   CD8_EffectorMemory      70    3         0          0        3      4    20
##   CD8_EarlyActiv          59    4         0          0       24      4     9
##   CD8_NaiveLike           63    1         0          0       31      1     3
##   CD4_NaiveLike            1    1         3          0       91      0     5
##   Tfh                      0    0        78          2       20      0     0
##   Th1                      3    4         4         29       37      0    23
##   Treg                     1    2         2         11       12      0    72
thrMaxPerBatch <- 90
removeLabels <- names(which(apply(batchLabels, 1, function(x) max(x) > thrMaxPerBatch)))

removeLabels

[1] “CD4_NaiveLike”

metricsLabels <- setdiff(unique(object@meta.data[[meta.label]]),removeLabels)

Set critical parameters

nfeatures <- 2000  
ndim <- 20

Standardize gene symbols, to ensure that the same genes are used across datasets.

do_standard_symbols <- TRUE

if (do_standard_symbols) {
  object@assays$RNA@counts <- object@assays$RNA@data
  data("EnsemblGeneTable.Mm")
  
  object <- StandardizeGeneSymbols(object, EnsemblGeneTable=EnsemblGeneTable.Mm)
}
## Number of genes in input object: 27231
## Number of genes with standard symbols: 24772 (90.97%)
## Examples of non-standard Gene.names:
## Hist1h4k,Fam49a,Gm9595,Hist1h2ak,Gm18192,Hist1h2bn
## Additional number of genes with accepted Gene name synonym:  1589
## Additional number of genes after replacing dots:  0
## Number of duplicated Gene.name: 939 (3.45%)
## Final number of genes: 25422 (93.36%)
## Renaming default assay from RNA to tmp

Uncorrected data

As a baseline, we can evaluate integration metrics on the uncorrected data prior to any batch effect correction.

First, run a standard Seurat pipeline for dimensionality reduction:

object <- object |> NormalizeData() |>
  FindVariableFeatures(nfeatures = nfeatures) |>
  ScaleData() |>
  RunPCA(npcs=ndim) |>
  RunUMAP(dims=1:ndim)
p1 <- DimPlot(object, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (uncorrected)")
p2 <- DimPlot(object, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (uncorrected)")

p1 | p2

Clusters are largely driven by study/batch of origin rather than subtype. We can quantify batch mixing and preservation of biological variability by applying some of the metrics implemented in the scIntegrationMetrics package.

To quantify batch mixing, CiLISI is a cell type-aware version of the integration LISI that accounts for different composition between datasets. On the other hand, preservation of biological variability can be quantified by how close to each other cells of the same type are, and how separated from each other cells of different types are in the joint integrated embedding. A useful metric for preservation of biological variability is Average Silhouette Width (ASW) of cell labels (celltype_ASW) in the corrected PCA space.

Integration metrics prior to integration

integrationMetrics <- list()
useMetrics <- c("CiLISI","celltype_ASW")

method = "uncorrected"

integrationMetrics[[method]] <- getIntegrationMetrics(object=object,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$uncorrected
## $CiLISI
## [1] 0.009686983
## 
## $celltype_ASW
## [1] 0.04816328

Consistently variable features

A simple approach to mitigate batch effects is to select HVGs that are consistently variable across datasets. For this we will split by dataset/batch, calculate HVG for each, then identify shared genes using SelectIntegrationFeatures from Seurat.

obj.list <- SplitObject(object, split.by = meta.batch)
for (i in 1:length(obj.list)) {
  obj.list[[i]] <- obj.list[[i]] |>
                  NormalizeData(assay="RNA") |>
                  FindVariableFeatures(nfeatures=nfeatures*2)
}

hvg <- SelectIntegrationFeatures(obj.list, nfeatures = nfeatures)

Re-calculate dimensionality reduction using hvg

object.sh <- object
object.sh@assays$RNA@var.features <- hvg
object.sh <- ScaleData(object.sh) |> RunPCA(npcs=ndim) |> RunUMAP(dims=1:ndim)
p1 <- DimPlot(object.sh, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (shared features)")
p2 <- DimPlot(object.sh, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (shared features)")

p1 | p2

method = "sharedFeatures"

integrationMetrics[[method]] <- getIntegrationMetrics(object=object.sh,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$sharedFeatures
## $CiLISI
## [1] 0.01246234
## 
## $celltype_ASW
## [1] 0.1071285

STACAS integration

Run STACAS with default parameters on the list of objects

object.stacas <- Run.STACAS(obj.list, dims = 1:ndim, anchor.features = hvg) |>
      RunUMAP(dims = 1:ndim) 
p1 <- DimPlot(object.stacas, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (STACAS)")
p2 <- DimPlot(object.stacas, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (STACAS)")

p1 | p2

integrationMetrics[["STACAS"]] <- getIntegrationMetrics(object=object.stacas,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$STACAS
## $CiLISI
## [1] 0.1282444
## 
## $celltype_ASW
## [1] 0.106899

Semi-supervised STACAS integration

When available, cell type annotations can be used to guide the alignment. STACAS will use this information to penalize anchors where cell types are inconsistent.

For this example, we will build a simple scGate model for CD4 vs. CD8 T cells. The predicted labels will be used to guide STACAS integration.

library(scGate)
my.genes.blocklist <- scGate::genes.blacklist.default$Mm

CD8T <- scGate::gating_model(name="CD8T", signature=c("Cd8a","Cd8b1"))
CD4T <- scGate::gating_model(name="CD4T", signature=c("Cd4","Cd40lg"))
models <- list("CD8T"=CD8T, "CD4T"=CD4T)
models
## $CD8T
##   levels   use_as name  signature
## 1 level1 positive CD8T Cd8a;Cd8b1
## 
## $CD4T
##   levels   use_as name  signature
## 1 level1 positive CD4T Cd4;Cd40lg

We can run these scGate models on individual datasets, to predict CD4 and CD8 T cells.

obj.list <- lapply(obj.list, function(x) {
  scGate(x, model=models, multi.asNA=TRUE)
})

We can see that some datasets contain only CD4 or CD8 T cells, others contain both:

table(obj.list[["Xiong"]]$scGate_multi)
## 
## CD4T CD8T 
##  539 3272
table(obj.list[["Ekiz"]]$scGate_multi)
## 
## CD8T 
##  488
table(obj.list[["Magen_TILs"]]$scGate_multi)
## 
## CD4T 
##  497

We can now apply ssSTACAS with the predicted cell type labels, which will be used to guide the integration:

object.ss <- Run.STACAS(obj.list, dims = 1:ndim,
                        anchor.features = hvg,
                        cell.labels = "scGate_multi") |>
              RunUMAP(dims=1:ndim)

Note that there is no need for ALL cells to be annotated: we recommend to set labels to NA or unknown for cells that cannot be confidently annotated, and they won’t be penalized for label inconsistency. In addition, you can decide how much weight to give to cell labels with the label.confidence parameter (from 0 to 1).

p1 <- DimPlot(object.ss, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (ssSTACAS CD4/CD8)")
p2 <- DimPlot(object.ss, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (ssSTACAS CD4/CD8)")

p1 | p2

integrationMetrics[["ssSTACAS"]] <- getIntegrationMetrics(object=object.ss,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$ssSTACAS
## $CiLISI
## [1] 0.1281298
## 
## $celltype_ASW
## [1] 0.1239857

Semi-supervised STACAS with full annotations

In general, cell type labels obtained by any manual or automatic annotation can be fed as input to ssSTACAS, by specifying the metadata annotation column with the cell.labels parameter. In a best-case scenario where all cell labels are known:

object.ss.full <- Run.STACAS(obj.list, dims = 1:ndim,
                             anchor.features = hvg,
                             cell.labels = meta.label) |>
                  RunUMAP(dims=1:ndim)

Note that there is no need for ALL cells to be annotated: we recommend to set labels to NA or unknown for cells that cannot be confidently annotated, and they won’t be penalized for label inconsistency. In addition, you can decide how much weight to give to cell labels with the label.confidence parameter (from 0 to 1).

p1 <- DimPlot(object.ss.full, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (ssSTACAS)")
p2 <- DimPlot(object.ss.full, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (ssSTACAS)")

p1 | p2

integrationMetrics[["ssSTACAS_full"]] <- getIntegrationMetrics(object=object.ss.full,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$ssSTACAS_full
## $CiLISI
## [1] 0.1231655
## 
## $celltype_ASW
## [1] 0.1430031

Seurat CCA

Canonical correlation analysis (CCA) is the default integration method implemented in Seurat. CCA aims to find linear combinations of features across data sets that are maximally correlated, identifying shared correlation structures across data sets (see Butler et al.). Note: CCA can be very slow!

# Find anchors
cca.anchors <- FindIntegrationAnchors(obj.list, anchor.features = hvg,
                                      reduction = "cca", dims = 1:ndim)

# Integrate data
object.cca <- IntegrateData(cca.anchors, dims=1:ndim)
rm(cca.anchors)


# Visualize on UMAP space
object.cca <- object.cca |> ScaleData() |>
  RunPCA(npcs=ndim) |> RunUMAP(dims=1:ndim)
p1 <- DimPlot(object.cca, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (Seurat CCA)")
p2 <- DimPlot(object.cca, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (Seurat CCA")

p1 | p2

integrationMetrics[["Seurat_CCA"]] <- getIntegrationMetrics(object=object.cca,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$Seurat_CCA
## $CiLISI
## [1] 0.1791048
## 
## $celltype_ASW
## [1] -0.01923617

By visual inspection and by the integration metrics, we can see the CCA achieve a high batch mixing, but performs poorly in terms of bio-conservation (low celltype_ASW). This combination is typical of methods that overcorrect batch effects, i.e. they mix well datasets, by they also mix cells of different types.

Seurat rPCA

Reciprocal PCA (rPCA) is presented in Seurat as an alternative to CCA that is faster and less prone to overcorrection.

# Find anchors
rpca.anchors <- FindIntegrationAnchors(obj.list, anchor.features = hvg,
                                       reduction = "rpca", dims = 1:ndim)

# Integrate data
object.rpca <- IntegrateData(anchorset = rpca.anchors, dims=1:ndim)

rm(rpca.anchors)

object.rpca <- object.rpca |> ScaleData() |>
  RunPCA(npcs=ndim) |> RunUMAP(dims=1:ndim)
p1 <- DimPlot(object.rpca, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (Seurat rPCA)")
p2 <- DimPlot(object.rpca, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (Seurat rPCA)")

p1 | p2

integrationMetrics[["Seurat_rPCA"]] <- getIntegrationMetrics(object=object.rpca,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$Seurat_rPCA
## $CiLISI
## [1] 0.1597844
## 
## $celltype_ASW
## [1] 0.04414

Harmony

Harmony is a very popular method for single-cell integration data, based on iterative clustering of cells in reduced dimensionality spaces (see Korsunsky et al.). Because it works directly in low dimensional spaces, this method is one of the fastest tools for batch effect correction.

object.harmony <- RunHarmony(object, group.by.vars =  meta.batch) 

object.harmony <- RunUMAP(object.harmony, reduction = "harmony", dims=1:ndim)
p1 <- DimPlot(object.harmony, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (Harmony)")
p2 <- DimPlot(object.harmony, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (Harmony)")

p1 | p2

integrationMetrics[["Harmony"]] <- getIntegrationMetrics(object=object.harmony,
                                                      metrics = useMetrics,
                                                      method.reduction = "harmony",
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$Harmony
## $CiLISI
## [1] 0.05410919
## 
## $celltype_ASW
## [1] -0.005395199

FastMNN

FastMNN is implemented in the batchelor Bioconductor package. It calculates mutual nearest neighbors (MNNs) in low-dimensional representations of the input data, and uses such MNNs to infer the magnitude and direction of the batch effect.

object.fastMNN <- RunFastMNN(obj.list, features = hvg, assay="RNA", d=ndim)

object.fastMNN <- RunUMAP(object.fastMNN, reduction = "mnn", dims=1:ndim)
p1 <- DimPlot(object.fastMNN, group.by = meta.batch) +
  theme(aspect.ratio = 1) + ggtitle("Dataset/batch (FastMNN)")
p2 <- DimPlot(object.fastMNN, group.by = meta.label,
              label=T, label.size = 5) + NoLegend() +
  theme(aspect.ratio = 1) + ggtitle("Cell labels (FastMNN)")

p1 | p2

integrationMetrics[["FastMNN"]] <- getIntegrationMetrics(object=object.fastMNN,
                                                      metrics = useMetrics,
                                                      method.reduction = "mnn",
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$FastMNN
## $CiLISI
## [1] 0.1229723
## 
## $celltype_ASW
## [1] 0.1083248

Summary of Integration Metrics

library(ggrepel)

integrationMetricsSummary <- data.frame(unlist(integrationMetrics)) |>
  tibble::rownames_to_column() |>
  dplyr::rename(value=unlist.integrationMetrics.) |>
  separate(rowname, c("Method","Metric"), sep="\\.")


integrationMetricsSummary |> filter(Metric %in% c("celltype_ASW","CiLISI")) |>
  pivot_wider(names_from = Metric, values_from = value ) |>
  ggplot(aes(x=CiLISI, y=celltype_ASW, label=Method)) +
  geom_point(aes(color=Method)) +
  geom_label_repel(aes(label = Method),
                     box.padding   = 0.35, 
                     point.padding = 0.5,
                     segment.color = 'grey50') + theme_light()

Semi-supervised STACAS is the best performing method when accounting for both batch mixing (high CiLISI) and preservation of biological variability (high celltype_ASW). Even limited and incomplete information about CD4 vs CD8 T cells can be beneficial for integration, as seen by increased performance of ssSTACAS compared to unsupervised STACAS. FastMNN also performed well, with similar performance to unsupervised STACAS. Popular methods such as Harmony, Seurat rPCA and Seurat CCA perform poorly on this integration task characterized by high cell type imbalance. In particular, Seurat CCA appears to overcorrect batch effects, indiscriminately mixing cells of different types. Importantly, a good selection of variable features to reduce the dimensionality of the data can alone improve the preservation of biological variability (see sharedFeatures vs. uncorrected).

Ideally, integration methods should be able mix well different datasets (high CiLISI) while keeping together cells of the same type (high celltype_ASW). For the integration of balanced datasets in terms of cell types (e.g. technical replicates of the same experiment), methods such as Seurat CCA or Harmony can perform well. However, in integration tasks characterized by high cell type imbalance, it is preferable to choose integration methods that preserve biological variability, such as STACAS and FastMNN. Additionally, when prior cell type information is available, it can be used as input to supervised methods such ssSTACAS to guide integration and avoid overcorrection of batch effects.

Further reading

The STACAS package and installation instructions are available at: STACAS package

The code for this demo can be found on GitHub

References

  • Andreatta A., Carmona S. J. (2021). STACAS: Sub-Type Anchor Correction for Alignment in Seurat to integrate single-cell RNA-seq data. - Bioinformatics

  • Andreatta M, Herault L, Gueguen P, Gfeller D, Berenstein AJ, Carmona SJ (2024) Semi-supervised integration of single-cell transcriptomics data - Nature Communications

  • Korsunsky, Ilya, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nature methods 16.12 (2019): 1289-1296.

  • Butler, Andrew, et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nature biotechnology 36.5 (2018): 411-420.

  • Hao, Y., Hao, S., Andersen-Nissen, E., Mauck III, W. M., Zheng, S., Butler, A., … & Satija, R. (2021). Integrated analysis of multimodal single-cell data. - Cell