## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, results = 'hide')
## ---- eval=FALSE--------------------------------------------------------------
# if(!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
# BiocManager::install("STdeconvolve")
## ---- load_library------------------------------------------------------------
library(STdeconvolve)
## ---- 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
## ---- 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)
## ---- 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)
## ---- 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
## ---- visualize_all, fig.width=8, fig.height=4--------------------------------
vizAllTopics(deconProp, pos,
groups = annot,
group_cols = rainbow(length(levels(annot))),
r=0.4)
## ---- 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")
## ---- 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
## ---- 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))
## ---- 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))
## ---- 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))
## ---- 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
}
## ---- annotations-------------------------------------------------------------
celltype_annotations <- annotateCellTypesGSEA(beta = results$beta,
gset = mobProxyLayerMarkers,
qval = 0.05)
## ---- annotation_result_1-----------------------------------------------------
celltype_annotations$results$`2`
## ---- annotation_result_2-----------------------------------------------------
celltype_annotations$predictions
## ---- 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"))
## ---- spatial_experiment_load, eval=FALSE-------------------------------------
# library(SpatialExperiment)
# library(TENxVisiumData)
#
# ## load the MouseBrainCoronal SpatialExperiment object from `TENxVisiumData`
# se <- TENxVisiumData::MouseBrainCoronal()
## ---- eval=FALSE--------------------------------------------------------------
# f <- "visiumTutorial/"
#
# if(!file.exists(f)){
# dir.create(f)
# }
## ---- 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)
## ---- eval=FALSE--------------------------------------------------------------
# se <- SpatialExperiment::read10xVisium(samples = f,
# type = "sparse",
# data = "filtered")
## ---- 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")
## ---- 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")
## ---- eval=FALSE--------------------------------------------------------------
# counts <- cleanCounts(cd, min.lib.size = 100, min.reads = 10)
## ---- eval=FALSE--------------------------------------------------------------
# corpus <- restrictCorpus(counts,
# removeAbove=1.0,
# removeBelow = 0.05,
# nTopOD = 1000)
## ---- eval=FALSE--------------------------------------------------------------
# ldas <- fitLDA(t(as.matrix(corpus)), Ks = c(15))
## ---- eval=FALSE--------------------------------------------------------------
# optLDA <- optimalModel(models = ldas, opt = 15)
#
# results <- getBetaTheta(optLDA, perc.filt = 0.05, betaScale = 1000)
# deconProp <- results$theta
# deconGexp <- results$beta
## ---- 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")
## -----------------------------------------------------------------------------
sessionInfo()