Visium: cell type deconvolution for multi-cellular spots

Highlights

Interpretable cell type deconvolution;
Healthy reference annotates cancerous query;
Insightful cell type colocalisation analysis.

Introduction

In this case study, we will apply PhiSpace to cell type deconvolution:

Query: Visium dataset consisting of samples from non-small cell lung cancer (NSCLC) patients (De Zuani et al., 2024).
Reference: the Human Lung Cell Atlas (HLCA) scRNA-seq data (Sikkema et al., 2023).

10x Visium is a supercellular spatial transcriptomics (ST) platform. The minimal measurement unit for the Visium platform are 55μm spots, where each spot potentially contains multiple cells. In this context, cell type deconvolution refers to the statistical task of estimating the cell type abundance in each spot.

Download the complete Visium dataset from the shared folder and preserve its directory layout. It includes Visium_utils.R, the qs2 data objects, and cached annotation results.

Set PHISPACE_VISIUM_DATA to your data folder before running this vignette, for example Sys.setenv(PHISPACE_VISIUM_DATA = "/path/to/Visium"). Place Visium_utils.R in that folder and the data files at the relative paths shown below. This vignette uses qs2::qs_read() and qs2::qs_save() with migrated .qs2 files. The old .qs files cannot be converted by renaming them.

# Name of the game
suppressPackageStartupMessages(library(PhiSpace))
suppressPackageStartupMessages(library(vizOmics))
# Tidyverse packages
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(magrittr))
suppressPackageStartupMessages(library(ggpubr))
suppressPackageStartupMessages(library(tidyr))
# Other utils
suppressPackageStartupMessages(library(qs2)) # Fast read and write of large R objects
suppressPackageStartupMessages(library(ComplexHeatmap)) # plot heatmap

dat_dir <- Sys.getenv("PHISPACE_VISIUM_DATA", unset = "data/Visium")
if (!file.exists(file.path(dat_dir, "Visium_utils.R"))) {
  stop("Set PHISPACE_VISIUM_DATA to the data folder containing Visium_utils.R.")
}
dat_dir <- paste0(normalizePath(dat_dir, winslash = "/"), "/")
dir.create(file.path(dat_dir, "output"), recursive = TRUE, showWarnings = FALSE)

source(paste0(dat_dir, "Visium_utils.R"), local = TRUE)

Data preparation

The shared folder contains:

A 10% subset of the HLCA reference.
A list of 40 healthy and cancerous lung Visium samples. We use 18 samples: two healthy samples and one background and one cancerous sample from each of eight patients.
PhiSpace feature-importance scores and cached annotation results.

impScores <- qs_read(paste0(dat_dir, "output/refImpScores.qs2"))
selectedFeat <- selectFeat(impScores, 500)$selectedFeat
reference <- qs_read(paste0(dat_dir, "data/LungRef/AzimuthLung2.0_sce_0.1sub.qs2"))
# Retain only the reference assay and genes used by the model to reduce RAM use.
assays(reference) <- assays(reference)["log1p"]
reference <- reference[intersect(rownames(reference), selectedFeat), ]
invisible(gc())
query_list <- qs_read(paste0(dat_dir, "data/Visium_NSCLC/allTissues.qs2"))
selectedTissues <- tissueNames_Visium[1:18]
stopifnot(all(selectedTissues %in% names(query_list)))
query_list <- lapply(query_list[selectedTissues], function(x) {
  # Counts are needed for marker analysis; log1p is used for annotation.
  assays(x) <- assays(x)[c("counts", "log1p")]
  x
})
invisible(gc())

Visualise a cancerous sample P11_T3. Here we look at two markers: KRT19 (also known as CYFRA21-1), a NSCLC marker and CD79A, a B cell marker.

VisTissueName <- "P11_T3"
qu_vis <- query_list[[VisTissueName]]
VizSpatial(qu_vis, feature = "KRT19", assay2use = "log1p")

VizSpatial(qu_vis, feature = "CD79A", assay2use = "log1p")

Run PhiSpace

Run PhiSpace using all the genes can be RAM-demanding, hence we select 500 genes that are most useful for predicting each cell type. This gives us 4613 genes. PhiSpace() annotates all 18 samples with one model and stores the normalised scores of each sample in reducedDim(query_list[[sample]], "PhiSpace"). The shared data folder also provides these scores (output/PhiSpaceScores.qs2) for systems without enough RAM to recompute the annotation.

