seurat_03_integration.knit

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(Seurat)
    library(cowplot)
    library(ggplot2)
})

alldata <- readRDS("data/results/covid_qc_dr.rds")
print(names(alldata@reductions))

## [1] "pca"               "umap"              "tsne"             
## [4] "UMAP10_on_PCA"     "UMAP_on_ScaleData" "UMAP_on_Graph"

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”).

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

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap10_on_pca_ to umap10onpca_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_scaledata_ to umaponscaledata_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_graph_ to umapongraph_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap10_on_pca_ to umap10onpca_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_scaledata_ to umaponscaledata_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_graph_ to umapongraph_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap10_on_pca_ to umap10onpca_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_scaledata_ to umaponscaledata_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_graph_ to umapongraph_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap10_on_pca_ to umap10onpca_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_scaledata_ to umaponscaledata_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_graph_ to umapongraph_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap10_on_pca_ to umap10onpca_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_scaledata_ to umaponscaledata_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_graph_ to umapongraph_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap10_on_pca_ to umap10onpca_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_scaledata_ to umaponscaledata_

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from umap_on_graph_ to umapongraph_

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)
}

hvgs_per_dataset <- lapply(alldata.list, function(x) {
    x@assays$RNA@var.features
})
# venn::venn(hvgs_per_dataset,opacity = .4,zcolor = scales::hue_pal()(3),cexsn
# = 1,cexil = 1,lwd=1,col='white',frame=F,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"))

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

alldata.anchors <- FindIntegrationAnchors(object.list = alldata.list, dims = 1:30,
    reduction = "cca")

## Computing 2000 integration features

## Scaling features for provided objects

## Finding all pairwise anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2041 anchors

## Filtering anchors

##  Retained 1712 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2111 anchors

## Filtering anchors

##  Retained 1644 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2682 anchors

## Filtering anchors

##  Retained 2167 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 1814 anchors

## Filtering anchors

##  Retained 1425 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2297 anchors

## Filtering anchors

##  Retained 1863 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2505 anchors

## Filtering anchors

##  Retained 2025 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2151 anchors

## Filtering anchors

##  Retained 1612 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2653 anchors

## Filtering anchors

##  Retained 2226 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2916 anchors

## Filtering anchors

##  Retained 2161 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2934 anchors

## Filtering anchors

##  Retained 2485 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2057 anchors

## Filtering anchors

##  Retained 1621 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2546 anchors

## Filtering anchors

##  Retained 2186 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2738 anchors

## Filtering anchors

##  Retained 2092 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 2750 anchors

## Filtering anchors

##  Retained 2283 anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 3043 anchors

## Filtering anchors

##  Retained 2765 anchors

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

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

## Merging dataset 1 into 2

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

## Merging dataset 6 into 5

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

## Merging dataset 2 1 into 5 6

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

## Merging dataset 3 into 5 6 2 1

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

## Merging dataset 4 into 5 6 2 1 3

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

We can observe that a new assay slot is now created under the name CCA.

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 ingegrated data.

alldata.int@active.assay

## [1] "RNA" "CCA"
## [1] "CCA"

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.

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

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 08:47:25 UMAP embedding parameters a = 0.9922 b = 1.112

## 08:47:25 Read 5532 rows and found 30 numeric columns

## 08:47:25 Using Annoy for neighbor search, n_neighbors = 30

## 08:47:25 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 08:47:26 Writing NN index file to temp file /var/folders/f_/vj_w4xx933z1rr95yf4rhphr0000gp/T//RtmpS9VPsy/file17a876667960c
## 08:47:26 Searching Annoy index using 1 thread, search_k = 3000
## 08:47:27 Annoy recall = 100%
## 08:47:27 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 08:47:27 Initializing from normalized Laplacian + noise (using irlba)
## 08:47:28 Commencing optimization for 500 epochs, with 248064 positive edges
## 08:47:37 Optimization finished

alldata.int <- RunTSNE(alldata.int, dims = 1:30)

We can now plot the un-integrated and the integrated space reduced dimensions.

plot_grid(ncol = 3,
  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")
)