Notebook

Label cell types using Seurat Label Transfer¶

To build our reference, we would like to start with labels that originate from published cell type references.

One of the approaches for this cell type labeling is Seurat, which labels cells by integration with a reference dataset.

Label transfer using Seurat is described on their website, and was introduced in this publication:

Stuart, T. et al. Comprehensive Integration of Single-Cell Data. Cell 177, 1888–1902.e21 (2019)

Here, we'll load in our cells in batches, and assign cell types based on the PBMC reference dataset provided by the Satija lab as part of their 2021 publication in Cell (described below).

Load libraries¶

dplyr: Data frame manipulation tools
H5weaver: A package for reading .h5 files generated by AIFI
hise: The HISE SDK
parallel: Parallelization of processes in R
purrr: Functional programing tools
Seurat: Single-cell data analysis tools
SeuratObject: Data structures for Seurat

We also set the timeout option to be high so that R waits for us to download the large reference dataset from Zenodo, below

In [1]:

quiet_library <- function(...) { suppressPackageStartupMessages(library(...)) }

quiet_library(dplyr)
quiet_library(H5weaver)
quiet_library(hise)
quiet_library(parallel)
quiet_library(purrr)
quiet_library(Seurat)
quiet_library(SeuratObject)

options(timeout = 10000)

Prepare Seurat PBMC reference¶

To use Seurat to label our cells, we'll utilize the PBMC reference provided by the Satija Lab for PBMCs generated using Seurat V5. This reference is derived from data in this publication from the Satija lab:

Hao, Y. and Hao, S. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587.e29 (2021)

The version of record for this reference dataset is provided in a Zenodo repository at this accession:
https://zenodo.org/records/7779017

Additional information about the cell type labels in this reference is available on the Azimuth website.

We'll download the reference from Zenodo for label transfer:

In [2]:

if(!dir.exists("reference")) {
    dir.create("reference")
}
download.file(
    "https://zenodo.org/records/7779017/files/pbmc_multimodal_2023.rds?download=1",
    "reference/pbmc_multimodal_2023.rds"
)

In [3]:

reference <- readRDS("reference/pbmc_multimodal_2023.rds")

Between Level 2 (L2) and Level 3 (L3) labels provided by the Satija lab, we like to add an additional level that we call L2.5, which separates the Treg Naive and Memory cells based on L3, and assigns a CD8 TEMRA cell label to cells with the L3 labels CD8 TEM_4 and CD8 TEM_5. All other cell types use their L2 assignments.

In [4]:

l3 <- as.character(reference@meta.data$celltype.l3)
l2 <- as.character(reference@meta.data$celltype.l2)
l2.5 <- l2
l2.5[l3 == "Treg Naive"] <- "Treg Naive"
l2.5[l3 == "Treg Memory"] <- "Treg Memory"
l2.5[l3 %in% c("CD8 TEM_4", "CD8 TEM_5")] <- "CD8 TEMRA"

reference <- AddMetaData(reference, metadata = l2.5, col.name = "celltype.l2.5")

Retreive sample metadata¶

In an earlier step, we assembled and stored sample metadata in HISE. We'll pull this file, and use it to retrieve file for our labeling process.

In [5]:

sample_meta_uuid <- "2da66a1a-17cc-498b-9129-6858cf639caf"

In [6]:

res <- cacheFiles(list(sample_meta_uuid))
sample_meta_file <- list.files(
    paste0("cache/", sample_meta_uuid), 
    pattern = ".csv",
    full.names = TRUE
)

In [7]:

hise_meta <- read.csv(sample_meta_file)

Cache input files¶

Next, we'll use the hise package to cache all of the input files. With this many files, there are occasional problems with transfer, so we'll also add a check for existing files in our function, and run a second pass to make sure we have everything we want to label.

In [ ]:

# Helper function with cache path check
cache_file <- function(h5_uuid) {
    cache_dir <- paste0("cache/",h5_uuid)
     if(!dir.exists(cache_dir)) {
         res <- cacheFiles(list(h5_uuid))
     }
}

# Walk file UUIDs to cache
file_ids <- hise_meta$file.id
walk(file_ids, cache_file)

Run a second pass to make sure we have everything

In [ ]:

walk(file_ids, cache_file)

Divide data into chunks for parallel processing¶

