In satijalab/seurat: Tools for Single Cell Genomics

all_times <- list()  # store the time for each chunk
knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      now <<- Sys.time()
    } else {
      res <- difftime(Sys.time(), now, units = "secs")
      all_times[[options$label]] <<- res
    }
  }
}))
knitr::opts_chunk$set(
  tidy = TRUE,
  tidy.opts = list(width.cutoff = 95),
  message = FALSE,
  warning = FALSE,
  fig.width = 10,
  time_it = TRUE,
  error = TRUE
)

The recent increase in publicly available single-cell datasets poses a significant challenge for integrative analysis. For example, multiple tissues have now been profiled across dozens of studies, representing hundreds of individuals and millions of cells. In Hao et al, 2023 proposed a dictionary learning based method, atomic sketch integration, could also enable efficient and large-scale integrative analysis. Our procedure enables the integration of large compendiums of datasets without ever needing to load the full scale of data into memory. In our manuscript we use atomic sketch integration to integrate millions of scRNA-seq from human lung and human PBMC.

In this vignette, we demonstrate how to use atomic sketch integration to harmonize scRNA-seq experiments 1M cells, though we have used this procedure to integrate datasets of 10M+ cells as well. We analyze a dataset from Parse Biosciences, in which PBMC from 24 human samples (12 healthy donors, 12 Type-1 diabetes donors), which is available here.

Sample a representative subset of cells ('atoms') from each dataset
Integrate the atoms from each dataset, and define a set of cell states
Reconstruct (integrate) the full datasets, based on the atoms
Annotate all cells in the full datasets
Identify cell-type specific differences between healthy and diabetic patients

Prior to running this vignette, please install Seurat v5, as well as the BPCells package, which we use for on-disk storage. You can read more about using BPCells in Seurat v5 here. We also recommend reading the Sketch-based analysis in Seurat v5 vignette, which introduces the concept of on-disk and in-memory storage in Seurat v5.

library(Seurat)
library(BPCells)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(patchwork)
# set this option when analyzing large datasets
options(future.globals.maxSize = 3e9)

Create a Seurat object containing data from 24 patients

We downloaded the original dataset and donor metadata from Parse Biosciences. While the BPCells package can work directly with h5ad files, for optimal performance, we converted the dataset to the compressed sparse format used by BPCells, as described here. We create a Seurat object for this dataset. Since the input to CreateSeuratObject is a BPCells matrix, the data remains on-disk and is not loaded into memory. After creating the object, we split the dataset into 24 layers, one for each sample (i.e. patient), to facilitate integration.

parse.mat <- open_matrix_dir(dir = "/brahms/hartmana/vignette_data/bpcells/parse_1m_pbmc")
# need to move
metadata <- readRDS("/brahms/haoy/vignette_data/ParseBio_PBMC_meta.rds")
object <- CreateSeuratObject(counts = parse.mat, meta.data = metadata)

object <- NormalizeData(object)
# split assay into 24 layers
object[['RNA']] <- split(object[['RNA']], f = object$sample)
object <- FindVariableFeatures(object, verbose = FALSE)

Sample representative cells from each dataset

Inspired by pioneering work aiming to identify 'sketches' of scRNA-seq data, our first step is to sample a representative set of cells from each dataset. We compute a leverage score (estimate of 'statistical leverage') for each cell, which helps to identify cells that are likely to be member of rare subpopulations and ensure that these are included in our representative sample. Importantly, the estimation of leverage scores only requires data normalization, can be computed efficiently for sparse datasets, and does not require any intensive computation or dimensional reduction steps. We load each object separately, perform basic preprocessing (normalization and variable feature selection), and select and store 5,000 representative cells from each dataset. Since there are 24 datasets, the sketched dataset now contains 120,000 cells. These cells are stored in a new sketch assay, and are loaded in-memory.

object <- SketchData(object = object, ncells = 5000, method = 'LeverageScore', sketched.assay = 'sketch')
object

Perform integration on the sketched cells across samples

Next we perform integrative analysis on the 'atoms' from each of the datasets. Here, we perform integration using the streamlined Seurat v5 integration worfklow, and utilize the reference-based RPCAIntegration method. The function performs all corrections in low-dimensional space (rather than on the expression values themselves) to further improve speed and memory usage, and outputs a merged Seurat object where all cells have been placed in an integrated low-dimensional space (stored as integrated.rpca). However, we emphasize that you can perform integration here using any analysis technique that places cells across datasets into a shared space. This includes CCA Integration, Harmony, and scVI. We demonstrate how to use these tools in Seurat v5 here.

