--- title: "STdeconvolve Vignette" author: "Brendan F. Miller" output: BiocStyle::html_document vignette: > %\VignetteIndexEntry{STdeconvolve Vignette} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, results = 'hide') ``` # Installation ```{r, eval=FALSE} if(!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("STdeconvolve") ``` # Introduction `STdeconvolve` is an unsupervised machine learning approach to deconvolve multi-cellular pixel-resolution spatial transcriptomics datasets in order to recover the putative transcriptomic profiles of cell-types and their proportional representation within spatially resolved pixels without reliance on external single-cell transcriptomics references. ## Deconvolution In this tutorial, we will walk through some of the main functionalities of `STdeconvolve`. ```{r, load_library} library(STdeconvolve) ``` Given a counts matrix from pixel-resolution spatial transcriptomics data where each spatially resolved measurement may represent mixtures from potentially multiple cell-types, `STdeconvolve` infers the putative transcriptomic profiles of cell-types and their proportional representation within each multi-cellular spatially resolved pixel. Such a pixel-resolution spatial transcriptomics dataset of the mouse olfactory bulb is built in and can be loaded. ```{r, data} data(mOB) pos <- mOB$pos ## x and y positions of each pixel cd <- mOB$counts ## matrix of gene counts in each pixel annot <- mOB$annot ## annotated tissue layers assigned to each pixel ``` `STdeconvolve` first feature selects for genes most likely to be relevant for distinguishing between cell-types by looking for highly overdispersed genes across ST pixels. Pixels with too few genes or genes with too few reads can also be removed. ```{r, feature_selection, fig.width=6, fig.height=3} ## remove pixels with too few genes counts <- cleanCounts(counts = cd, min.lib.size = 100, min.reads = 1, min.detected = 1, verbose = TRUE) ## feature select for genes corpus <- restrictCorpus(counts, removeAbove = 1.0, removeBelow = 0.05, alpha = 0.05, plot = TRUE, verbose = TRUE) ``` `STdeconvolve` then applies latent Dirichlet allocation (LDA), a generative statistical model commonly used in natural language processing, to discover `K` latent cell-types. `STdeconvolve` fits a range of LDA models to inform the choice of an optimal `K`. ```{r, fit_ldas, fig.width=6, fig.height=4} ## Note: the input corpus needs to be an integer count matrix of pixels x genes ldas <- fitLDA(t(as.matrix(corpus)), Ks = seq(2, 9, by = 1), perc.rare.thresh = 0.05, plot=TRUE, verbose=TRUE) ``` In this example, we will use the model with the lowest model perplexity. The shaded region indicates where a fitted model for a given K had an `alpha` \> 1. `alpha` is an LDA parameter that is solved for during model fitting and corresponds to the shape parameter of a symmetric Dirichlet distribution. In the model, this Dirichlet distribution describes the cell-type proportions in the pixels. A symmetric Dirichlet with `alpha` \> 1 would lead to more uniform cell-type distributions in the pixels and difficulty identifying distinct cell-types. Instead, we want models with `alpha` \< 1, resulting in sparse distributions where only a few cell-types are represented in a given pixel. The resulting `theta` matrix can be interpreted as the proportion of each deconvolved cell-type across each spatially resolved pixel. The resulting `beta` matrix can be interpreted as the putative gene expression profile for each deconvolved cell-type normalized to a library size of 1. This `beta` matrix can be scaled by a depth factor (ex. 1000) for interpretability. ```{r, model} ## select model with minimum perplexity optLDA <- optimalModel(models = ldas, opt = "min") ## Extract pixel cell-type proportions (theta) and cell-type gene expression ## profiles (beta) for the given dataset. ## We can also remove cell-types from pixels that contribute less than 5% of the ## pixel proportion and scale the deconvolved transcriptional profiles by 1000 results <- getBetaTheta(optLDA, perc.filt = 0.05, betaScale = 1000) deconProp <- results$theta deconGexp <- results$beta ``` ## Visualization We can now visualize the proportion of each deconvolved cell-type across the original spatially resolved pixels. ```{r, visualize_all, fig.width=8, fig.height=4} vizAllTopics(deconProp, pos, groups = annot, group_cols = rainbow(length(levels(annot))), r=0.4) ``` For faster plotting, we can visualize the pixel proportions of a single cell-type separately using `vizTopic()`: ```{r, visualize_one, fig.width=8, fig.height=4} vizTopic(theta = deconProp, pos = pos, topic = "5", plotTitle = "X5", size = 5, stroke = 1, alpha = 0.5, low = "white", high = "red") ``` With deconvolved cell-types in hand, we will now go over two different strategies for annotating the deconvolved cell-types. Recall that `STdeconvolve` does not require a reference to deconvolve cell-types in multi-cellular spatially-resolved pixels. However, we would still like some way to identify the deconvolved cell-types to determine if they may represent known cell-types. In addition to the predicted pixel proportions, `STdeconvolve` also returns predicted transcriptional profiles of the deconvolved cell-types as the `beta` matrix. We can use these transcriptional profiles to compare to known cell-type transcriptional profiles and see if we can annotate them. For demonstration purposes, let's use the 5 annotated tissue layer labels (i.e. "Granular Cell Layer", "Mitral Cell Layer", etc) assigned to each pixel and use these to make transcriptional profiles for each of the annotated tissue layers in the MOB. ```{r, proxy_annotations} # proxy theta for the annotated layers mobProxyTheta <- model.matrix(~ 0 + annot) rownames(mobProxyTheta) <- names(annot) # fix names colnames(mobProxyTheta) <- unlist(lapply(colnames(mobProxyTheta), function(x) { unlist(strsplit(x, "annot"))[2] })) mobProxyGexp <- counts %*% mobProxyTheta ``` ## Annotation Strategy 1: Transcriptional correlations First, we can find the Pearson's correlation between the transcriptional profiles of the deconvolved cell-types and those of a ground truth reference. ```{r, correlation_beta, fig.width=6, fig.height=5} corMtx_beta <- getCorrMtx(# the deconvolved cell-type `beta` (celltypes x genes) m1 = as.matrix(deconGexp), # the reference `beta` (celltypes x genes) m2 = t(as.matrix(mobProxyGexp)), # "b" = comparing beta matrices, "t" for thetas type = "b") ## row and column names need to be characters rownames(corMtx_beta) <- paste0("decon_", seq(nrow(corMtx_beta))) correlationPlot(mat = corMtx_beta, # colLabs (aka x-axis, and rows of matrix) colLabs = "Deconvolved cell-types", # rowLabs (aka y-axis, and columns of matrix) rowLabs = "Ground truth cell-types", title = "Transcriptional correlation", annotation = TRUE) + ## this function returns a `ggplot2` object, so can add additional aesthetics ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0)) ``` Notice that cell-type 1, 4, and 5 correlate the strongest with the Granular cell layer, cell-type 2 with the Olfactory nerve layer, etc. These agree with there predicted spatial proportions in the MOB dataset. We can also confirm this by also computing the correlation between the predicted and ground truth cell-type proportions via comparing the `theta` matrices. ```{r, correlation_theta, fig.width=6, fig.height=5} corMtx_theta <- getCorrMtx(# deconvolved cell-type `theta` (pixels x celltypes) m1 = as.matrix(deconProp), # the reference `theta` (pixels x celltypes) m2 = as.matrix(mobProxyTheta), # "b" = comparing beta matrices, "t" for thetas type = "t") ## row and column names need to be characters rownames(corMtx_theta) <- paste0("decon_", seq(nrow(corMtx_theta))) correlationPlot(mat = corMtx_theta, # colLabs (aka x-axis, and rows of matrix) colLabs = "Deconvolved cell-types", # rowLabs (aka y-axis, and columns of matrix) rowLabs = "Ground truth cell-types", title = "Proportional correlation", annotation = TRUE) + ## this function returns a `ggplot2` object, so can add additional aesthetics ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0)) ``` Finally, we can also pair up each reference cell-type with the deconvolved cell-type that has the highest correlation. ```{r, correlation_pairing, fig.width=6, fig.height=5} ## Order the cell-types rows based on best match (highest correlation) with ## each community. ## Cannot have more rows than columns for this pairing, so transpose pairs <- lsatPairs(t(corMtx_theta)) m <- t(corMtx_theta)[pairs$rowix, pairs$colsix] correlationPlot(mat = t(m), # transpose back # colLabs (aka x-axis, and rows of matrix) colLabs = "Deconvolved cell-types", # rowLabs (aka y-axis, and columns of matrix) rowLabs = "Ground truth cell-types", title = "Transcriptional correlation", annotation = TRUE) + ## this function returns a `ggplot2` object, so can add additional aesthetics ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0)) ``` Note that only the paired deconvolved cell-types remain. Ones that paired less strongly with a given ground truth are dropped after assigning pairs. ## Annotation Strategy 2: GSEA Next, given a list of reference gene sets for different cell types, we can performed gene set enrichment analysis on the deconvolved transcriptional profiles to test for significant enrichment of any known ground truth cell-types. First, let's identify marker genes for each tissue layer based on log2(fold-change) compared to the other tissue layers. This will be our list of gene sets for each tissue layer. ```{r, marker_genes} mobProxyLayerMarkers <- list() ## make the tissue layers the rows and genes the columns gexp <- t(as.matrix(mobProxyGexp)) for (i in seq(length(rownames(gexp)))){ celltype <- i ## log2FC relative to other cell-types ## highly expressed in cell-type of interest highgexp <- names(which(gexp[celltype,] > 10)) ## high log2(fold-change) compared to other deconvolved cell-types and limit ## to the top 200 log2fc <- sort( log2(gexp[celltype,highgexp]/colMeans(gexp[-celltype,highgexp])), decreasing=TRUE)[1:200] ## for gene set of the ground truth cell-type, get the genes ## with log2FC > 1 (so FC > 2 over the mean exp of the other cell-types) markers <- names(log2fc[log2fc > 1]) mobProxyLayerMarkers[[ rownames(gexp)[celltype] ]] <- markers } ``` ```{r, annotations} celltype_annotations <- annotateCellTypesGSEA(beta = results$beta, gset = mobProxyLayerMarkers, qval = 0.05) ``` `annotateCellTypesGSEA` returns a list where the first entry, `$results`, contains a list of matrices that show any reference cell-types that had a significant positive enrichment score in each of the deconvolved cell-types. For example, here are the reference cell-types that were significantly enriched in deconvolved cell-type 2: ```{r, annotation_result_1} celltype_annotations$results$`2` ``` Note that the "5: Olfactory Nerve Layer" is significantly positively enriched in the transcriptional profiles of cell-type 2 whereas "1: Granular Cell Layer" is negatively enriched. `annotateCellTypesGSEA` also contains `$predictions`, which is a named vector of the most significant matched reference cell-type with the highest positive enrichment score for each deconvolved cell-type. Note that if there were no significant reference cell-types with positive enrichment, then the deconvolved cell-type will have no matches. ```{r, annotation_result_2} celltype_annotations$predictions ``` Note how the best matches are closely associated with the transcriptional and pixel proportion correlations. # `SpatialExperiment` inputs Spatial datasets stored as a (`SpatialExperiment`)[https://bioconductor.org/packages/release/bioc/html/SpatialExperiment.html] object, and as an extension, 10X Visium datasets made available through (`TENxVisiumData`)[https://bioconductor.org/packages/release/data/experiment/html/TENxVisiumData.html], can be accessed the following way: ```{r, spatial_experiment_install, eval=FALSE} ## install `SpatialExperiment` and `TENxVisiumData` if not already if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install(c("SpatialExperiment", "TENxVisiumData")) ``` ```{r, spatial_experiment_load, eval=FALSE} library(SpatialExperiment) library(TENxVisiumData) ## load the MouseBrainCoronal SpatialExperiment object from `TENxVisiumData` se <- TENxVisiumData::MouseBrainCoronal() ``` Alternatively one can also create a `SpatialExperiment` object from the downloaded [Visium mouse brain section (coronal) dataset](https://www.10xgenomics.com/resources/datasets/mouse-brain-section-coronal-1-standard-1-1-0) directly. In particular, we are interested in the filtered count matrix and the spatial positions of the barcodes. First, make a directory to store the downloaded files: ```{r, eval=FALSE} f <- "visiumTutorial/" if(!file.exists(f)){ dir.create(f) } ``` Download and unzip the `Feature / barcode matrix (filtered)` and the `Spatial imaging data`: ```{r, eval=FALSE} if(!file.exists(paste0(f, "V1_Adult_Mouse_Brain_filtered_feature_bc_matrix.tar.gz"))){ tar_gz_file <- "http://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Adult_Mouse_Brain/V1_Adult_Mouse_Brain_filtered_feature_bc_matrix.tar.gz" download.file(tar_gz_file, destfile = paste0(f, "V1_Adult_Mouse_Brain_filtered_feature_bc_matrix.tar.gz"), method = "auto") } untar(tarfile = paste0(f, "V1_Adult_Mouse_Brain_filtered_feature_bc_matrix.tar.gz"), exdir = f) if(!file.exists(paste0(f, "V1_Adult_Mouse_Brain_spatial.tar.gz"))){ spatial_imaging_data <- "http://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Adult_Mouse_Brain/V1_Adult_Mouse_Brain_spatial.tar.gz" download.file(spatial_imaging_data, destfile = paste0(f, "V1_Adult_Mouse_Brain_spatial.tar.gz"), method = "auto") } untar(tarfile = paste0(f, "V1_Adult_Mouse_Brain_spatial.tar.gz"), exdir = f) ``` Load the filtered counts and spatial barcode information into a `SpatialExperiment`: ```{r, eval=FALSE} se <- SpatialExperiment::read10xVisium(samples = f, type = "sparse", data = "filtered") ``` From here, the count matrix can be accessed and setup for feature selection in `STdeconvolve` via: ```{r, eval=FALSE} ## this is the genes x barcode sparse count matrix ## make sure that SpatialExperiment is loaded because `assay` isn't an exported ## object into the namespace cd <- assay(se, "counts") ``` "x" and "y" coordinates of the barcodes can be obtained via: ```{r, eval=FALSE} pos <- SpatialExperiment::spatialCoords(se) ## change column names to x and y ## for this dataset, we will visualize barcodes using ## "pxl_col_in_fullres" = "y" coordinates, ## and "pxl_row_in_fullres" = "x" coordinates colnames(pos) <- c("y", "x") ``` Poor genes and barcodes will be removed from the count matrix ```{r, eval=FALSE} counts <- cleanCounts(cd, min.lib.size = 100, min.reads = 10) ``` And then we can feature select for overdispersed genes that are present in less than 100% of the barcodes and more than 5%. We will also use the top 1000 most significant overdispersed genes by default. ```{r, eval=FALSE} corpus <- restrictCorpus(counts, removeAbove=1.0, removeBelow = 0.05, nTopOD = 1000) ``` Now, we can fit an LDA model to the feature selected corpus, with K=15 ```{r, eval=FALSE} ldas <- fitLDA(t(as.matrix(corpus)), Ks = c(15)) ``` Next, select the LDA model of interest and get the **beta** (cell-type transcriptional profiles) and **theta** (cell-type barcode proportions) matrices. ```{r, eval=FALSE} optLDA <- optimalModel(models = ldas, opt = 15) results <- getBetaTheta(optLDA, perc.filt = 0.05, betaScale = 1000) deconProp <- results$theta deconGexp <- results$beta ``` Now, we can visualize the deconvolved cell-type proprotions in each spot. ```{r, fig.height=8, fig.width=7, eval=FALSE} plt <- vizAllTopics(theta = deconProp, pos = pos, r = 45, lwd = 0, showLegend = TRUE, plotTitle = NA) + ggplot2::guides(fill=ggplot2::guide_legend(ncol=2)) + ## outer border ggplot2::geom_rect(data = data.frame(pos), ggplot2::aes(xmin = min(x)-90, xmax = max(x)+90, ymin = min(y)-90, ymax = max(y)+90), fill = NA, color = "black", linetype = "solid", size = 0.5) + ggplot2::theme( plot.background = ggplot2::element_blank() ) + ## remove the pixel "groups", which are color aesthetics for the pixel borders ggplot2::guides(colour = "none") ``` # Additional tutorials and analyses For additional tutorials and commands to reproduce the preprocessing of certain datasets used in the manuscript, check out: # Session information ```{r} sessionInfo() ```