For labeling, we'll take files in batches of up to 10 files. We'll label those files, then output the results for each sample in the batch.

In [8]:

hise_meta <- hise_meta %>% 
  arrange(file.batchID)

In [9]:

b <- rep(1:11, each = 10)[1:nrow(hise_meta)]
df_chunk_list <- split(hise_meta, b)

In [10]:

map_int(df_chunk_list, nrow)

1: 10
2: 10
3: 10
4: 10
5: 10
6: 10
7: 10
8: 10
9: 10
10: 10
11: 8

Prepare output directory¶

We'll store results for each sample in output/Hao_PBMC/.

In [11]:

out_dir <- "output/Hao_PBMC/"
if(!dir.exists(out_dir)) {
    dir.create(out_dir, recursive = TRUE)
}

Functions for parallel label transfer of chunks¶

For each chunk, we'll retrieve the data from HISE, perform label transfer and ADT imputation using Seurat, and then store the labels and imputed ADT matrices to allow us to assess cell type identity.

The main function to perform these steps is label_chunk(), provided below. There are also 4 helper functions that assist in performing these steps for each sample:
read_so() Reads the .h5 files stored in HISE as Seurat Objects for analysis
write_labels() Writes the labeling results to .csv files for each sample
get_sample_adt() subsets the ADT matrix for each sample and returns ADT values per sample
write_adt() Writes the imputed ADT matrix values to a .h5 file for later use

In [12]:

read_so <- function(h5_uuid) {
    cache_dir <- paste0("cache/",h5_uuid)
    h5_file <- list.files(
        cache_dir, 
        pattern = ".h5", 
        full.names = TRUE
    )
    
    counts <- read_h5_dgCMatrix(h5_file)
    rownames(counts) <- make.unique(rownames(counts))
    meta <- read_h5_cell_meta(h5_file)
    rownames(meta) <- meta$barcodes

    so <- CreateSeuratObject(
      counts,
      meta.data = meta,
      assay = "RNA")

    so
}

write_labels <- function(sample_labels, sample_id) {
    out_file <- paste0("output/Hao_PBMC/", sample_id, "_Hao_PBMC.csv")
    write.csv(
        sample_labels,
        out_file,
        row.names = FALSE,
        quote = FALSE
    )
}

get_sample_adt <- function(sample_meta, adt_data) {
    adt_data[,sample_meta$barcodes]
}

write_adt <- function(sample_adt, sample_id) {
    out_file <- paste0("output/Hao_PBMC/", sample_id, "_ADT.h5")
    list_mat <- list(
        i = sample_adt@i,
        p = sample_adt@p,
        x = sample_adt@x,
        Dim = dim(sample_adt),
        rownames = rownames(sample_adt),
        colnames = colnames(sample_adt)
    )

    h5createFile(out_file)
    h5write(list_mat, out_file, "mat")
}

check_output <- function(sample_id) {
    out_file <- paste0("output/Hao_PBMC/", sample_id, "_Hao_PBMC.csv")
    file.exists(out_file)
}

label_chunk <- function(meta_data) {
    # check for existing labels
    sample_ids <- unique(meta_data$pbmc_sample_id)
    out_check <- map_lgl(sample_ids, check_output)
    if(sum(out_check) == length(sample_ids)) {
        return(invisible(NULL))
    }
    
    # Read and combine individual samples
    so_list <- map(meta_data$file.id, read_so)
    combined <- Reduce(merge, so_list)
    rm(so_list)

    # Transform data to match the reference
    combined <- SCTransform(
        combined,
        method = "glmGamPoi", 
        verbose = FALSE
    )
    
    # find anchors
    anchors <- FindTransferAnchors(
      reference = reference,
      query = combined,
      normalization.method = "SCT",
      reference.reduction = "spca",
      dims = 1:50
    )  
        
    #perform projection to get labels
    combined <- MapQuery(
      anchorset = anchors,
      query = combined,
      reference = reference,
      refdata = list(
        celltype.l1 = "celltype.l1",
        celltype.l2 = "celltype.l2",
        celltype.l3 = "celltype.l3",
        celltype.l2.5 = "celltype.l2.5",
        predicted_ADT = "ADT"
      ),
      reference.reduction = "spca", 
      reduction.model = "wnn.umap"
    )

    # Split the metadata by sample for output
    sample_meta_list <- split(combined@meta.data, combined@meta.data$pbmc_sample_id)

    # Write labels
    walk2(
        sample_meta_list, names(sample_meta_list),
        write_labels
    )

    # Subset and write projected ADT data
    sample_adt_list <- map(
        sample_meta_list, 
        get_sample_adt, 
        adt_data = combined@assays$predicted_ADT@data)

    walk2(
        sample_adt_list, names(sample_meta_list),
        write_adt
    )
}