YtrainName <- "ann_finest_level"
length(selectedFeat)

## [1] 4613

PhiScPath <- paste0(dat_dir, "output/PhiSpaceScores.qs2")
if(!file.exists(PhiScPath)){
  query_list <- PhiSpace(
    reference,
    query_list,
    phenotypes = YtrainName,
    refAssay = "log1p",
    selectedFeat = selectedFeat,
    regMethod = "PLS",
    center = TRUE,
    scale = FALSE
  )
  qs_save(lapply(query_list, reducedDim, "PhiSpace"), PhiScPath)
} else {
  PhiScores <- qs_read(PhiScPath)
  for(tissue in names(query_list)){
    reducedDim(query_list[[tissue]], "PhiSpace") <- PhiScores[[tissue]]
  }
}

As an example, we visualise the predicted B cell abundance in sample P11_T3.

qu_vis <- query_list[["P11_T3"]]
VizSpatial(qu_vis, reducedDim = "B cells", reducedDim2use = "PhiSpace", censor = T)

How about the performance of other deconvolution methods? How to compare their performances when there is no ground truth? See the Visium case study in our PhiSpace ST paper for more details (spoiler alert: PhiSpace was the best).

Multi-sample cell type co-localisation analysis

Obtaining interpretable cell type deconvolution results is only the first steps towards discovery of spatial biology. We showcase one type of analysis that can be done based on PhiSpace deconvolution - cell type colocalisation. The main idea is that, from each sample, we calculate a cell type correlation matrix (called co-presence matrix in our PhiSpace ST paper), representing how cell types tend to colocalised. Then we apply dimension reduction methods to these matrices (vectorised) to visualise them.

First we calcualte the sample-specific correlation matrices and vectorise them.

VisTissueNames <- names(query_list)
corMat_list <- vector("list", length(VisTissueNames)) %>% `names<-`(VisTissueNames)
for(ii in 1:length(VisTissueNames)){
  VisTissueName <- VisTissueNames[ii]
  corMat <- cor(reducedDim(query_list[[VisTissueName]], "PhiSpace"))
  corMat_list[[ii]] <- corMat
}
cTypeCombo <- outer(
  rownames(corMat_list[[1]]), colnames(corMat_list[[1]]),
  function(x, y) paste0(x, "<->", y)
)
corMatVec <- sapply(corMat_list, function(x) as.vector(x[lower.tri(x)])) %>% 
  `colnames<-`(VisTissueNames) %>%
  `rownames<-`(cTypeCombo[lower.tri(cTypeCombo)]) %>% t()

Next we simply compute a PCA of all correlation matrices. It turns out that a simple PCA suffices here as the signal is very strong.

set.seed(92394)
pca_res <- getPC(corMatVec, ncomp = 2, center = T, scale = F)

We visualise the first two PCs. As we can see, all samples were well separated into groups according to disease conditions:

Healthy: samples from healthy donors;
Background: non-cancerous tissue from cancerous samples;
LUAD: lung adenocarcinoma;
LUSC: lung squamous cell carcinoma.

tempPCAplot(pca_res) + expand_limits(x = c(-5.5, 7.5))

plotLoadings(pca_res$loadings, nfeat = 20) + theme(axis.text.x = element_text(angle = 45))

In addition to the separation of patients, we can also observe that LUSC samples tend to have the largest PC1 values, healthy and background the lowest, and LUAD samples having intermediate PC1 values. Given that LUSC tends to have worse prognosis than LUAD, layout of the PC plot about seems to suggest that PC1 represent prognosis to some degree. Lacking detailed information about patient prognosis from the original publication, we opt for a data-driven way to find out.

According to Chen & Dhahbi (2021), CELSR2 is a marker for better prognosis in LUSC; KRT17, KRT6A, S100A2, TRIM29, REPS1, and GPC1 are markers for worse prognosis in LUAD; PERP, ELFN2, ARHGAP12 and QSOX1 are markers for better prognosis in LUSC. Based on how the favourable prognosis markers are upregulated (see definition of meanEsts in Visium_utils.R), we calculated the prognosis scores for all patients.

