10X单细胞空间联合分析之SPOTlight进阶版

hello，大家好，不知道大家发现了没有，穷和忙是一对双胞胎，总是在一个人身上同时出现，越穷越忙，越忙越觉得穷，这个循环什么时候能跳出去？？有时候想想，自己的奋斗目标很多人唾手可得，有点心酸啊~~~

今天周日，我们分享一个简单的，SPOTlight，但是几日不见，软件更新了，这就是老外敬业的地方，文章发表了，软件还在不断的优化，文章在SPOTlight: seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes,我们来梳理一下新的运用示例。

Spatially resolved gene expression profiles are key to understand tissue organization and function. However, novel spatial transcriptomics (ST) profiling techniques lack single-cell resolution and require a combination with single-cell RNA sequencing (scRNA-seq) information to deconvolute the spatially indexed datasets. Leveraging the strengths of both data types, we developed SPOTlight, a computational tool that enables the integration of ST with scRNA-seq data to infer the location of cell types and states within a complex tissue. SPOTlight is centered around a seeded non-negative matrix factorization (NMF) regression, initialized using cell-type marker genes and non-negative least squares (NNLS) to subsequently deconvolute ST capture locations (spots).

加载

library(ggplot2)
library(SPOTlight)
library(SingleCellExperiment)
library(SpatialExperiment)
library(scater)
library(scran)

SPOTlight is a tool that enables the deconvolution of cell types and cell type proportions present within each capture location comprising mixtures of cells. Originally developed for 10X’s Visium - spatial transcriptomics - technology, it can be used for all technologies returning mixtures of cells.

SPOTlight is based on learning topic profile signatures, by means of an NMFreg model, for each cell type and finding which combination of cell types fits best the spot we want to deconvolute. Find below a graphical abstract visually summarizing the key steps.

载入文件

• ST (sparse) matrix with the expression, raw or normalized, where rows = genes and columns = capture locations.
• Single cell (sparse) matrix with the expression, raw or normalized, where rows = genes and columns = cells.
• Vector indicating the cell identity for each column in the single cell expression matrix.

加载示例数据

####sc
library(ExperimentHub)

# initialize a Hub instance which stores a complete set of recordd
eh <- ExperimentHub()

# retrieve any records that match our keyword(s) of interest
query(eh, "Tabula Muris Senis droplet Kidney")

空间数据

library(TENxVisiumData)
spe <- MouseKidneyCoronal()
# Use symbols instead of Ensembl IDs as feature names
rownames(spe) <- rowData(spe)$symbol
library(TabulaMurisSenisData)
sce <- TabulaMurisSenisDroplet(tissues = "Kidney")$Kidney
# Keep cells from 18m mice
sce <- subset(sce, , age == "18m")
# Keep cells with clear cell type annotations
sce <- subset(sce, , !free_annotation %in% c("nan", "CD45"))

Preprocessing

sce <- logNormCounts(sce)

Variance modelling

We aim to identify highly variable genes that drive biological heterogeneity. By feeding these genes to the model we improve the resolution of the biological structure and reduce the technical noise.

dec <- modelGeneVar(sce)
plot(dec$mean, dec$total, xlab = "Mean log-expression", ylab = "Variance")
curve(metadata(dec)$trend(x), col = "blue", add = TRUE)

图片.png

# Get the top 3000 genes.
hvg <- getTopHVGs(dec, n = 3000)

colLabels(sce) <- colData(sce)$free_annotation

# Get vector indicating which genes are neither ribosomal or mitochondrial
genes <- !grepl(pattern = "^Rp[l|s]|Mt", x = rownames(sce))

# Compute marker genes
mgs <- scoreMarkers(sce, subset.row = genes)

mgs_fil <- lapply(names(mgs), function(i) {
    x <- mgs[[i]]
    # Filter and keep relevant marker genes, those with AUC > 0.8
    x <- x[x$mean.AUC > 0.8, ]
    # Sort the genes from highest to lowest weight
    x <- x[order(x$mean.AUC, decreasing = TRUE), ]
    # Add gene and cluster id to the dataframe
    x$gene <- rownames(x)
    x$cluster <- i
    data.frame(x)
})
mgs_df <- do.call(rbind, mgs_fil)

Cell Downsampling

Next, we randomly select at most 100 cells per cell identity. If a cell type is comprised of <100 cells, all the cells will be used. If we have very biologically different cell identities (B cells vs. T cells vs. Macrophages vs. Epithelial) we can use fewer cells since their transcriptional profiles will be very different. In cases when we have more transcriptionally similar cell identities we need to increase our N to capture the biological heterogeneity between them.

