## ----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()