Apply label transfer to all chunks¶

We'll use mclapply() to perform our label transfer steps in parallel for multiple chunks.

Because of the extremely memory-intensive use of SCTransform(), we're limited in the number of chunks that we can process simultaneously.

In [13]:

res <- mclapply(
    df_chunk_list,
    label_chunk,
    mc.cores = 3
)

Assemble results¶

Next, we'll assemble all of the labeling results in a single file for storage in HISE and downstream utilization.

In [14]:

label_files <- list.files(
    out_dir,
    pattern = ".csv",
    full.names = TRUE
)

In [15]:

all_labels <- map_dfr(label_files, read.csv)

In [16]:

all_labels <- all_labels %>%
  select(pbmc_sample_id, barcodes, starts_with("predicted"))

In [18]:

out_csv <- paste0(
    "output/ref_seurat_labels_PBMC_",
    Sys.Date(),
    ".csv"
)
write.csv(
    all_labels,
    out_csv,
    row.names = FALSE,
    quote = FALSE
)

In [19]:

adt_files <- list.files(
    out_dir,
    pattern = ".h5",
    full.names = TRUE
)

In [20]:

read_adt_mat <- function(adt_file) {
    mat_list <- rhdf5::h5read(adt_file, "/mat")
    mat_list <- lapply(mat_list, as.vector)
    mat <- Matrix::sparseMatrix(
        i = mat_list$i,
        p = mat_list$p,
        x = mat_list$x,
        dims = mat_list$Dim,
        dimnames = list(
            mat_list$rownames,
            mat_list$colnames
        ),
        index1 = FALSE
    )
    mat
}

In [21]:

adt_list <- map(adt_files, read_adt_mat)

In [22]:

all_adt <- do.call(cbind, adt_list)

In [23]:

h5_list <- list(
    matrix_dgCMatrix = all_adt
)

In [24]:

h5_list <- h5_list_convert_from_dgCMatrix(h5_list)

In [25]:

str(h5_list)

List of 1
 $ matrix:List of 6
  ..$ data    : num [1:475230846] 0.406 0.7 0.973 2.481 1.206 ...
  ..$ indices : int [1:475230846] 0 1 2 3 4 5 6 7 8 9 ...
  ..$ indptr  : int [1:2093079] 0 227 454 681 908 1135 1362 1589 1817 2044 ...
  ..$ shape   : int [1:2] 228 2093078
  ..$ barcodes: chr [1:2093078] "cf71f47048b611ea8957bafe6d70929e" "cf71f54248b611ea8957bafe6d70929e" "cf71fa1048b611ea8957bafe6d70929e" "cf71fb7848b611ea8957bafe6d70929e" ...
  ..$ features:List of 1
  .. ..$ id: chr [1:228] "CD39" "Rat-IgG1-1" "CD107a" "CD62P" ...

In [26]:

out_h5 <- paste0(
    "output/ref_seurat_projected-ADT_PBMC_",
    Sys.Date(),
    ".h5"
)
write_h5_list(
    h5_list,
    out_h5
)

Store results in HISE¶

Finally, we store the output file in our Collaboration Space for later retrieval and use. We need to provide the UUID for our Collaboration Space (aka studySpaceId), as well as a title for this step in our analysis process.

The hise function uploadFiles() also requires the FileIDs from the original fileset for reference.

In [30]:

study_space_uuid <- "64097865-486d-43b3-8f94-74994e0a72e0"
title <- paste("Ref. Seurat Label Predictions", Sys.Date())

In [28]:

in_list <- as.list(c(sample_meta_uuid, hise_meta$file.id))

In [29]:

out_list <- as.list(c(out_csv, out_h5))

In [31]:

uploadFiles(
    files = out_list,
    studySpaceId = study_space_uuid,
    title = title,
    inputFileIds = in_list,
    store = "project",
    doPrompt = FALSE
)

