---
description: Combining and harmonizing samples or datasets from
  different batches such as experiments or conditions to enable
  meaningful cross-sample comparisons.
subtitle:  Seurat Toolkit
title:  Data Integration
---

<div>

> **Note**
>
> Code chunks run R commands unless otherwise specified.

</div>

In this tutorial we will look at different ways of integrating multiple
single cell RNA-seq datasets. We will explore a few different methods to
correct for batch effects across datasets. 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 some methods for single data
integration:

  -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Markdown          Language          Library           Ref
  ----------------- ----------------- ----------------- -----------------------------------------------------------------------------------------------------------------------------------
  CCA               R                 Seurat            [Cell](https://www.sciencedirect.com/science/article/pii/S0092867419305598?via%3Dihub)

  MNN               R/Python          Scater/Scanpy     [Nat. Biotech.](https://www.nature.com/articles/nbt.4091)

  Conos             R                 conos             [Nat.
                                                        Methods](https://www.nature.com/articles/s41592-019-0466-z?error=cookies_not_supported&code=5680289b-6edb-40ad-9934-415dac4fdb2f)

  Conos             R                 Harmony           [Nat. Methods](https://www.nature.com/articles/s41592-019-0619-0)

  Scanorama         Python            scanorama         [Nat. Biotech.](https://www.nature.com/articles/s41587-019-0113-3)
  -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

## Data preparation

Let's first load necessary libraries and the data saved in the previous
lab.

``` {r}
suppressPackageStartupMessages({
    library(Seurat)
    library(ggplot2)
    library(patchwork)
    library(reticulate)
})

# Activate scanorama Python venv
reticulate::use_virtualenv("/opt/venv/scanorama")
reticulate::py_discover_config()
```

``` {r}
# download pre-computed data if missing or long compute
fetch_data <- TRUE

# url for source and intermediate data
path_data <- "https://export.uppmax.uu.se/naiss2023-23-3/workshops/workshop-scrnaseq"
path_file <- "data/covid/results/seurat_covid_qc_dr.rds"
if (!dir.exists(dirname(path_file))) dir.create(dirname(path_file), recursive = TRUE)
if (fetch_data && !file.exists(path_file)) download.file(url = file.path(path_data, "covid/results/seurat_covid_qc_dr.rds"), destfile = path_file)
alldata <- readRDS(path_file)
print(names(alldata@reductions))
```

We 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**).

``` {r}
#| fig-height: 6
#| fig-width: 8

alldata.list <- SplitObject(alldata, split.by = "orig.ident")

for (i in 1:length(alldata.list)) {
    alldata.list[[i]] <- NormalizeData(alldata.list[[i]], verbose = FALSE)
    alldata.list[[i]] <- FindVariableFeatures(alldata.list[[i]], selection.method = "vst", nfeatures = 2000,verbose = FALSE)
}

# get the variable genes from all the datasets.
hvgs_per_dataset <- lapply(alldata.list, function(x) { x@assays$RNA@var.features })

# also add in the variable genes that was selected on the whole dataset
hvgs_per_dataset$all = VariableFeatures(alldata)

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"))
```

As you can see, there are a lot of genes that are variable in just one
dataset. There are also some genes in the gene set that was selected
using all the data that are not variable in any of the individual
datasets. These are most likely genes driven by batch effects.

A better way to select features for integration is to combine the
information on variable genes across the dataset. This can be done with
the function `SelectIntegrationFeatures` that combines the ranks of the
variable features in the different datasets.

``` {r}
hvgs_all = SelectIntegrationFeatures(alldata.list)
hvgs_per_dataset$all_ranks = hvgs_all

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"))
```

For all downstream integration we will use this set of genes.

## CCA

We identify anchors using the `FindIntegrationAnchors()` function, which
takes a list of Seurat objects as input.

``` {r}
alldata.anchors <- FindIntegrationAnchors(object.list = alldata.list, dims = 1:30,reduction = "cca", anchor.features = hvgs_all)
```

We then pass these anchors to the `IntegrateData()` function, which
returns a Seurat object.

``` {r}
alldata.int <- IntegrateData(anchorset = alldata.anchors, dims = 1:30, new.assay.name = "CCA")
```

We can observe that a new assay slot is now created under the name
`CCA`. If you do not specify the assay name, the default will be
`integrated`.

``` {r}
names(alldata.int@assays)

# by default, Seurat now sets the integrated assay as the default assay, so any operation you now perform will be on the integrated data.
alldata.int@active.assay
```

After running `IntegrateData()`, the Seurat object will contain a new
Assay with the integrated (or **batch-corrected**) expression matrix.
Note that the original (uncorrected values) are still stored in the
object in the "RNA" assay, so you can switch back and forth. We can then
use this new integrated matrix for downstream analysis and
visualization. Here we scale the integrated data, run PCA, and visualize
the results with UMAP and TSNE. The integrated datasets cluster by cell
type, instead of by technology.

As `CCA` is the `active.assay` now, the functions will by default run on
the data in that assay. But you could also specify in each of the
functions to run them in a specific assay with the parameter
`assay = "CCA"`.

``` {r}
#Run Dimensionality reduction on integrated space
alldata.int <- ScaleData(alldata.int, verbose = FALSE)
alldata.int <- RunPCA(alldata.int, npcs = 30, verbose = FALSE)
alldata.int <- RunUMAP(alldata.int, dims = 1:30)
alldata.int <- RunTSNE(alldata.int, dims = 1:30)
```

We can now plot the unintegrated and the integrated space reduced
dimensions.

``` {r}
#| fig-height: 6
#| fig-width: 10

wrap_plots(
  DimPlot(alldata, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA raw_data"),
  DimPlot(alldata, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE raw_data"),
  DimPlot(alldata, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP raw_data"),
  
  DimPlot(alldata.int, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"),
  DimPlot(alldata.int, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"),
  DimPlot(alldata.int, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated"),
  ncol = 3
) + plot_layout(guides = "collect")
```

### Clean memory

Again we have a lot of large objects in the memory. We have the original
data `alldata` but also the integrated data in `alldata.int`. We also
have the split objects in `alldata.list` and the anchors in
`alldata.anchors`. In principle we only need the integrated object for
now, but we will also keep the list for running Scanorama further down
in the tutorial.

We also want to keep the original umap for visualization purposes, so we
copy it over to the integrated object.

``` {r}
alldata.int@reductions$umap_raw = alldata@reductions$umap

# remove all objects that will not be used.
rm(alldata,  alldata.anchors)
# run garbage collect to free up memory
gc()
```

Let's plot some marker genes for different cell types onto the
embedding.

  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

``` {r}
#| fig-height: 8
#| fig-width: 10

myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7", "FCGR3A", "CST3", "FCER1A")
FeaturePlot(alldata.int, reduction = "umap", dims = 1:2, features = myfeatures, ncol = 4, order = T) + NoLegend() + NoAxes() + NoGrid()
```

## Harmony

An alternative method for integration is Harmony, for more details on
the method, please se their paper [Nat.
Methods](https://www.nature.com/articles/s41592-019-0619-0). This method
runs the integration on a dimensionality reduction, in most applications
the PCA. So first, we will rerun scaling and PCA with the same set of
genes that were used for the CCA integration.

OBS! Make sure to revert back to the `RNA` assay.

``` {r}
alldata.int@active.assay = "RNA"
VariableFeatures(alldata.int) = hvgs_all
alldata.int = ScaleData(alldata.int, vars.to.regress = c("percent_mito", "nFeature_RNA"))
alldata.int = RunPCA(alldata.int, reduction.name = "pca_harmony")
```

Now we are ready to run Harmony.

``` {r}
#| fig-height: 10
#| fig-width: 13

library(harmony)

alldata.int <- RunHarmony(
  alldata.int,
  group.by.vars = "orig.ident",
  reduction.use = "pca_harmony",
  dims.use = 1:50,
  assay.use = "RNA")
```

Harmony will create another reduction slot in your seurat object with
the name **harmony**, so now we can use that reduction instead of PCA to
run UMAP.

``` {r}
alldata.int <- RunUMAP(alldata.int, dims = 1:50, reduction = "harmony", reduction.name = "umap_harmony")

DimPlot(alldata.int, reduction = "umap_harmony", group.by = "orig.ident") + NoAxes() + ggtitle("Harmony UMAP")
```

## Scanorama

<div>

> **Important**
>
> If you are running locally using Docker and you have a Mac with ARM
> chip, the Scanorama reticulate module is known to crash. In this case,
> you might want to skip this section.

</div>

Another integration method is Scanorama (see [Nat.
Biotech.](https://www.nature.com/articles/s41587-019-0113-3)). This
method is implemented in python, but we can run it through the
Reticulate package.

We will run it with the same set of variable genes, but first we have to
create a list of all the objects per sample.

``` {r}
#| fig-height: 5
#| fig-width: 16

assaylist <- list()
genelist <- list()
for(i in 1:length(alldata.list)) {
  assaylist[[i]] <- t(as.matrix(GetAssayData(alldata.list[[i]], "data")[hvgs_all,]))
  genelist[[i]] <- hvgs_all
}

lapply(assaylist,dim)
```

Then, we use the `scanorama` function through reticulate. The integrated
data is added back into the Seurat object as a new Reduction.

``` {r}
#| fig-height: 5
#| fig-width: 16

# Activate scanorama Python venv
scanorama <- reticulate::import("scanorama")

integrated.data <- scanorama$integrate(datasets_full = assaylist,
                                       genes_list = genelist )

# Now we create a new dim reduction object in the format that Seurat uses
intdimred <- do.call(rbind, integrated.data[[1]])
colnames(intdimred) <- paste0("PC_", 1:100)
rownames(intdimred) <- colnames(alldata.int)

# Add standard deviations in order to draw Elbow Plots in Seurat
stdevs <- apply(intdimred, MARGIN = 2, FUN = sd)

# Create a new dim red object.
alldata.int[["scanorama"]] <- CreateDimReducObject(
  embeddings = intdimred,
  stdev      = stdevs,
  key        = "PC_",
  assay      = "RNA")
```

``` {r}
#Here we use all PCs computed from Scanorama for UMAP calculation
alldata.int <- RunUMAP(alldata.int, dims = 1:100, reduction = "scanorama",reduction.name = "umap_scanorama")

DimPlot(alldata.int, reduction = "umap_scanorama", group.by = "orig.ident") + NoAxes() + ggtitle("Harmony UMAP")
```

## Overview all methods

Now we will plot UMAPS with all three integration methods side by side.

``` {r}
#| fig-height: 8
#| fig-width: 9

p1 <- DimPlot(alldata.int, reduction = "umap_raw", group.by = "orig.ident") + ggtitle("UMAP raw_data")
p2 <- DimPlot(alldata.int, reduction = "umap", group.by = "orig.ident") + ggtitle("UMAP CCA")
p3 <- DimPlot(alldata.int, reduction = "umap_harmony", group.by = "orig.ident") + ggtitle("UMAP Harmony")
p4 <- DimPlot(alldata.int, reduction = "umap_scanorama", group.by = "orig.ident")+ggtitle("UMAP Scanorama")

wrap_plots(p1, p2, p3, p4, nrow = 2) + plot_layout(guides = "collect")
```

<div>

> **Discuss**
>
> Look at the different integration results, which one do you think
> looks the best? How would you motivate selecting one method over the
> other? How do you think you could best evaluate if the integration
> worked well?

</div>

Let's save the integrated data for further analysis.

``` {r}
saveRDS(alldata.int,"data/covid/results/seurat_covid_qc_dr_int.rds")
```

## Extra task

You have now done the Seurat integration with CCA which is quite slow.
There are other options in the `FindIntegrationAnchors()` function. Try
rerunning the integration with `rpca` and/or `rlsi` and create a new
UMAP. Compare the results.

## Session info

```{=html}
<details>
```
```{=html}
<summary>
```
Click here
```{=html}
</summary>
```
``` {r}
sessionInfo()
```

```{=html}
</details>
```