Comparing totalVI and OSCA book CITE-seq analyses (original) (raw)

1 Overview
2 Technical steps to facilitate comparison
3 Key outputs of the chapter 12 analysis
4 Analogs to the chapter 12 findings, based on totalVI quantifications
5 Conclusions
6 Session information

Overview

This vignette endeavors to put Bioconductor and scvi-tools together to help understand how different data structures and methods relevant to CITE-seq analysis contribute to interpretation of CITE-seq exeperiments.

The scvi-tools tutorial (for version 0.20.0) analyzes a pair of 10x PBMC CITE-seq experiments (5k and 10k cells). Chapter 12 of the OSCA book analyzes only the 10k dataset.

Technical steps to facilitate comparison

The following subsections are essentially “code-only”. We exhibit steps necessary to assemble a SingleCellExperiment instance with the subset of the totalVI quantifications produced for the cells from the “10k” dataset.

Acquire software

library(SingleCellExperiment)
library(scater)
library(scviR)

Obtain key data representations

ch12sce = getCh12Sce(clear_cache=FALSE)
ch12sce

## class: SingleCellExperiment 
## dim: 33538 7472 
## metadata(2): Samples se.averaged
## assays(2): counts logcounts
## rownames(33538): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
##   ENSG00000268674
## rowData names(3): ID Symbol Type
## colnames(7472): AAACCCAAGATTGTGA-1 AAACCCACATCGGTTA-1 ...
##   TTTGTTGTCGAGTGAG-1 TTTGTTGTCGTTCAGA-1
## colData names(3): Sample Barcode sizeFactor
## reducedDimNames(0):
## mainExpName: Gene Expression
## altExpNames(1): Antibody Capture

options(timeout=3600)
fullvi = getTotalVI5k10kAdata()
fullvi

## AnnData object with n_obs × n_vars = 10849 × 4000
##     obs: 'n_genes', 'percent_mito', 'n_counts', 'batch', '_scvi_labels', '_scvi_batch', 'leiden_totalVI'
##     var: 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm', 'highly_variable_nbatches'
##     uns: '_scvi_manager_uuid', '_scvi_uuid', 'hvg', 'leiden', 'log1p', 'neighbors', 'umap'
##     obsm: 'X_totalVI', 'X_umap', 'denoised_protein', 'protein_expression', 'protein_foreground_prob'
##     layers: 'counts', 'denoised_rna'
##     obsp: 'connectivities', 'distances'

Assemble a SingleCellExperiment with totalVI outputs

Acquire cell identities and batch labels

totvi_cellids = rownames(fullvi$obs)
totvi_batch = fullvi$obs$batch

Acquire quantifications and latent space positions

totvi_latent = fullvi$obsm$get("X_totalVI")
totvi_umap = fullvi$obsm$get("X_umap")
totvi_denoised_rna = fullvi$layers$get("denoised_rna")
totvi_denoised_protein = fullvi$obsm$get("denoised_protein")
totvi_leiden = fullvi$obs$leiden_totalVI

Drop 5k data from all

is5k = which(totvi_batch == "PBMC5k")
totvi_cellids = totvi_cellids[-is5k]
totvi_latent = totvi_latent[-is5k,]
totvi_umap = totvi_umap[-is5k,]
totvi_denoised_rna = totvi_denoised_rna[-is5k,]
totvi_denoised_protein = totvi_denoised_protein[-is5k,]
totvi_leiden = totvi_leiden[-is5k]

Label the rows of components

rownames(totvi_latent) = totvi_cellids
rownames(totvi_umap) = totvi_cellids
rownames(totvi_denoised_rna) = totvi_cellids
rownames(totvi_denoised_protein) = totvi_cellids
names(totvi_leiden) = totvi_cellids

Find common cell ids

In this section we reduce the cell collections to cells common to the chapter 12 and totalVI datasets.

comm = intersect(totvi_cellids, ch12sce$Barcode)

Build the totalVI SingleCellExperiment

# select and order
totvi_latent = totvi_latent[comm,]
totvi_umap = totvi_umap[comm,]
totvi_denoised_rna = totvi_denoised_rna[comm,]
totvi_denoised_protein = totvi_denoised_protein[comm,]
totvi_leiden = totvi_leiden[comm]
 
# organize the totalVI into SCE with altExp

totsce = SingleCellExperiment(SimpleList(logcounts=t(totvi_denoised_rna))) # FALSE name
rowData(totsce) = S4Vectors::DataFrame(fullvi$var)
rownames(totsce) = rownames(fullvi$var)
rowData(totsce)$Symbol = rownames(totsce)
nn = SingleCellExperiment(SimpleList(logcounts=t(totvi_denoised_protein))) # FALSE name
reducedDims(nn) = list(UMAP=totvi_umap)
altExp(totsce) = nn
altExpNames(totsce) = "denoised_protein"
totsce$leiden = totvi_leiden
altExp(totsce)$leiden = totvi_leiden
altExp(totsce)$ch12.clusters = altExp(ch12sce[,comm])$label