[1] "Authorization token invalid or expired."
[1] "Retrying..."

$files

'output/ref_seurat_labels_PBMC_2024-02-18.csv'
'output/ref_seurat_projected-ADT_PBMC_2024-02-18.h5'

$traceId

'e9ecdbe0-8f06-4c88-ac02-457f26c7abe1'

In [32]:

sessionInfo()

R version 4.3.2 (2023-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS/LAPACK: /opt/conda/lib/libopenblasp-r0.3.25.so;  LAPACK version 3.11.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] Seurat_5.0.1       SeuratObject_5.0.1 sp_2.1-2           purrr_1.0.2       
 [5] hise_2.16.0        H5weaver_1.2.0     rhdf5_2.46.1       Matrix_1.6-4      
 [9] data.table_1.14.10 dplyr_1.1.4       

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3     jsonlite_1.8.8         magrittr_2.0.3        
  [4] spatstat.utils_3.0-4   vctrs_0.6.5            ROCR_1.0-11           
  [7] spatstat.explore_3.2-5 RCurl_1.98-1.14        base64enc_0.1-3       
 [10] htmltools_0.5.7        curl_5.1.0             Rhdf5lib_1.24.1       
 [13] sctransform_0.4.1      parallelly_1.36.0      KernSmooth_2.23-22    
 [16] htmlwidgets_1.6.4      ica_1.0-3              plyr_1.8.9            
 [19] plotly_4.10.3          zoo_1.8-12             uuid_1.2-0            
 [22] igraph_1.6.0           mime_0.12              lifecycle_1.0.4       
 [25] pkgconfig_2.0.3        R6_2.5.1               fastmap_1.1.1         
 [28] fitdistrplus_1.1-11    future_1.33.1          shiny_1.8.0           
 [31] digest_0.6.34          colorspace_2.1-0       patchwork_1.1.3       
 [34] tensor_1.5             RSpectra_0.16-1        irlba_2.3.5.1         
 [37] progressr_0.14.0       fansi_1.0.6            spatstat.sparse_3.0-3 
 [40] httr_1.4.7             polyclip_1.10-6        abind_1.4-5           
 [43] compiler_4.3.2         withr_3.0.0            fastDummies_1.7.3     
 [46] MASS_7.3-60            tools_4.3.2            lmtest_0.9-40         
 [49] httpuv_1.6.13          future.apply_1.11.1    goftest_1.2-3         
 [52] glue_1.7.0             nlme_3.1-164           rhdf5filters_1.14.1   
 [55] promises_1.2.1         grid_4.3.2             pbdZMQ_0.3-10         
 [58] Rtsne_0.17             cluster_2.1.6          reshape2_1.4.4        
 [61] generics_0.1.3         gtable_0.3.4           spatstat.data_3.0-3   
 [64] tidyr_1.3.0            utf8_1.2.4             spatstat.geom_3.2-7   
 [67] RcppAnnoy_0.0.21       ggrepel_0.9.4          RANN_2.6.1            
 [70] pillar_1.9.0           stringr_1.5.1          spam_2.10-0           
 [73] IRdisplay_1.1          RcppHNSW_0.5.0         later_1.3.2           
 [76] splines_4.3.2          lattice_0.22-5         survival_3.5-7        
 [79] deldir_2.0-2           tidyselect_1.2.0       miniUI_0.1.1.1        
 [82] pbapply_1.7-2          gridExtra_2.3          scattermore_1.2       
 [85] matrixStats_1.2.0      stringi_1.8.3          lazyeval_0.2.2        
 [88] evaluate_0.23          codetools_0.2-19       tibble_3.2.1          
 [91] cli_3.6.2              uwot_0.1.16            IRkernel_1.3.2        
 [94] xtable_1.8-4           reticulate_1.34.0      repr_1.1.6.9000       
 [97] munsell_0.5.0          Rcpp_1.0.12            globals_0.16.2        
[100] spatstat.random_3.2-2  png_0.1-8              ellipsis_0.3.2        
[103] ggplot2_3.4.4          assertthat_0.2.1       dotCall64_1.1-1       
[106] bitops_1.0-7           listenv_0.9.0          viridisLite_0.4.2     
[109] scales_1.3.0           ggridges_0.5.5         leiden_0.4.3.1        
[112] crayon_1.5.2           rlang_1.1.3            cowplot_1.1.2