Combining and harmonizing samples or datasets from different batches such as experiments or conditions to enable meaningful cross-sample comparisons.
Authors
Åsa Björklund
Paulo Czarnewski
Susanne Reinsbach
Roy Francis
Published
09-Feb-2024
Note
Code chunks run R commands unless otherwise specified.
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:
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")for (i in1: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 datasethvgs_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.
hvgs_all =SelectIntegrationFeatures(alldata.list)hvgs_per_dataset$all_ranks = hvgs_alltemp <-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.
2 CCA
We identify anchors using the FindIntegrationAnchors() function, which takes a list of Seurat objects as input.
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.
names(alldata.int@assays)
[1] "RNA" "CCA"
# 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
[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.
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".
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.
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 memorygc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 3414522 182.4 4989418 266.5 4989418 266.5
Vcells 203287569 1551.0 566105528 4319.1 882313251 6731.6
Let’s plot some marker genes for different cell types onto the embedding.
An alternative method for integration is Harmony, for more details on the method, please se their paper Nat. Methods. 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.
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.
Another integration method is Scanorama (see Nat. Biotech.). 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.
Then, we use the scanorama function through reticulate. The integrated data is added back into the Seurat object as a new Reduction.
# Activate scanorama Python venvscanorama <- 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 usesintdimred <-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 Seuratstdevs <-apply(intdimred, MARGIN =2, FUN = sd)# Create a new dim red object.alldata.int[["scanorama"]] <-CreateDimReducObject(embeddings = intdimred,stdev = stdevs,key ="PC_",assay ="RNA")
#Here we use all PCs computed from Scanorama for UMAP calculationalldata.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")
5 Overview all methods
Now we will plot UMAPS with all three integration methods side by side.
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?
Let’s save the integrated data for further analysis.
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.