# add average ADT abundance to metadata, for adt_profiles

tot.se.averaged <- sumCountsAcrossCells(altExp(totsce), altExp(totsce)$leiden,
    exprs_values="logcounts", average=TRUE)
rownames(tot.se.averaged) = gsub("_TotalSeqB", "", rownames(tot.se.averaged))
metadata(totsce)$se.averaged = tot.se.averaged

Reduce the chapter 12 dataset to the cells held in common

colnames(ch12sce) = ch12sce$Barcode
ch12sce_matched = ch12sce[, comm]

Key outputs of the chapter 12 analysis

Clustering and projection based on the ADT quantifications

The TSNE projection of the normalized ADT quantifications and the walktrap cluster assignments are produced for the cells common to the two approaches.

plotTSNE(altExp(ch12sce_matched), color_by="label", text_by="label")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text_repel()`).

Cluster profiles based on averaging ADT quantities within clusters

This heatmap uses precomputed cluster averages that are lodged in the metadata element of the SingleCellExperiment. Colors represent the log2-fold change from the grand average across all clusters.

adtProfiles(ch12sce_matched)

Marker expression patterns in mRNA-based sub-clusters of ADT-based clusters

We enhance the annotation of the list of subclusters retrieved using getCh12AllSce and then drill into mRNA-based subclusters of ADT-based cluster 3 to compare expression levels of three genes.

ch12_allsce = getCh12AllSce() 
ch12_allsce = lapply(ch12_allsce, function(x) {
   colnames(x)= x$Barcode; 
   cn = colnames(x); 
   x = x[,intersect(cn,comm)]; x})
of.interest <- "3"
markers <- c("GZMH", "IL7R", "KLRB1")
plotExpression(ch12_allsce[[of.interest]], x="subcluster",
    features=markers, swap_rownames="Symbol", ncol=3)

There is a suggestion of a boolean basis for subcluster identity, depending on low or high expression of the selected genes.

Graduated relationships between mRNA and surface protein expression

Following the exploration in OSCA chapter 12, cluster 3 is analyzed for a regression association between expression measures of three genes and the ADT-based abundance of CD127.

plotExpression(ch12_allsce[["3"]], x="CD127", show_smooth=TRUE, show_se=FALSE, 
    features=c("IL7R", "TPT1", "KLRB1", "GZMH"), swap_rownames="Symbol")

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Analogs to the chapter 12 findings, based on totalVI quantifications

The Leiden clustering in UMAP projection

plotUMAP(altExp(totsce), color_by="leiden", text_by="leiden")

Cluster profiles based on average ADT abundance, using denoised protein quantifications

The approach to profiling the ADT abundances used in the totalVI tutorial employs scaling to (0,1).

tav = S4Vectors::metadata(totsce)$se.averaged
ata = assay(tav)
uscale = function(x) (x-min(x))/max(x)
scmat = t(apply(ata,1,uscale))
pheatmap::pheatmap(scmat, cluster_rows=FALSE)

Concordance in ADT-based clustering between OSCA and totalVI

A quick view of the concordance of the two clustering outcomes is

atot = altExp(totsce)
ach12 = altExp(ch12sce_matched)
tt = table(ch12=ach12$label, VI=atot$leiden)
pheatmap::pheatmap(log(tt+1))

With this we can pick out some clusters with many cells in common:

lit = tt[c("9", "12", "5", "3"), c("0", "1", "2", "8", "6", "5")]
rownames(lit) = sQuote(rownames(lit))
colnames(lit) = sQuote(colnames(lit))
lit

##       VI
## ch12    '0'  '1'  '2'  '8'  '6'  '5'
##   '9'  1334    0    0    0    0    0
##   '12'    0  993    0    0   15    0
##   '5'     0  102  671    8   44    2
##   '3'     0    0    5  322  297   67

Subcluster assessment for OSCA cluster “3”

Let’s examine the distributions of marker mRNAs in the Leiden totalVI clusters corresponding to OSCA cluster “3”:

tsub = totsce[,which(altExp(totsce)$leiden %in% c("5", "6", "8"))]
markers <- c("GZMH", "IL7R", "KLRB1")
altExp(tsub)$leiden = factor(altExp(tsub)$leiden) # squelch unused levels
tsub$leiden = factor(tsub$leiden) # squelch unused levels
plotExpression(tsub, x="leiden",
    features=markers, swap_rownames="Symbol", ncol=3)

Note that the y axis label is incorrect – we are plotting the denoised expression values from totalVI.

The display seems roughly consistent with the “boolean basis” observed above with the mRNA-based subclustering.

Graduated relationships between ADT and mRNA abundance as measured by totalVI

The same approach is taken as above. We don’t have TPT1 in the 4000 genes retained in the totalVI exercise.

rn = rownames(altExp(tsub))
rn = gsub("_TotalSeqB", "", rn)
rownames(altExp(tsub)) = rn
rowData(altExp(tsub)) = DataFrame(Symbol=rn)
plotExpression(tsub, x="CD127", show_smooth=TRUE, show_se=FALSE,
   features=c("IL7R", "KLRB1", "GZMH"), swap_rownames="Symbol")

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Conclusions

We have shown how rudimentary programming and data organization can be used to make outputs of OSCA and totalVI methods amenable to comparison in the Bioconductor framework.

The scviR package includes a shiny app in the functionexplore_subcl that should be expanded to facilitate exploration of totalVI subclusters. Much more work remains to be done in the area of exploring

additional approaches to integrative interpretation of ADT and mRNA abundance patterns, such as intersection and concatenation methods in the feature selection materials in OSCA ch. 12
effects of tuning and architecture details for the totalVI VAE

Session information

## R version 4.5.0 RC (2025-04-04 r88126)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.2 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] scater_1.36.0               scuttle_1.18.0             
##  [3] reshape2_1.4.4              ggplot2_3.5.2              
##  [5] scviR_1.8.0                 SingleCellExperiment_1.30.0
##  [7] SummarizedExperiment_1.38.0 Biobase_2.68.0             
##  [9] GenomicRanges_1.60.0        GenomeInfoDb_1.44.0        
## [11] IRanges_2.42.0              S4Vectors_0.46.0           
## [13] BiocGenerics_0.54.0         generics_0.1.3             
## [15] MatrixGenerics_1.20.0       matrixStats_1.5.0          
## [17] shiny_1.10.0                basilisk_1.20.0            
## [19] reticulate_1.42.0           BiocStyle_2.36.0           
## 
## loaded via a namespace (and not attached):
##  [1] DBI_1.2.3               gridExtra_2.3           rlang_1.1.6            
##  [4] magrittr_2.0.3          compiler_4.5.0          RSQLite_2.3.9          
##  [7] mgcv_1.9-3              dir.expiry_1.16.0       png_0.1-8              
## [10] vctrs_0.6.5             stringr_1.5.1           pkgconfig_2.0.3        
## [13] crayon_1.5.3            fastmap_1.2.0           magick_2.8.6           
## [16] dbplyr_2.5.0            XVector_0.48.0          labeling_0.4.3         
## [19] promises_1.3.2          rmarkdown_2.29          UCSC.utils_1.4.0       
## [22] ggbeeswarm_0.7.2        tinytex_0.57            purrr_1.0.4            
## [25] bit_4.6.0               xfun_0.52               cachem_1.1.0           
## [28] beachmat_2.24.0         jsonlite_2.0.0          blob_1.2.4             
## [31] later_1.4.2             DelayedArray_0.34.0     BiocParallel_1.42.0    
## [34] irlba_2.3.5.1           parallel_4.5.0          R6_2.6.1               
## [37] stringi_1.8.7           bslib_0.9.0             RColorBrewer_1.1-3     
## [40] limma_3.64.0            jquerylib_0.1.4         Rcpp_1.0.14            
## [43] bookdown_0.43           knitr_1.50              splines_4.5.0          
## [46] httpuv_1.6.15           Matrix_1.7-3            tidyselect_1.2.1       
## [49] abind_1.4-8             yaml_2.3.10             viridis_0.6.5          
## [52] codetools_0.2-20        curl_6.2.2              lattice_0.22-7         
## [55] tibble_3.2.1            plyr_1.8.9              basilisk.utils_1.20.0  
## [58] withr_3.0.2             evaluate_1.0.3          BiocFileCache_2.16.0   
## [61] pillar_1.10.2           BiocManager_1.30.25     filelock_1.0.3         
## [64] munsell_0.5.1           scales_1.3.0            xtable_1.8-4           
## [67] glue_1.8.0              pheatmap_1.0.12         tools_4.5.0            
## [70] BiocNeighbors_2.2.0     ScaledMatrix_1.16.0     cowplot_1.1.3          
## [73] grid_4.5.0              colorspace_2.1-1        nlme_3.1-168           
## [76] GenomeInfoDbData_1.2.14 beeswarm_0.4.0          BiocSingular_1.24.0    
## [79] vipor_0.4.7             cli_3.6.4               rsvd_1.0.5             
## [82] S4Arrays_1.8.0          viridisLite_0.4.2       dplyr_1.1.4            
## [85] gtable_0.3.6            sass_0.4.10             digest_0.6.37          
## [88] SparseArray_1.8.0       ggrepel_0.9.6           farver_2.1.2           
## [91] memoise_2.0.1           htmltools_0.5.8.1       lifecycle_1.0.4        
## [94] httr_1.4.7              statmod_1.5.0           mime_0.13              
## [97] bit64_4.6.0-1