In our experience we have found that for this step it is better to select the cells from each cell type from the same batch if you have a joint dataset from multiple runs. This will ensure that the model removes as much signal from the batch as possible and actually learns the biological signal.

For the purpose of this vignette and to speed up the analysis, we are going to use 20 cells per cell identity:

# split cell indices by identity
idx <- split(seq(ncol(sce)), sce$free_annotation)
# downsample to at most 20 per identity & subset
n_cells <- 2
cs_keep <- lapply(idx, function(i) {
    n <- length(i)
    if (n < n_cells)
        n_cells <- n
    sample(i, n_cells)
})
sce <- sce[, unlist(cs_keep)]

Deconvolution

You are now set to run SPOTlight to deconvolute the spots!

Briefly, here is how it works:

1. NMF is used to factorize a matrix into two lower dimensionality matrices without negative elements. We first have an initial matrix V (SCE count matrix), which is factored into W and H. Unit variance normalization by gene is performed in V and in order to standardize discretized gene expression levels, ‘counts-umi’. Factorization is then carried out using the non-smooth NMF method, implemented in the R package NMF NMF (14). This method is intended to return sparser results during the factorization in W and H, thus promoting cell-type-specific topic profile and reducing overfitting during training. Before running factorization, we initialize each topic, column, of W with the unique marker genes for each cell type with weights. In turn, each topic of H in SPOTlight is initialized with the corresponding belongance of each cell for each topic, 1 or 0. This way, we seed the model with prior information, thus guiding it towards a biologically relevant result. This initialization also aims at reducing variability and improving the consistency between runs.

2. NNLS regression is used to map each capture location’s transcriptome in V’ (SE count matrix) to H’ using W as the basis. We obtain a topic profile distribution over each capture location which we can use to determine its composition.

3. we obtain Q, cell-type specific topic profiles, from H. We select all cells from the same cell type and compute the median of each topic for a consensus cell-type-specific topic signature. We then use NNLS to find the weights of each cell type that best fit H’ minimizing the residuals.

res <- SPOTlight(
    x = sce,
    y = spe,
    groups = sce$free_annotation,
    mgs = mgs_df,
    hvg = hvg,
    weight_id = "mean.AUC",
    group_id = "cluster",
    gene_id = "gene")

可视化

plotTopicProfiles(
    x = mod,
    y = sce$free_annotation,
    facet = FALSE,
    min_prop = 0.01,
    ncol = 1) +
    theme(aspect.ratio = 1)

图片.png

plotTopicProfiles(
    x = mod,
    y = sce$free_annotation,
    facet = TRUE,
    min_prop = 0.01,
    ncol = 6)

图片.png

Lastly we can take a look at which genes the model learned for each topic. Higher values indicate that the gene is more relevant for that topic. In the below table we can see how the top genes for Topic1 are characteristic for B cells (i.e. Cd79a, Cd79b, Ms4a1…).

library(NMF)
sign <- basis(mod)
colnames(sign) <- paste0("Topic", seq_len(ncol(sign)))
head(sign)

Spatial Correlation Matrix

plotCorrelationMatrix(mat)

图片.png

Co-localization

plotInteractions(mat, "heatmap")

图片.png

plotInteractions(mat, "network")

图片.png

Scatterpie

ct <- colnames(mat)
mat[mat < 0.1] <- 0

# Define color palette
# (here we use 'paletteMartin' from the 'colorBlindness' package)
paletteMartin <- c(
    "#000000", "#004949", "#009292", "#ff6db6", "#ffb6db", 
    "#490092", "#006ddb", "#b66dff", "#6db6ff", "#b6dbff", 
    "#920000", "#924900", "#db6d00", "#24ff24", "#ffff6d")

pal <- colorRampPalette(paletteMartin)(length(ct))
names(pal) <- ct

plotSpatialScatterpie(
    x = spe,
    y = mat,
    cell_types = colnames(y),
    img = FALSE,
    scatterpie_alpha = 1,
    pie_scale = 0.4) +
    scale_fill_manual(
        values = pal,
        breaks = names(pal))

图片.png

spe$res_ss <- res[[2]][colnames(spe)]
xy <- spatialCoords(spe)
spe$x <- xy[, 1]
spe$y <- xy[, 2]
ggcells(spe, aes(x, y, color = res_ss)) +
    geom_point() +
    scale_color_viridis_c() +
    coord_fixed() +
    theme_bw()

图片.png

生活很好，有你更好