STACAS.compareMetrics.Tcells.Rmd

---
title: "Comparing integration tools on heterogeneous T cell datasets"
author: "M. Andreatta and S. Carmona"
date: "29/01/2024"
output:
  rmdformats::readthedown:
    self-contained: true
    highlight: haddock
    thumbnails: false
    css: styles.css
knit: (function(input_file, encoding) {
  out_dir <- 'docs';
  rmarkdown::render(input_file, encoding=encoding,
 output_file=file.path(dirname(input_file), out_dir, 'Tcell.demo.html'))})
---

# 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](https://www.nature.com/articles/s41467-021-23324-4). 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](https://figshare.com/account/projects/123709/articles/12478571)

# R environment

Get and load some useful packages
```{r message=F, warning=F,results=F, eval=T}
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

```


```{r message=F, warning=F,results=F}
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.
```{r}
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.

```{r}
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%).

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

batchLabels <- round( batchLabels.table / rowSums(batchLabels.table) * 100) 
batchLabels
thrMaxPerBatch <- 90
```

```{r results='asis'}
removeLabels <- names(which(apply(batchLabels, 1, function(x) max(x) > thrMaxPerBatch)))

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



Set critical parameters
```{r}
nfeatures <- 2000  
ndim <- 20 
```


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

```{r, eval=T, results='asis', collapse=T}
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)
}
```

# 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:
```{r message=F, warning=F,results=F}
object <- object |> NormalizeData() |>
  FindVariableFeatures(nfeatures = nfeatures) |>
  ScaleData() |>
  RunPCA(npcs=ndim) |>
  RunUMAP(dims=1:ndim)
```
```{r fig.height=6, fig.width=14}
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](https://github.com/carmonalab/scIntegrationMetrics).

To quantify batch mixing, CiLISI is a cell type-aware version of the [integration LISI](https://github.com/immunogenomics/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
```{r, warning=F, message=F, collapse=T }
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
```


# 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.

```{r}
obj.list <- SplitObject(object, split.by = meta.batch)
```

```{r message=F, warning=F,results=F, eval=T}

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`
```{r message=F, warning=F,results=F}
object.sh <- object
object.sh@assays$RNA@var.features <- hvg
object.sh <- ScaleData(object.sh) |> RunPCA(npcs=ndim) |> RunUMAP(dims=1:ndim)
```

```{r fig.height=6, fig.width=14}
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
```

```{r message=F, warning=F,collapse=T} 
method = "sharedFeatures"

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

integrationMetrics$sharedFeatures

```


# STACAS integration

Run STACAS with default parameters on the list of objects
```{r, message=F, warning=F,results=F}
object.stacas <- Run.STACAS(obj.list, dims = 1:ndim, anchor.features = hvg) |>
      RunUMAP(dims = 1:ndim) 
```

```{r fig.height=6, fig.width=14}
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
```

```{r message=F, warning=F,collapse=T} 
integrationMetrics[["STACAS"]] <- getIntegrationMetrics(object=object.stacas,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$STACAS

```

# 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](https://github.com/carmonalab/scGate) model for CD4 vs. CD8 T cells. The predicted labels will be used to guide STACAS integration.

```{r warning=F, message=F}
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
```

We can run these scGate models on individual datasets, to predict CD4 and CD8 T cells.
```{r results=F, warning=F, message=F}
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:
```{r collapse=T}
table(obj.list[["Xiong"]]$scGate_multi)
table(obj.list[["Ekiz"]]$scGate_multi)
table(obj.list[["Magen_TILs"]]$scGate_multi)
```

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

```{r, results=F, warning=F, message=F }
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).


```{r fig.height=6, fig.width=14}
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
```


```{r message=F, warning=F,collapse=T} 
integrationMetrics[["ssSTACAS"]] <- getIntegrationMetrics(object=object.ss,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$ssSTACAS

```


# 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:

```{r, results=F, warning=F, message=F }
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).

```{r fig.height=6, fig.width=14}
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
```

```{r message=F, warning=F, collapse=T} 
integrationMetrics[["ssSTACAS_full"]] <- getIntegrationMetrics(object=object.ss.full,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$ssSTACAS_full

```

# 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.](https://www.nature.com/articles/nbt.4096)). Note: CCA can be very slow!

```{r, results=F, warning=F, message=F,  eval=T}
# 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)
```

```{r fig.height=6, fig.width=14}
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
```

```{r message=F, warning=F, collapse=T} 
integrationMetrics[["Seurat_CCA"]] <- getIntegrationMetrics(object=object.cca,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$Seurat_CCA

```

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](https://satijalab.org/seurat/articles/integration_rpca.html). 

```{r, results=F, warning=F, message=F }
# 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)

```


```{r fig.height=6, fig.width=14}
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
```

```{r message=F, warning=F, collapse=T} 
integrationMetrics[["Seurat_rPCA"]] <- getIntegrationMetrics(object=object.rpca,
                                                      metrics = useMetrics,
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$Seurat_rPCA
```


# 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.](https://www.nature.com/articles/s41592-019-0619-0)). Because it works directly in low dimensional spaces, this method is one of the fastest tools for batch effect correction.

```{r, message=F, warning=F,results=F}
object.harmony <- RunHarmony(object, group.by.vars =  meta.batch) 

object.harmony <- RunUMAP(object.harmony, reduction = "harmony", dims=1:ndim)
```

```{r fig.height=6, fig.width=14}
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
```

```{r message=F, warning=F, collapse=T} 
integrationMetrics[["Harmony"]] <- getIntegrationMetrics(object=object.harmony,
                                                      metrics = useMetrics,
                                                      method.reduction = "harmony",
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$Harmony
```
# FastMNN

FastMNN is implemented in the [batchelor](https://bioconductor.org/packages/release/bioc/html/batchelor.html) 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.

```{r, message=F, warning=F,results=F}
object.fastMNN <- RunFastMNN(obj.list, features = hvg, assay="RNA", d=ndim)

object.fastMNN <- RunUMAP(object.fastMNN, reduction = "mnn", dims=1:ndim)
```

```{r fig.height=6, fig.width=14}
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
```

```{r message=F, warning=F, collapse=T} 
integrationMetrics[["FastMNN"]] <- getIntegrationMetrics(object=object.fastMNN,
                                                      metrics = useMetrics,
                                                      method.reduction = "mnn",
                                                      meta.label = meta.label,
                                                      meta.batch = meta.batch,
                                                      metricsLabels = metricsLabels)

integrationMetrics$FastMNN
```


# Summary of Integration Metrics

```{r fig.height=4, fig.width=6}
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](https://github.com/carmonalab/STACAS)

The code for this demo can be found on [GitHub](https://github.com/carmonalab/STACAS.demo)

# 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