DefaultAssay(object) <- 'sketch'
object <- FindVariableFeatures(object, verbose = F)
object <- ScaleData(object, verbose = F)
object <- RunPCA(object, verbose = F)
# integrate the datasets
object <- IntegrateLayers(object, method = RPCAIntegration, orig = 'pca',
                            new.reduction = 'integrated.rpca', dims = 1:30, k.anchor = 20,
                            reference = which(Layers(object, search = 'data') %in% c( 'data.H_3060')), 
                            verbose = F)
# cluster the integrated data
object <- FindNeighbors(object, reduction = 'integrated.rpca', dims = 1:30)
object <- FindClusters(object, resolution = 2)
object <- RunUMAP(object,  reduction = 'integrated.rpca', dims = 1:30, return.model = T, verbose = F)

# you can now rejoin the layers in the sketched assay
# this is required to perform differential expression
object[['sketch']] <- JoinLayers(object[['sketch']])
c10_markers <- FindMarkers(object = object, ident.1 = 10, max.cells.per.ident = 500, only.pos = TRUE)
head(c10_markers)

# You can now annotate clusters using marker genes. 
# We performed this step, and include the results in the 'sketch.celltype' metadata column

plot.s1 <- DimPlot(object, group.by = 'sample', reduction = 'umap')
plot.s2 <- DimPlot(object, group.by = 'celltype.manual', reduction = 'umap')

plot.s1 + plot.s2 + plot_layout(ncol = 1)

Integrate the full datasets

Now that we have integrated the subset of atoms of each dataset, placing them each in an integrated low-dimensional space, we can now place each cell from each dataset in this space as well. We load the full datasets back in individually, and use the ProjectIntegration function to integrate all cells. After this function is run, the integrated.rpca.full space now embeds all cells in the dataset.Even though all cells in the dataset have been integrated together, the non-sketched cells are not loaded into memory. Users can still switch between the sketch (sketched cells, in-memory) and RNA (full dataset, on disk) for analysis. After integration, we can also project cell type labels from the sketched cells onto the full dataset using ProjectData.

# resplit the sketched cell assay into layers
# this is required to project the integration onto all cells
object[['sketch']] <- split(object[['sketch']], f = object$sample)

object <- ProjectIntegration(object = object,
                                    sketched.assay =  'sketch',
                                    assay = 'RNA',
                                    reduction = 'integrated.rpca'
                             )


object <- ProjectData(object = object,
                      sketched.assay =  'sketch',
                      assay = 'RNA',
                      sketched.reduction = 'integrated.rpca.full',
                      full.reduction = 'integrated.rpca.full',
                      dims = 1:30,
                      refdata = list(celltype.full = 'celltype.manual')
                      )

object <- RunUMAP(object,  reduction = 'integrated.rpca.full', dims = 1:30 , reduction.name = 'umap.full', reduction.key = 'UMAP_full_')

p1 <- DimPlot(object, reduction = 'umap.full', group.by = 'sample',alpha = 0.1)
p2 <- DimPlot(object, reduction = 'umap.full', group.by = 'celltype.full', alpha = 0.1)
p1 + p2 + plot_layout(ncol = 1)

Compare healthy and diabetic samples

By integrating all samples together, we can now compare healthy and diabetic cells in matched cell states. To maximize statistical power, we want to use all cells - not just the sketched cells - to perform this analysis. As recommended by Soneson et all. and Crowell et al., we use an aggregation-based (pseudobulk) workflow. We aggregate all cells within the same cell type and sample using the AggregateExpression function. This returns a Seurat object where each 'cell' represents the pseudobulk profile of one cell type in one individual.

After we aggregate cells, we can perform celltype-specific differential expression between healthy and diabetic samples using DESeq2. We demonstrate this for CD14 monocytes.

bulk <- AggregateExpression(object, return.seurat = T, slot = 'counts',
            assays  = 'RNA', group.by = c("celltype.full","sample", 'disease'))

# each sample is an individual-specific celltype-specific pseudobulk profile
tail(Cells(bulk))

cd14.bulk <- subset(bulk,celltype.full == "CD14 Mono")
Idents(cd14.bulk) <- 'disease'
de_markers <- FindMarkers(cd14.bulk, ident.1 = 'D',ident.2 = 'H', slot = 'counts', test.use = 'DESeq2', verbose = F )
de_markers$gene <- rownames(de_markers)    
ggplot(de_markers, aes(avg_log2FC, -log10(p_val))) + geom_point(size=0.5, alpha=0.5)   + theme_bw() + ylab("-log10(unadjusted p-value)")+geom_text_repel(aes(label = ifelse(p_val_adj<0.01, gene, "")),colour = 'red', size = 3)