Note that since the markers are different for LUSC and LUAD, we cannot directly compare prognosis scores between LUSC and LUAD patients. However, the scores are still informative. For example, P24 was predicted to have the worst prognosis (lowest favourable prognosis score) among LUAD patients, which probably explained why it was located closer to LUSC patients on the PCA scatter plot above.

# Common gene of all Visium slides
commonGeneNames <- Reduce(intersect, lapply(query_list, function(x) rownames(x)))
mk_list <- list(
  LUSC_good = c("CELSR2"),
  LUAD_bad = c("KRT17", "KRT6A", "S100A2", "TRIM29", "REPS1", "GPC1"),
  LUAD_good = c("PERP", "ELFN2", "ARHGAP12", "QSOX1")
)
mk_list <- lapply(mk_list, function(x) intersect(x, commonGeneNames))
# Function for calculate mean expression level of marker genes in all tissues
assay2use <- "counts"
meanGeneLvls <- sapply(query_list[1:18], function(x) mean(rowSums(assay(x, assay2use))))
mkLvls_LUSC <- data.frame(favourable = meanEsts(mk_list[["LUSC_good"]])) 
mkLvls_LUAD_good <- data.frame(favourable = meanEsts(mk_list[["LUAD_good"]]))
mkLvls_LUAD_bad <- data.frame(favourable = meanEsts(mk_list[["LUAD_bad"]]))
mkLvls_LUAD <- mkLvls_LUAD_good/mkLvls_LUAD_bad
data.frame(
  favourable = c(
    mkLvls_LUSC[Visium_cancerTypes == "LUSC",],
    mkLvls_LUAD[Visium_cancerTypes == "LUAD",]
  ),
  sample = c(
    VisTissueNames[Visium_cancerTypes == "LUSC"],
    VisTissueNames[Visium_cancerTypes == "LUAD"]
  ), 
  cancerType = c(
    Visium_cancerTypes[Visium_cancerTypes == "LUSC"],
    Visium_cancerTypes[Visium_cancerTypes == "LUAD"]
  )
) %>% 
  ggplot() + geom_bar(aes(sample, favourable, fill = cancerType), stat = "identity") +
  scale_fill_manual(values = VisiumCancerCols) + theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 90), axis.title.x = element_blank()) + 
  facet_wrap(~cancerType, scales = "free", ncol = 1) +
  theme(legend.position = "none", axis.title.y = element_blank())

To focus on the role of immune cell types, we repeat the above analysis using only immune cell types. The separation of disease conditions was largely retained. Note how plasmacytoid dendritic cells (pDCs) played an important role in explaining the separation of samples along PC1.

immuneDF <- data.frame(
  celltype = rownames(corMat_list[[1]]),
  immuneOrNot = c(
    "immune", "immune", "other",
    "immune", "other", "other",
    "other", "immune", "other",
    "other", "immune", "other",
    "other", "immune", "other",
    "other", "other", "other",
    "immune", "immune", "immune",
    "other", "other", "other",
    "immune", "other", "other",
    "other", "other", "other",
    "other", "other", "other",
    "other", "other", "other",
    "immune", "immune", "other",
    "immune", "immune", "immune",
    "other", "other", "other",
    "other", "other", "other",
    "other", "other", "other",
    "immune", "other", "other",
    "other", "immune", "immune",
    "other", "other", "other",
    "immune"
  )
)
immuneCellTypes <- immuneDF$celltype[immuneDF$immuneOrNot == "immune"]
intDF <- do.call("rbind", stringr::str_split(colnames(corMatVec), "<->"))
intIdx <- apply(
  intDF, 1,
  function(x) (x[1] %in% immuneCellTypes) | (x[2] %in% immuneCellTypes)
)
set.seed(92399)
corMatVec_immune <- corMatVec[,intIdx]
pcaImmune_res <- getPC(corMatVec_immune, ncomp = 2, center = T, scale = F)
tempPCAplot(pcaImmune_res) + expand_limits(x = c(-4, 4))

plotLoadings(pcaImmune_res$loadings, "comp1", nfeat = 20)

Jiadong Mao

2026-10-13

Highlights

Introduction

Data preparation

Run PhiSpace

Multi-sample cell type co-localisation analysis