Motivation

The spatialHeatmap package provides functionalities for visualizing cell-, tissue- and organ-specific data of biological assays by coloring the corresponding spatial features defined in anatomical images according to quantitative abundance levels of measured biomolecules, such as transcripts, proteins or metabolites (Zhang et al. 2024). A color key is used to represent the quantitative assay values and can be customized by the user. This core functionality of the package is called a spatial heatmap (SHM) plot. Additional important functionalities include Spatial Enrichment (SE), Spatial Co-Expression (SCE), and Single Cell to SHM Co-Visualization (SC2SHM-CoViz). These extra utilities are useful for identifying biomolecules with spatially selective abundance patterns (SE), groups of biomolecules with related abundance profiles using hierarchical clustering, K-means clustering, or network analysis (SCE), and to co-visualize single cell embedding results with SHMs (SCSHM-CoViz). The latter co-visualization functionality (Figure @ref(fig:illus)E) is described in a separate vignette called SCSHM-CoViz.

The functionalities of spatialHeatmap can be used either in a command-driven mode from within R or a graphical user interface (GUI) provided by a Shiny App that is part of this project. While the R-based mode provides flexibility to customize and automate analysis routines, the Shiny App includes a variety of convenience features that will appeal to experimentalists and users less familiar with R. The Shiny App can be used on both local computers as well as centralized server-based deployments (e.g. cloud-based or custom servers). This way it can be used for both hosting community data on a public server or running on a private system. The core functionalities of the spatialHeatmap package are illustrated in Figure @ref(fig:illus).

Overview of spatialHeatmap functionality. (A) The spatialHeatmap package plots numeric assay data onto spatially annotated images. The assay data can be provided as numeric vectors, tabular data, SummarizedExperiment, or SingleCellExperiment objects. The latter two are widely used data containers for organizing both assays as well as associated annotation data and experimental designs. (B) Anatomical and other spatial images need to be provided as annotated SVG (aSVG) files where the spatial features and the corresponding components of the assay data have matching labels (e.g. tissue labels). To work with SVG data efficiently, the SVG S4 class container has been developed. The assay data are used to color the matching spatial features in aSVG images according to a color key. (C)-(D) The result is called a spatial heatmap (SHM). (E) Large-scale data mining such as hierarchical clustering and network analysis can be integrated to facilitate the identification of biomolecules with similar abundance profiles. Moreover, (E) Single cell embedding results can be co-visualized with SHMs.

The package supports anatomical images from public repositories or those provided by users. In general any type of image can be used as long as it can be provided in SVG (Scalable Vector Graphics) format and the corresponding spatial features, such as organs, tissues, cellular compartments, are annotated (see aSVG below). The numeric values plotted onto an SHM are usually quantitative measurements from a wide range of profiling technologies, such as microarrays, next generation sequencing (e.g. RNA-Seq and scRNA-Seq), proteomics, metabolomics, or many other small- or large-scale experiments. For convenience, several preprocessing and normalization methods for the most common use cases are included that support raw and/or preprocessed data. Currently, the main application domains of the spatialHeatmap package are numeric data sets and spatially mapped images from biological, agricultural and biomedical areas. Moreover, the package has been designed to also work with many other spatial data types, such as population data plotted onto geographic maps. This high level of flexibility is one of the unique features of spatialHeatmap. Related software tools for biological applications in this field are largely based on pure web applications (Maag 2018; Lekschas et al. 2015; Papatheodorou et al. 2018; Winter et al. 2007; Waese et al. 2017) or local tools (Muschelli, Sweeney, and Crainiceanu 2014) that typically lack customization functionalities. These restrictions limit users to utilizing pre-existing expression data and/or fixed sets of anatomical image collections. To close this gap for biological use cases, we have developed spatialHeatmap as a generic R/Bioconductor package for plotting quantitative values onto any type of spatially mapped images in a programmable environment and/or in an intuitive to use GUI application.

Design

The core feature of spatialHeatmap is to map assay values (e.g. gene expression data) of one or many biomolecules (e.g. genes) measured under different conditions in form of numerically graded colors onto the corresponding cell types or tissues represented in a chosen SVG image. In the gene profiling field, this feature supports comparisons of the expression values among multiple genes by plotting their SHMs next to each other. Similarly, one can display the expression values of a single or multiple genes across multiple conditions in the same plot (Figure @ref(fig:mul)). This level of flexibility is very efficient for visualizing complicated expression patterns across genes, cell types and conditions. In case of more complex anatomical images with overlapping multiple layer tissues, it is important to visually expose the tissue layer of interest in the plots. To address this, several default and customizable layer viewing options are provided. They allow to hide features in the top layers by making them transparent in order to expose features below them. This transparency viewing feature is highlighted below in the mouse example (Figure @ref(fig:musshm)). Moreover, one can plot multiple distinct aSVGs in a single SHM plot as shown in Figure @ref(fig:arabshm). This is particularly useful for displaying abundance trends across multiple development stages, where each is represented by its own aSVG image. In addition to static SHM representations, one can visualize them in form of interactive HTML files or videos.

To maximize reusability and extensibility, the package organizes large-scale omics assay data along with the associated experimental design information in a SummarizedExperiment object [Figure @ref(fig:illus)A; Morgan et al. (2018)]. The latter is one of the core S4 classes within the Bioconductor ecosystem that has been widely adapted by many other software packages dealing with gene-, protein- and metabolite-level profiling data. In case of gene expression data, the assays slot of the SummarizedExperiment container is populated with a gene expression matrix, where the rows and columns represent the genes and tissue/conditions, respectively. The colData slot contains experimental design definitions including replicate and treatment information. The tissues and/or cell type information in the object maps via colData to the corresponding features in the SVG images using unique identifiers for the spatial features (e.g. tissues or cell types). This allows to color the features of interest in an SVG image according to the numeric data stored in a SummarizedExperiment object. For simplicity the numeric data can also be provided as numeric vectors or data.frames. This can be useful for testing purposes and/or the usage of simple data sets that may not require the more advanced features of the SummarizedExperiment class, such as measurements with only one or a few data points. The details about how to access the SVG images and properly format the associated expression data are provided in the Supplementary Section of this vignette.

Image Format: SVG

SHMs are images where colors encode numeric values in features of any shape. For plotting SHMs, Scalable Vector Graphics (SVG) has been chosen as image format since it is a flexible and widely adapted vector graphics format that provides many advantages for computationally embedding numerical and other information in images. SVG is based on XML formatted text describing all components present in images, including lines, shapes and colors. In case of biological images suitable for SHMs, the shapes often represent anatomical or cell structures. To assign colors to specific features in SHMs, annotated SVG (aSVG) files are used where the shapes of interest are labeled according to certain conventions so that they can be addressed and colored programmatically. One or multiple aSVG files can be parsed and stored in the SVG S4 container with utilities provided by the spatialHeatmap package (Figure @ref(fig:illus)B). The main slots of SVG include coordinate, attribute, dimension, svg, and raster. They correspond to feature coordinates, styling attributes (color, line width, etc.), width and heigth, original aSVG instances, and raster image paths, respectively. Raster images are required only when including photographic image components in SHMs (Figure @ref(fig:overCol)), which is optional. Detailed instruction for creating custom aSVGs is provied in a separate tutorial.

SVGs and aSVGs of anatomical structures can be downloaded from many sources including the repositories described below. Alternatively, users can generate them themselves with vector graphics software such as Inkscape. Typically, in aSVGs one or more shapes of a feature of interest, such as the cell shapes of an organ, are grouped together by a common feature identifier. Via these group identifiers one or many feature types can be colored simultaneously in an aSVG according to biological experiments assaying the corresponding feature types with the required spatial resolution. To color spatial features according to numeric assay values, common identifiers are required for spatial features between the assay data and aSVGs. The color gradient used to visually represent the numeric assay values is controlled by a color gradient parameter. To visually interpret the meaning of the colors, the corresponding color key is included in the SHM plots. Additional details for properly formatting and annotating both aSVG images and assay data are provided in the Supplementary Section section of this vignette.

Data Repositories

If not generated by the user, SHMs can be generated with data downloaded from various public repositories. This includes gene, protein and metabolic profiling data from databases, such as GEO, BAR and Expression Atlas from EMBL-EBI (Papatheodorou et al. 2018). A particularly useful resource, when working with spatialHeatmap, is the EBI Expression Atlas. This online service contains both assay and anatomical images. Its assay data include mRNA and protein profiling experiments for different species, tissues and conditions. The corresponding anatomical image collections are also provided for a wide range of species including animals and plants. In spatialHeatmap several import functions are provided to work with the expression and aSVG repository from the Expression Atlas directly. The aSVG images developed by the spatialHeatmap project are available in its own repository called spatialHeatmap aSVG Repository, where users can contribute their aSVG images that are formatted according to our guidlines.

Tutorial Overview

The following sections of this vignette showcase the most important functionalities of the spatialHeatmap package using as initial example a simple to understand testing data set, and then more complex mRNA profiling data from the Expression Atlas and GEO databases. The co-visualization functionality is explained in a separate vignette (see SCSHM-CoViz). First, SHM plots are generated for both the testing and mRNA expression data. The latter include gene expression data sets from RNA-Seq and microarray experiments of Human Brain, Mouse Organs, Chicken Organs, and Arabidopsis Shoots. The first three are RNA-Seq data from the Expression Atlas, while the last one is a microarray data set from GEO. Second, gene context analysis tools are introduced, which facilitate the visualization of gene modules sharing similar expression patterns. This includes the visualization of hierarchical clustering results with traditional matrix heatmaps (Matrix Heatmap) as well as co-expression network plots (Network). Third, the Spatial Enrichment functionality is illustrated with mouse RNA-seq data. Lastly, an overview of the corresponding Shiny App is presented that provides access to the same functionalities as the R functions, but executes them in an interactive GUI environment (Chang et al. 2021; Chang and Borges Ribeiro 2018).

Installation

The spatialHeatmap package should be installed from an R (version ≥ 3.6) session with the BiocManager::install command.

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("spatialHeatmap")

Packages and Documentation

Next, the packages required for running the sample code in this vignette need to be loaded.

library(spatialHeatmap); library(SummarizedExperiment); library(ExpressionAtlas); library(GEOquery); library(igraph); library(BiocParallel); library(kableExtra); library(org.Hs.eg.db); library(org.Mm.eg.db); library(ggplot2)

The following lists the vignette(s) of this package in an HTML browser. Clicking the corresponding name will open this vignette.

browseVignettes('spatialHeatmap')

To reduce runtime, intermediate results can be cached under ~/.cache/shm.

cache.pa <- '~/.cache/shm' # Path of the cache directory.

A temporary directory is created to save output files.

tmp.dir <- normalizePath(tempdir(check=TRUE), winslash="/", mustWork=FALSE)

Quick Start

Spatial Heatmaps (SHMs) are plotted with the shm function. To provide a quick and intuitive overview how these plots are generated, the following uses a generalized tesing example where a small vector of random numeric values is generated that are used to color features in an aSVG image. The image chosen for this example is an aSVG depicting the human brain. The corresponding image file homo_sapiens.brain.svg is included in this package for testing purposes.

svg.hum.pa <- system.file("extdata/shinyApp/data", 'homo_sapiens.brain.svg', package="spatialHeatmap")
svg.hum <- read_svg(svg.hum.pa)

feature.hum <- attribute(svg.hum)[[1]]
tail(feature.hum[, 1:6], 3) # Partial features and respective attributes
## # A tibble: 3 × 6
##   feature                 id             fill  stroke sub.feature          index
##   <chr>                   <chr>          <chr>  <dbl> <chr>                <int>
## 1 cerebellar.hemisphere   UBERON_0002245 none   0.016 cerebellar.hemisphe…   247
## 2 nucleus.accumbens       UBERON_0001882 none   0.016 nucleus.accumbens      248
## 3 telencephalic.ventricle UBERON_0002285 none   0.016 telencephalic.ventr…   249

coord.df <- coordinate(svg.hum)[[1]]
tail(coord.df, 3) # Partial features and respective coordinates
## # A tibble: 3 × 4
##       x     y feature                 index
##   <dbl> <dbl> <chr>                   <int>
## 1  194.  326. telencephalic.ventricle   249
## 2  194.  326. telencephalic.ventricle   249
## 3  194.  326. telencephalic.ventricle   249

set.seed(20) # To obtain reproducible results, a fixed seed is set.
unique(feature.hum$feature)[1:10]
##  [1] "g4320"                 "locus.ceruleus"        "diencephalon"         
##  [4] "medulla.oblongata"     "middle.temporal.gyrus" "caudate.nucleus"      
##  [7] "middle.frontal.gyrus"  "occipital.lobe"        "parietal.lobe"        
## [10] "pineal.gland"
my_vec <- setNames(sample(1:100, 4), c('substantia.nigra', 'putamen', 'prefrontal.cortex', 'notMapped'))
my_vec
##  substantia.nigra           putamen prefrontal.cortex         notMapped 
##                38                63                 2                29

Plot SHM

Before plotting SHMs, the aSVG instance and numeric data are stored in an SPHM object for the sake of efficient data management and reusability.

dat.quick <- SPHM(svg=svg.hum, bulk=my_vec)

Next, the SHM is plotted with the shm function (Figure @ref(fig:testshm)). Internally, the numbers in my_vec are translated into colors based on the color key and then painted onto the corresponding features in the aSVG. In the given example (Figure @ref(fig:testshm)) only three features (‘substantia.nigra’, ‘putamen’, and ‘prefrontal.cortex’) in the aSVG have matching entries in the data my_vec.

shm.res <- shm(data=dat.quick, ID='testing', ncol=1, sub.title.size=20, legend.nrow=3)

SHM of human brain with testing data. The plots from the left to the right represent the color key for the numeric data, followed by four SHM plots and the legend of the spatial features. The numeric values provided in my_vec are used to color the corresponding features in the SHM plots according to the color key while the legend plot identifies the spatial regions.

The named numeric values in my_vec, that have name matches with the features in the chosen aSVG, are stored in the mapped_feature slot.

# Mapped features
spatialHeatmap::output(shm.res)$mapped_feature
##        ID           feature value    fill                    SVG
## 1 testing  substantia.nigra    38 #FF8700 homo_sapiens.brain.svg
## 2 testing           putamen    63 #FF0000 homo_sapiens.brain.svg
## 3 testing prefrontal.cortex     2 #FFFF00 homo_sapiens.brain.svg

Human Brain

This subsection introduces how to query and download cell- and tissue-specific assay data in the Expression Atlas database using the ExpressionAtlas package (Keays 2019). After the choosen data is downloaded directly into a user's R session, the expression values for selected genes can be plotted onto a chosen aSVG image with or without prior preprocessing steps (e.g. normalization).

all.hum <- read_cache(cache.pa, 'all.hum') # Retrieve data from cache.
if (is.null(all.hum)) { # Save downloaded data to cache if it is not cached.
  all.hum <- searchAtlasExperiments(properties="cerebellum", species="Homo sapiens")
  save_cache(dir=cache.pa, overwrite=TRUE, all.hum)
}

all.hum[2, ]

## DataFrame with 1 row and 4 columns
##      Accession      Species                  Type                  Title
##    <character>  <character>           <character>            <character>
## 1 E-GEOD-75852 Homo sapiens RNA-seq of coding RNA Human iPSC-Derived C..

rse.hum <- read_cache(cache.pa, 'rse.hum') # Read data from cache.
if (is.null(rse.hum)) { # Save downloaded data to cache if it is not cached.
  rse.hum <- getAtlasData('E-GEOD-67196')[[1]][[1]]
  save_cache(dir=cache.pa, overwrite=TRUE, rse.hum)
}

colData(rse.hum)[1:2, c(2, 4)]
## DataFrame with 2 rows and 2 columns
##                organism  organism_part
##             <character>    <character>
## SRR1927019 Homo sapiens     cerebellum
## SRR1927020 Homo sapiens frontal cortex

# Remote aSVG repos.
data(aSVG.remote.repo)
tmp.dir.ebi <- file.path(tmp.dir, 'ebi.zip')
tmp.dir.shm <- file.path(tmp.dir, 'shm.zip')
# Download the remote aSVG repos as zip files.
download.file(aSVG.remote.repo$ebi, tmp.dir.ebi)
download.file(aSVG.remote.repo$shm, tmp.dir.shm)
remote <- list(tmp.dir.ebi, tmp.dir.shm)

tmp.dir <- file.path(tempdir(check=TRUE), 'shm') # Empty directory. 

feature.df.hum <- return_feature(feature=c('frontal cortex', 'cerebellum'), species=c('homo sapiens', 'brain'), dir=tmp.dir, remote=remote) # Query aSVGs
feature.df.hum[1:8, ] # Return first 8 rows for checking
unique(feature.df.hum$SVG) # Return all matching aSVGs

svg.dir <- system.file("extdata/shinyApp/data", package="spatialHeatmap") # Directory of the aSVG collection in spatialHeatmap
feature.df.hum <- return_feature(feature=c('frontal cortex', 'cerebellum'), species=c('homo sapiens', 'brain'), keywords.any=TRUE, return.all=FALSE, dir=svg.dir, remote=NULL)

tail(feature.df.hum[, c('feature', 'stroke', 'SVG')], 3)
## # A tibble: 3 × 3
##   feature             stroke SVG                      
##   <chr>                <dbl> <chr>                    
## 1 cerebral.cortex      0.1   mus_musculus.brain_sp.svg
## 2 cerebellum           0.1   mus_musculus.brain_sp.svg
## 3 somatosensor.cortex  0.100 mus_musculus.brain_sp.svg

hum.tar <- system.file('extdata/shinyApp/data/target_human.txt', package='spatialHeatmap')
target.hum <- read.table(hum.tar, header=TRUE, row.names=1, sep='\t') # Importing
colData(rse.hum) <- DataFrame(target.hum) # Loading to "colData"
colData(rse.hum)[c(1:2, 41:42), 4:5]
## DataFrame with 4 rows and 2 columns
##             organism_part     disease
##               <character> <character>
## SRR1927019     cerebellum         ALS
## SRR1927020 frontal cortex         ALS
## SRR1927059     cerebellum      normal
## SRR1927060 frontal cortex      normal

se.nor.hum <- norm_data(data=rse.hum, norm.fun='ESF', log2.trans=TRUE)

se.aggr.hum <- aggr_rep(data=se.nor.hum, sam.factor='organism_part', con.factor='disease', aggr='mean')
assay(se.aggr.hum)[c(120, 49939, 49977), ]

se.fil.hum <- filter_data(data=se.aggr.hum, sam.factor='organism_part', con.factor='disease', pOA=c(0.01, 5), CV=c(0.3, 100))
se.fil.hum <- cvt_id(db='org.Hs.eg.db', data=se.fil.hum, from.id='ENSEMBL', to.id='SYMBOL')

# Subsetting aSVG features.
svg.hum.sub <- sub_sf(svg.hum, show=c(64:132, 162:163, 164, 190:218))
tail(attribute(svg.hum.sub)[[1]][, 1:6], 3)
## # A tibble: 3 × 6
##   feature    id             fill  stroke sub.feature    index
##   <chr>      <chr>          <chr>  <dbl> <chr>          <int>
## 1 cerebellum UBERON_0002037 none   0.016 cerebellum_2_4   216
## 2 cerebellum UBERON_0002037 none   0.016 cerebellum_2_3   217
## 3 cerebellum UBERON_0002037 none   0.016 cerebellum_2_2   218
# Storing assay data and subsetted aSVG in an 'SPHM' object.
dat.hum <- SPHM(svg=svg.hum.sub, bulk=se.fil.hum)

res.hum <- shm(data=dat.hum, ID=c('OLFM4', 'PRSS22'), lay.shm='con', legend.r=1.5, legend.nrow=3, h=0.6)

shm(data=dat.hum, ID=c('OLFM4', 'PRSS22'), lay.shm='con', legend.r=1.5, legend.nrow=3, h=0.6, aspr=2.3, animation.scale=0.7, bar.width=0.1, bar.value.size=4, out.dir=tmp.dir)

res <- shm(data=dat.hum, ID=c('OLFM4', 'PRSS22'), lay.shm='con', legend.r=1.5, legend.nrow=3, h=0.5, aspr=2.3, animation.scale=0.7, bar.width=0.08, bar.value.size=12)
ggsave(file="./shm_hum.svg", plot=output(res)$spatial_heatmap, width=10, height=8) 

write_svg(input=res, out.dir=tmp.dir)

Rscript -e "spatialHeatmap::plot_meta(svg.path=system.file('extdata/shinyApp/data', 'mus_musculus.brain.svg', package='spatialHeatmap'), bulk=system.file('extdata/shinyApp/data', 'bulk_mouse_cocluster.rds', package='spatialHeatmap'), sam.factor='tissue', aggr='mean', ID=c('AI593442', 'Adora1'), ncol=1, bar.width=0.1, legend.nrow=5, h=0.6)"

Mouse Organs

This section generates an SHM plot for mouse data from the Expression Atlas. The code components are very similar to the previous Human Brain example. For brevity, the corresponding text explaining the code has been reduced to a minimum.

all.mus <- read_cache(cache.pa, 'all.mus') # Retrieve data from cache.
if (is.null(all.mus)) { # Save downloaded data to cache if it is not cached.
  all.mus <- searchAtlasExperiments(properties="kidney", species="Mus musculus")
  save_cache(dir=cache.pa, overwrite=TRUE, all.mus)
}

all.mus[7, ]
rse.mus <- read_cache(cache.pa, 'rse.mus') # Read data from cache.
if (is.null(rse.mus)) { # Save downloaded data to cache if it is not cached.
  rse.mus <- getAtlasData('E-MTAB-2801')[[1]][[1]]
  save_cache(dir=cache.pa, overwrite=TRUE, rse.mus)
}

colData(rse.mus)[1:3, ]
## DataFrame with 3 rows and 4 columns
##           AtlasAssayGroup     organism organism_part      strain
##               <character>  <character>   <character> <character>
## SRR594393              g7 Mus musculus         brain      DBA/2J
## SRR594394             g21 Mus musculus         colon      DBA/2J
## SRR594395             g13 Mus musculus         heart      DBA/2J

feature.df.mus <- return_feature(feature=c('heart', 'kidney'), species=c('Mus musculus'), dir=tmp.dir, remote=remote)

feature.df.mus <- return_feature(feature=c('heart', 'kidney'), species=NULL, dir=svg.dir, remote=NULL) 

unique(feature.df.mus$SVG)
## [1] "gallus_gallus.svg"     "mus_musculus.male.svg"

svg.mus.pa <- system.file("extdata/shinyApp/data", "mus_musculus.male.svg", package="spatialHeatmap")
svg.mus <- read_svg(svg.mus.pa)

mus.tar <- system.file('extdata/shinyApp/data/target_mouse.txt', package='spatialHeatmap')
target.mus <- read.table(mus.tar, header=TRUE, row.names=1, sep='\t') # Importing
colData(rse.mus) <- DataFrame(target.mus) # Loading
target.mus[1:3, ]
##           AtlasAssayGroup     organism organism_part strain
## SRR594393              g7 Mus musculus         brain DBA.2J
## SRR594394             g21 Mus musculus         colon DBA.2J
## SRR594395             g13 Mus musculus         heart DBA.2J

rse.mus <- cvt_id(db='org.Mm.eg.db', data=rse.mus, from.id='ENSEMBL', to.id='SYMBOL', desc=TRUE) # Convert Ensembl ids to UniProt ids.  
se.nor.mus <- norm_data(data=rse.mus, norm.fun='CNF', log2.trans=TRUE) # Normalization
se.aggr.mus <- aggr_rep(data=se.nor.mus, sam.factor='organism_part', con.factor='strain', aggr='mean') # Aggregation of replicates
se.fil.mus <- filter_data(data=se.aggr.mus, sam.factor='organism_part', con.factor='strain', pOA=c(0.01, 5), CV=c(0.6, 100)) # Filtering of genes with low counts and variance 

SHM: Transparency

The pre-processed expression data for gene Scml2 is plotted in form of an SHM. In this case the plot includes expression data for 8 tissues across 3 mouse strains.

dat.mus <- SPHM(svg=svg.mus, bulk=se.fil.mus)
shm(data=dat.mus, ID=c('H19'), legend.width=0.7, legend.text.size=10, sub.title.size=9, ncol=3, ft.trans=c('skeletal muscle'), legend.ncol=2, line.size=0.2, line.color='grey70')

SHM of mouse organs. This is a multiple-layer image where the shapes of the ‘skeletal muscle’ is set transparent to expose ‘lung’ and ‘heart’.

The SHM plots in Figures @ref(fig:musshm) and below demonstrate the usage of the transparency feature via the ft.trans parameter. The corresponding mouse organ aSVG image includes overlapping tissue layers. In this case the skelectal muscle layer partially overlaps with lung and heart tissues. To view lung and heart in Figure @ref(fig:musshm), the skelectal muscle tissue is set transparent with ft.trans=c('skeletal muscle').

To fine control the visual effects in feature rich aSVGs, the line.size and line.color parameters are useful. This way one can adjust the thickness and color of complex structures.

gg <- shm(data=dat.mus, ID=c('H19'), legend.text.size=10, sub.title.size=9, ncol=3, legend.ncol=2, line.size=0.1, line.color='grey70')

SHM of mouse organs. This is a multiple-layer image where the view onto ‘lung’ and ‘heart’ is obstructed by displaying the ‘skeletal muscle’ tissue.

A third example on real data from Expression Atlas is SHMs of time series across chicken organs. Since the procedures are the same with the examples above, this example is illustrated in the Supplementary Section.

Arabidopsis Shoot

This section generates an SHM for Arabidopsis thaliana tissues with gene expression data from the Affymetrix microarray technology. The chosen experiment used ribosome-associated mRNAs from several cell populations of shoots and roots that were exposed to hypoxia stress (Mustroph et al. 2009). In this case the expression data will be downloaded from GEO with utilites from the GEOquery package (Davis and Meltzer 2007). The data preprocessing routines are specific to the Affymetrix technology. The remaining code components for generating SHMs are very similar to the previous examples. For brevity, the text in this section explains mainly the steps that are specific to this data set.

gset <- read_cache(cache.pa, 'gset') # Retrieve data from cache.
if (is.null(gset)) { # Save downloaded data to cache if it is not cached.
  gset <- getGEO("GSE14502", GSEMatrix=TRUE, getGPL=TRUE)[[1]]
  save_cache(dir=cache.pa, overwrite=TRUE, gset)
}
se.sh <- as(gset, "SummarizedExperiment")

rownames(se.sh) <- make.names(rowData(se.sh)[, 'Gene.Symbol'])

colData(se.sh)[60:63, 1:2]
## DataFrame with 4 rows and 2 columns
##                            title geo_accession
##                      <character>   <character>
## GSM362227 shoot_hypoxia_pGL2_r..     GSM362227
## GSM362228 shoot_hypoxia_pGL2_r..     GSM362228
## GSM362229 shoot_control_pRBCS_..     GSM362229
## GSM362230 shoot_control_pRBCS_..     GSM362230

svg.sh.pa <- system.file("extdata/shinyApp/data", "arabidopsis.thaliana_shoot_shm.svg", package="spatialHeatmap")
svg.sh <- read_svg(svg.sh.pa)

sh.tar <- system.file('extdata/shinyApp/data/target_arab.txt', package='spatialHeatmap')
target.sh <- read.table(sh.tar, header=TRUE, row.names=1, sep='\t') # Importing
colData(se.sh) <- DataFrame(target.sh) # Loading
target.sh[60:63, ]
##                           col.name     samples conditions
## shoot_hypoxia_pGL2_rep1  GSM362227  shoot_pGL2    hypoxia
## shoot_hypoxia_pGL2_rep2  GSM362228  shoot_pGL2    hypoxia
## shoot_control_pRBCS_rep1 GSM362229 shoot_pRBCS    control
## shoot_control_pRBCS_rep2 GSM362230 shoot_pRBCS    control

se.aggr.sh <- aggr_rep(data=se.sh, sam.factor='samples', con.factor='conditions', aggr='mean') # Replicate agggregation using mean
se.fil.arab <- filter_data(data=se.aggr.sh, sam.factor='samples', con.factor='conditions', pOA=c(0.03, 6), CV=c(0.30, 100)) # Filtering of genes with low intensities and variance

SHM: Microarray

The expression profile for the HRE2 gene is plotted for the control and the hypoxia treatment across six cell types (Figure @ref(fig:shshm)).

dat.sh <- SPHM(svg=svg.sh, bulk=se.fil.arab)
shm(data=dat.sh, ID=c("HRE2"), legend.ncol=2, legend.text.size=10, legend.key.size=0.02)

SHM of Arabidopsis shoots. The expression profile of the HRE2 gene is plotted for control and hypoxia treatment across six cell types.

Superimposing raster and vector graphics

spatialHeatmap allows to superimpose raster images with vector-based SHMs. This way one can generate SHMs that resemble photographic representations of tissues, organs or entire organisms. For this to work the shapes represented in the vector-graphics need to be an aligned carbon copy of the raster image. Supported file formats for the raster image are JPG/JPEG and PNG, and for the vector image it is SVG. Matching raster and vector graphics are indicated by identical base names in their file names (e.g. imageA.png and imageA.svg). The layout order in SHMs composed of multiple independent images can be controlled by numbering the corresponding file pairs accordingly such as imageA_1.png and imageA_1.svg, imageA_2.png and imageA_2.svg, etc.

In the following example, the required image pairs have been pre-generated from a study on abaxial bundle sheath (ABS) cells in maize leaves (Bezrutczyk et al. 2021). Their file names are labeled 1 and 2 to indicate two developmental stages.

Import paths of first png/svg image pair:

raster.pa1 <- system.file('extdata/shinyApp/data/maize_leaf_shm1.png', package='spatialHeatmap')
svg.pa1 <- system.file('extdata/shinyApp/data/maize_leaf_shm1.svg', package='spatialHeatmap')

Import paths of second png/svg image pair:

raster.pa2 <- system.file('extdata/shinyApp/data/maize_leaf_shm2.png', package='spatialHeatmap')
svg.pa2 <- system.file('extdata/shinyApp/data/maize_leaf_shm2.svg', package='spatialHeatmap')

The two pairs of png/svg images are imported in the SVG container svg.overlay.

svg.overlay <- read_svg(svg.path=c(svg.pa1, svg.pa2), raster.path=c(raster.pa1, raster.pa2))

A slice of attributes in the first aSVG instance is shown.

attribute(svg.overlay)[[1]][1:3, ]
## # A tibble: 3 × 10
##   feature id      fill    stroke sub.feature index element parent    index.all
##   <chr>   <chr>   <chr>    <dbl> <chr>       <int> <chr>   <chr>         <int>
## 1 rect817 rect817 none    0.0843 rect817         1 rect    container         1
## 2 cell1   cell1   #98f0aa 0      path819         2 g       container         2
## 3 cell1   cell1   #98f0aa 0      path821         3 g       container         2
## # ℹ 1 more variable: index.sub <int>

Create random quantitative assay data.

df.ovl <- data.frame(matrix(runif(6, min=0, max=5), nrow=3))
colnames(df.ovl) <- c('cell1', 'cell2') # Assign column names.
rownames(df.ovl) <- paste0('gene', seq_len(3)) # Assign row names 
df.ovl[1:2, ]
##          cell1       cell2
## gene1 1.637970 3.788981771
## gene2 1.850373 0.009639693

To minimize masking of the features in the SHMs by dense regions in the raster images, the alpha.overlay argument allows to adjust the transparency level. In Figure @ref(fig:overCol), the spatial features of interest are superimposed onto the raster image.

dat.over <- SPHM(svg=svg.overlay, bulk=df.ovl)
shm(data=dat.over, charcoal=FALSE, ID=c('gene1'), alpha.overlay=0.5, bar.width=0.09, sub.title.vjust=4, legend.r=0.2)

Superimposing raster images with SHMs (colorful backaground). The expression profiles of gene1 are plotted on ABS cells.

Another option for reducing masking effects is to display the raster image in black and white by setting charcoal=TRUE (Figure @ref(fig:overChar)).

shm(data=dat.over, charcoal=TRUE, ID=c('gene1'), alpha.overlay=0.5, bar.width=0.09, sub.title.vjust=4, legend.r=0.2)

Superimposing raster images with SHMs (black and white background). The expression profiles of gene1 are plotted on ABS cells.

SHMs of Multiple Variables

The SHM plots shown so far are restricted to two variables, here spatial features (e.g. tissues) and treatments. In theory, the complexity of experimental designs scales to any number of variables in spatialHeatmap. This section extends to experiments with three or more variables, such as experiments with spatiotemporal resolution and geographical locations, genotypes, treatments, etc.

To visualize multi-variable assay data, the variables are reduced to two by keeping the spatial feature unchanged and combining all other variables into a composite one. For instance, the following example contains four variables including spatial features, time points, drug treatments and injury conditions. The latter three are combined to produce a composite variable.

se.mus.vars <- readRDS(system.file('extdata/shinyApp/data/mus_brain_vars_se.rds', package='spatialHeatmap'))

colData(se.mus.vars)[1:3, ]
## DataFrame with 3 rows and 5 columns
##                                  tissue        time   treatment      injury
##                             <character> <character> <character> <character>
## hippocampus__3DPI.Veh.NoTBI hippocampus        3DPI         Veh       NoTBI
## hippocampus__3DPI.Veh.NoTBI hippocampus        3DPI         Veh       NoTBI
## hippocampus__3DPI.Veh.TBI   hippocampus        3DPI         Veh         TBI
##                                     comVar
##                                <character>
## hippocampus__3DPI.Veh.NoTBI 3DPI.Veh.NoTBI
## hippocampus__3DPI.Veh.NoTBI 3DPI.Veh.NoTBI
## hippocampus__3DPI.Veh.TBI     3DPI.Veh.TBI
unique(colData(se.mus.vars)$comVar)
## [1] "3DPI.Veh.NoTBI"   "3DPI.Veh.TBI"     "3DPI.Drug.NoTBI"  "3DPI.Drug.TBI"   
## [5] "29DPI.Veh.NoTBI"  "29DPI.Veh.TBI"    "29DPI.Drug.NoTBI" "29DPI.Drug.TBI"

se.mus.vars.nor <- norm_data(data=se.mus.vars, norm.fun='ESF', log2.trans=TRUE) # Normalization.
## Normalising: ESF 
##    type 
## "ratio"
se.mus.vars.aggr <- aggr_rep(data=se.mus.vars.nor, sam.factor='tissue', con.factor='comVar', aggr='mean') # Aggregate replicates.
assay(se.mus.vars.aggr)[1:3, 1:3]
##         hippocampus__3DPI.Veh.NoTBI hippocampus__3DPI.Veh.TBI
## Zscan10                    2.815413                0.09246384
## Zbp1                       1.757163                1.87394752
## Xlr4a                      1.412010                0.98498887
##         hippocampus__3DPI.Drug.NoTBI
## Zscan10                     1.089357
## Zbp1                        1.434792
## Xlr4a                       1.460883

svg.mus.brain.pa <- system.file("extdata/shinyApp/data", "mus_musculus.brain.svg", package="spatialHeatmap")
svg.mus.brain <- read_svg(svg.mus.brain.pa)

tail(attribute(svg.mus.brain)[[1]], 3)
## # A tibble: 3 × 10
##   feature          id    fill  stroke sub.feature index element parent index.all
##   <chr>            <chr> <chr>  <dbl> <chr>       <int> <chr>   <chr>      <int>
## 1 hypothalamus     UBER… none    0.05 hypothalam…    22 path    LAYER…        15
## 2 nose             UBER… none    0.05 nose           23 path    LAYER…        16
## 3 corpora.quadrig… UBER… none    0.05 corpora.qu…    24 path    LAYER…        17
## # ℹ 1 more variable: index.sub <int>

SHM: Multiple Variables

The expression values of gene Acnat1 in hippocampus and hypothalamus across the composite variable are mapped to matching aSVG features. The output SHM plots of each composite variable are organized in a composite plot in Figure @ref(fig:dimshm).

dat.mul.dim <- SPHM(svg=svg.mus.brain, bulk=se.mus.vars.aggr)
shm(data=dat.mul.dim, ID=c('Acnat1'), legend.r=1.5, legend.key.size=0.02, legend.text.size=12, legend.nrow=3)

SHM plots of multiple variable. Gene expression values of Acnat1 in hippocampus and hypothalamus under composite variables are mapped to corresponding aSVG features.

Multiple aSVGs

In a spatiotemporal application, different development stages may need to be represented in separate aSVG images. In such a case, the shm function is able to arrange multiple aSVGs in a single SHM plot. To organize the subplots, the names of the separate aSVG files are expected to include the following suffixes: *_shm1.svg, *_shm2.svg, etc.
As a simple testing example, the following stores random numbers as expression values in a data.frame.

df.random <- data.frame(matrix(runif(50, min=0, max=10), nrow=10))
colnames(df.random) <- c('shoot_totalA__treatment1', 'shoot_totalA__treatment2', 'shoot_totalB__treatment1', 'shoot_totalB__treatment2', 'notMapped') # Assign column names
rownames(df.random) <- paste0('gene', 1:10) # Assign row names 
df.random[1:2, ]
##       shoot_totalA__treatment1 shoot_totalA__treatment2
## gene1                 1.924185                2.6496046
## gene2                 4.520996                0.6959067
##       shoot_totalB__treatment1 shoot_totalB__treatment2 notMapped
## gene1                 5.112577                 7.649234  1.257528
## gene2                 4.749874                 4.356737  5.105836

The paths to the aSVG files are obtained, here for younger and older plants using *_shm1 and *_shm1, respectively, which are generated from Mustroph et al. (2009). Subsequently, the two aSVG files are loaded with the read_svg function.

svg.sh1 <- system.file("extdata/shinyApp/data", "arabidopsis.thaliana_organ_shm1.svg", package="spatialHeatmap")
svg.sh2 <- system.file("extdata/shinyApp/data", "arabidopsis.thaliana_organ_shm2.svg", package="spatialHeatmap")
svg.sh.mul <- read_svg(c(svg.sh1, svg.sh2))

The following generates the corresponding SHMs plot for gene2. The orginal image dimensions can be preserved by assigning TRUE to the preserve.scale argument.

dat.mul.svg <- SPHM(svg=svg.sh.mul, bulk=df.random)
shm(data=dat.mul.svg, ID=c('gene2'), width=0.7, legend.r=0.2, legend.width=1, preserve.scale=TRUE, bar.width=0.09, line.color='grey50') 

Note in Figure @ref(fig:arabshm) shoots are drawn with thicker outlines than roots. This is another useful feature of shm, i.e. preserving the outline thicknesses defined in aSVG files. This feature is particularly useful in cellular SHMs where different cell types may have different cell-wall thicknesses. The outline widths can be updated with update_feature programatically, or within Inkscape manually. The former is illustrated in the Supplementary Section.

Spatial Enrichment

The Spatial Enrichment (SpEn) functionality is an extension of SHMs for identifying groups of biomolecules (e.g. RNAs, proteins, metabolites) that are particularly abundant or enriched in certain spatial regions, such as tissue-specific transcripts. Given a group of spatial features, SpEn identifies biomolecules significantly up- or down-regulated in each spatial feature relative to all other features (references). These biomolecules are classified as spatially enriched or depleted respectively. Then by querying a feature in the enrichment results, the corresponding enriched and depleted biomolecules will be returned, and their abundance values are subsequently visualized as enrichment SHMs. Similarly, biomolecules enriched or depleted in one experimental variable relative to reference variables can be detected and visualized as well. SpEn utilizes differential expression (DE) analysis methods to detect enriched or depleted biomolecules, including edgeR (McCarthy et al. 2012), limma (Ritchie et al. 2015), DESeq2 (Love, Huber, and Anders 2014).

The application of SpEn is illustrated with the above mouse organ data. The function sf_var is used to subset spatial features and experiment variables of interest in the assay data. In the following, five features (‘brain’, ‘liver’, ‘lung’, ‘colon’, and ‘kidney’) and three experiment variables (mouse strains ‘DBA.2J’, ‘C57BL.6’, and ‘CD1’) are subsetted. The com.by argument specifies whether the enrichment will be performed on spatial features (ft) or variables (var). In the following, com.by is set ft, so SpEn is performed for spatial features and the variables under each spatial feature are treated as replicates.

sub.mus <- sf_var(data=rse.mus, feature='organism_part', ft.sel=c('brain', 'lung', 'colon', 'kidney', 'liver'), variable='strain', var.sel=c('DBA.2J', 'C57BL.6', 'CD1'), com.by='ft')
colData(sub.mus)[1:3, c('organism_part', 'strain')]
## DataFrame with 3 rows and 2 columns
##                organism_part      strain
##                  <character> <character>
## brain__DBA.2J          brain      DBA.2J
## colon__DBA.2J          colon      DBA.2J
## kidney__DBA.2J        kidney      DBA.2J

The subsetted data are filtered. Details about this step are given under the human brain section.

sub.mus.fil <- filter_data(data=sub.mus, pOA=c(0.5, 15), CV=c(0.8, 100), verbose=FALSE)

The SpEn is implemented in the function spatial_enrich. In the following, the method edgeR is chosen and the count data are internally normalized by the TMM method from edgeR. The enrichment results are selected by log2 fold change (log2.fc) and FDR (fdr). The outliers argument specifies a number of outliers allowed in references.

enr.res <- spatial_enrich(sub.mus.fil, method=c('edgeR'), norm='TMM', log2.fc=1, fdr=0.05, outliers=1)

The overlap of enriched biomolecules (type='up') across spatial features are presented in an UpSet plot (plot='upset'). Assigning matrix or venn to the plot argument will present the overlaps in form of a matrix plot or Venn diagram respectively.

ovl_enrich(enr.res, type='up', plot='upset')

Overlap of enriched biomolecules across spatial features.

The enriched and depleted genes in brain are queried with the function query_enrich. In the type column, ‘up’ and ‘down’ refer to ‘enriched’ and ‘depleted’ respectively, while the total column shows the total reference features excluding outliers.

en.brain <- query_enrich(enr.res, 'brain')
## Done!
up.brain <- subset(rowData(en.brain), type=='up' & total==4)
up.brain[1:2, 1:3] # Enriched.
## DataFrame with 2 rows and 3 columns
##               type     total      method
##        <character> <numeric> <character>
## Kif5c           up         4       edgeR
## Mapk10          up         4       edgeR
dn.brain <- subset(rowData(en.brain), type=='down' & total==4)
dn.brain[1:2, 1:3] # Depleted. 
## DataFrame with 2 rows and 3 columns
##                 type     total      method
##          <character> <numeric> <character>
## Cdhr5           down         4       edgeR
## Slc22a18        down         4       edgeR

One enriched (Kif5c) and one depleted (Cdhr5) gene in mouse brain are chosen to plot SHMs. The resulting SHMs are termed as SHMs of spatially-enriched/depleted biomolecules (enrichment SHMs) respectively.

dat.enrich <- SPHM(svg=svg.mus, bulk=en.brain)
shm(data=dat.enrich, ID=c('Kif5c', 'Cdhr5'), legend.r=1, legend.nrow=7, sub.title.size=15, ncol=3, bar.width=0.09, lay.shm='gene')

## ggplots/grobs: mus_musculus.male.svg ... 
## ggplot: Kif5c, DBA.2J  C57BL.6  CD1  
## ggplot: Cdhr5, DBA.2J  C57BL.6  CD1  
## Legend plot ... 
## CPU cores: 1 
## Converting "ggplot" to "grob" ... 
## Kif5c_DBA.2J_1  Kif5c_C57BL.6_1  Kif5c_CD1_1  Cdhr5_DBA.2J_1  Cdhr5_C57BL.6_1  Cdhr5_CD1_1  
## Converting "ggplot" to "grob" ... 
## 

Enrichment SHMs. Top row: SHMs of spatially-enriched gene. Bottom row: SHMs of spatially-depleted gene.

The expression profiles of the two chosen genes (Figure @ref(fig:enr)) are also presented in line graphs.

graph_line(assay(en.brain)[c('Kif5c', 'Cdhr5'), ], lgd.pos='right')

Line graph of gene expression profiles.

Hierarchical Clustering

SHMs are suitable for comparing assay profiles among small number of biomolecules (e.g. few genes or proteins) across cell types and conditions. To also support analysis routines of larger number of biomolecules, spatialHeatmap integrates functionalities for identifying groups of biomolecules with similar and/or dissimilar assay profiles, and subsequently analyzing the results with data mining methods that scale to larger numbers of biomolecules than SHMs, such as hierarchical clustering and network analysis methods.

To identify similar biomolecules, the submatrix function can be used. It identifies biomolecules sharing similar profiles with one or more query biomolecules of interest. The given example uses correlation coefficients as similarity metric. The p argument allows to restrict the number of similar biomolecules to return based on a percentage cutoff relative to the number of biomolecules in the assay data set (e.g. 1% of the top most similar biomolecules). If several query biomolecules are provided then the function returns the similar genes for each query, while assuring uniqueness among biomolecules in the result.

In a typical scenario, a spatial feature of interest is chosen in the first place, then a query biomolecule is chosen for the chosen feature, such as a tissue-specific gene. The following introduces the large-scale data mining using as sample data the preprocessed gene expression data (se.fil.mus) from the Mouse Organs section. The brain is selected as the query feature and ‘Kif5c’ is selected as the query gene, which is spatially enriched in brain and known to play important roles in motor neurons (Kanai et al. 2000).

sub.mat <- submatrix(data=se.fil.mus, ID='Kif5c', p=0.15)

The result from the previous step is the assay matrix subsetted to the genes sharing similar assay profiles with the query gene ‘Kif5c’.

assay(sub.mat)[1:2, 1:3] # Subsetted assay matrix

##           brain__DBA.2J colon__DBA.2J heart__DBA.2J
## Akap6          8.551867      2.331506      7.042128
## Adcyap1r1      7.962620      1.204208      3.775402

Subsequently, hierarchical clustering is applied to the subsetted assay matrix containing only the genes that share profile similarities with the query gene ‘Kif5c’. The clustering result is displayed as a matrix heatmap where the rows and columns are sorted by the corresponding hierarchical clustering dendrograms (Figure @ref(fig:static)). The position of the query gene (‘Kif5c’) is indicated in the heatmap by black lines. Setting static=FALSE will launch the interactive mode, where users can zoom into the heatmap by selecting subsections in the image or zoom out by double clicking.

res.hc <- matrix_hm(ID=c('Kif5c'), data=sub.mat, angleCol=60, angleRow=60, cexRow=0.8, cexCol=0.8, margin=c(10, 6), static=TRUE, arg.lis1=list(offsetRow=0.01, offsetCol=0.01))

Matrix Heatmap. Rows are genes and columns are samples. The query gene is tagged by black lines.

The most important information of Figure @ref(fig:static) is returned in res.hc. The row dendrogram is saved in res.hc$rowDendrogram. By using the function cut_dendro, it can be cut at a certain height (here h=15) to obtain the cluster containing the query gene.

cut_dendro(res.hc$rowDendrogram, h=15, 'Kif5c')

## [1] "Aplp1"  "Fbxl16" "Kif1a"  "Kif5a"  "Kif5c"  "Nptxr"  "Nrxn1"  "Ntm"   
## [9] "Stx1b"

Network Graphs

adj.mod <- adj_mod(data=sub.mat)

adj.mod[['adj']][1:3, 1:3]

##               Akap6  Adcyap1r1     Erich6
## Akap6     1.0000000 0.49348330 0.10614679
## Adcyap1r1 0.4934833 1.00000000 0.06216853
## Erich6    0.1061468 0.06216853 1.00000000

adj.mod[['mod']][1:3, ] 

##           0  1  2  3
## Akap6     0 30 30 23
## Adcyap1r1 0  0  0 44
## Erich6    7  5  4 30

Module Visualization

Network modules can be visualized with the network function. To plot a module containing a biomolecule (gene) of interest, its ID needs to be provided under the corresponding argument. Typically, this could be one of the biomolecules chosen for the above SHM plots. There are two modes to visualize the selected module: static or interactive. Figure @ref(fig:inter) was generated with static=TRUE. Setting static=FALSE will generate the interactive version. In the network plot shown below the nodes and edges represent genes and their interactions, respectively. The thickness of the edges denotes the adjacency levels, while the size of the nodes indicates the degree of connectivity of each biomolecule in the network. The adjacency and connectivity levels are also indicated by colors (Figure @ref(fig:inter)). The gene of interest assigned under ID is labeled in the plot with the suffix tag: *_target.

network(ID="Kif5c", data=sub.mat, adj.mod=adj.mod, adj.min=0, vertex.label.cex=1.2, vertex.cex=3, edge.cex=3, mar.key=c(3, 10, 1, 10), static=TRUE)

Static network. Node size denotes gene connectivity while edge thickness stands for co-expression similarity.

Setting static=FALSE launches the interactive network. In this mode there is an interactive color bar that denotes the gene connectivity. To modify it, the color labels need to be provided in a comma separated format, e.g.: yellow, orange, red. The latter would indicate that the gene connectivity increases from yellow to red.

If the expression matrix contains gene/protein annotation information in the rowData slot and specified through desc, then it will be shown when moving the cursor over a network node.

network(ID="Kif5c", data=sub.mat, adj.mod=adj.mod, desc='desc', static=FALSE)

Local System

The Shiny App of spatialHeatmap can be launched from an R session with the following function call.

shiny_shm()

The main dashboard panels of the Shiny App are organized as follows:

1. Datasets: options for uploading custom datasets, selecting default datsets, or downloading example datasets.
   
2. Spatial Heatmap: interactive tables of assay data/metadata, plotting SHMs with single or multiple genes, settings to customize the SHMs.
   
3. Spatial Enrichment: overlap plots of enrichment results across spatial features, data table of enriched and depleted biomolecules.
   
4. Data Mining: hierarchical clustering, K-means clustering, network analysis.
   
5. Co-visualization: co-visualizing single-cell and bulk data in embedding plots (PCA, UMAP, or TSNE) and SHMs repectively. This functionality is described in a separate vignette.

A screenshot is shown below depicting SHM plots generated with the spatialHeatmap Shiny App (Figure @ref(fig:shiny)).

Screenshot of spatialHeatmap’s Shiny App.

The assay data along with metatdata are uploaded to the Shiny App as tabular (e.g. in CSV or TSV format) or ‘.rds’ files. The latter is a SummarizedExperiment object saved with the saveRDS function. In addition, the images are uploaded as aSVG files. The filter_data function can be used to export assay data from SummarizedExperiment to tabular files (TSV format). In the following example, the file argument specifies the the output file name, while the sam.factor and con.factor specifies spatial features and experiment variables in the colData slot respectively, which will be retained in the column names of the output tabular file. An example of the output format is shown in Table @ref(tab:humtab). In addition, the desc argument can be optionally used to specify gene annotations in the rowData slot, which will be appended to the output tabular file. More details of assay data formats are provided in the Supplement.

se.fil.arab <- filter_data(data=se.aggr.sh, desc="Target.Description", sam.factor='samples', con.factor='conditions', pOA=c(0.03, 6), CV=c(0.30, 100), file='./filtered_data.txt')

Server Deployment

As most Shiny Apps, spatialHeatmap can be deployed as a centralized web service. A major advantage of a web server deployment is that the functionalities can be accessed remotely by anyone on the internet without the need to use R on the user system. For deployment one can use custom web servers or cloud services, such as AWS, GCP or shinysapps.io. An example web instance for testing spatialHeatmap online is available here.

Custom Shiny App

The spatialHeatmap package also allows users to create custom Shiny App instances using the custom_shiny function. This function provides options to include custom assay and image data, and define default settings (e.g. color schemes). For details users want to consult the help file. To maximize flexibility, the default settings are stored in a yaml file in the App. This makes it easy to refine and optimize default settings simply by changing this yaml file.

Database Backend

To maintain scalability, the Shiny App is designed to work with backend databases that comprise assay datasets, aSVG files, and the pairing between assay data and aSVGs. This allows users to manage large amounts of assay data and aSVGs in a batch. Both SummarizedExperiment and SinleCellExperiment formats are supported data structures to create the backend database. Assay datasets are saved in the same HDF5 database (Fischer, Smith, and Pau 2020). Then the assay data, aSVG files, and the pairing information (stored as a nested list) between assay data and aSVGs are compressed into a ‘.tar’ file as the final database. More details are referred to the help file of write_hdf5.

Numeric Data

The numceric data used to color the features in aSVG images can be provided as three different object types including vector, data.frame and SummerizedExperiment. When working with complex omics-based assay data then the latter provides the most flexibility, and thus should be the preferred container class for managing numeric data in spatialHeatmap. Both data.frame and SummarizedExperiment can hold data from many measured biomolecules, such as many genes or proteins. In contrast to this, the vector class is only suitable for data from single biomolecules. Due to its simplicity this less complex container is often useful for testing or when dealing with simple data sets.

vec <- sample(x=1:100, size=5) # Random numeric values
names(vec) <- c('putamen__variable1', 'putamen__variable2', 'prefrontal.cortex__variable1', 'prefrontal.cortex__variable2', 'notMapped') # Assign unique names to random values
vec

##           putamen__variable1           putamen__variable2 
##                           47                           38 
## prefrontal.cortex__variable1 prefrontal.cortex__variable2 
##                           62                           70 
##                    notMapped 
##                           10

dat.vec <- SPHM(svg=svg.hum, bulk=vec)
shm(data=dat.vec, ID='testing', ncol=1, legend.r=1.2, sub.title.size=14, ft.trans='g4320', legend.nrow=3)

df.test <- data.frame(matrix(sample(x=1:1000, size=100), nrow=20)) # Create numeric data.frame
colnames(df.test) <- names(vec) # Assign column names
rownames(df.test) <- paste0('gene', 1:20) # Assign row names
df.test[1:3, ]

##       putamen__variable1 putamen__variable2 prefrontal.cortex__variable1
## gene1                902                924                          816
## gene2                564                931                          123
## gene3                602                143                          883
##       prefrontal.cortex__variable2 notMapped
## gene1                          116        58
## gene2                          673       571
## gene3                          380       331

dat.df <- SPHM(svg=svg.hum, bulk=df.test)
shm(data=dat.df, ID=c('gene1'), ncol=1, legend.r=1.2, sub.title.size=14, legend.nrow=3)

df.test$ann <- paste0('ann', 1:20)
df.test[1:3, ]

##       putamen__variable1 putamen__variable2 prefrontal.cortex__variable1
## gene1                902                924                          816
## gene2                564                931                          123
## gene3                602                143                          883
##       prefrontal.cortex__variable2 notMapped  ann
## gene1                          116        58 ann1
## gene2                          673       571 ann2
## gene3                          380       331 ann3

spft <- c(rep('putamen', 4), rep('prefrontal.cortex', 4))
vars <- rep(c('variable1', 'variable1', 'variable2', 'variable2'), 2)
target.test <- data.frame(spFeature=spft, variable=vars, row.names=paste0('assay', 1:8))
target.test

##                spFeature  variable
## assay1           putamen variable1
## assay2           putamen variable1
## assay3           putamen variable2
## assay4           putamen variable2
## assay5 prefrontal.cortex variable1
## assay6 prefrontal.cortex variable1
## assay7 prefrontal.cortex variable2
## assay8 prefrontal.cortex variable2

df.se <- data.frame(matrix(sample(x=1:1000, size=160), nrow=20))
rownames(df.se) <- paste0('gene', 1:20)
colnames(df.se) <- row.names(target.test)
df.se[1:3, ]

##       assay1 assay2 assay3 assay4 assay5 assay6 assay7 assay8
## gene1    213    794    297    129    269    149    521    617
## gene2    961    331    673    636    598    980    839     84
## gene3    133    559    290    317    372    177    844    916

se <- SummarizedExperiment(assays=df.se, colData=target.test)
se

## class: SummarizedExperiment 
## dim: 20 8 
## metadata(0):
## assays(1): ''
## rownames(20): gene1 gene2 ... gene19 gene20
## rowData names(0):
## colnames(8): assay1 assay2 ... assay7 assay8
## colData names(2): spFeature variable

rowData(se) <- df.test['ann']

se.aggr <- aggr_rep(data=se, sam.factor='spFeature', con.factor='variable', aggr='mean')
assay(se.aggr)[1:2, ]

##       putamen__variable1 putamen__variable2 prefrontal.cortex__variable1
## gene1              503.5              213.0                          209
## gene2              646.0              654.5                          789
##       prefrontal.cortex__variable2
## gene1                        569.0
## gene2                        461.5

dat.se <- SPHM(svg=svg.hum, bulk=se.aggr)
shm(data=dat.se, ID=c('gene1'), ncol=1, legend.r=1.2, sub.title.size=14, ft.trans=c('g4320'), legend.nrow=3)

aSVGs files

tmp.dir1 <- file.path(tempdir(check=TRUE), 'shm1') 
if (!dir.exists(tmp.dir1)) dir.create(tmp.dir1)
svg.hum.pa <- system.file("extdata/shinyApp/data", 'homo_sapiens.brain.svg', package="spatialHeatmap") 
file.copy(from=svg.hum.pa, to=tmp.dir1, overwrite=TRUE) # Copy "homo_sapiens.brain.svg" file into 'tmp.dir1'

feature.df <- return_feature(feature=c('frontal cortex', 'prefrontal cortex'), species=c('homo sapiens', 'brain'), dir=tmp.dir1, remote=NULL, keywords.any=FALSE)
feature.df

f.new <- c('prefrontalCortex', 'frontalCortex')

s.new <- c('0.05', '0.1') # New strokes.
c.new <- c('red', 'green') # New colors.

feature.df.new <- cbind(featureNew=f.new, strokeNew=s.new, colorNew=c.new, feature.df)
feature.df.new

update_feature(df.new=feature.df.new, dir=tmp.dir1)

Advanced Functionalities

SHM: Re-matching

spatialHeatmap supports re-matching one spatial feature in the assay data to one or multiple counterparts in the aSVG. The re-matching is defined in a named list. In the list, a name slot refers to a spatial feature in assay data and the corresponding list elements represent aSVG features for re-matching.

The example below takes assay data from the Human Brain section, such as svg.hum.sub and se.fil.hum. The spatial feature frontal.cortex in assay data is re-matched to frontal.cortex and prefrontal.cortex in the aSVG through a list (remat.hum). Figure @ref(fig:shmNoMatch) is the SHM before re-matching while Figure @ref(fig:remat) is the SHM after re-matching.

remat.hum <- list(frontal.cortex=c('frontal.cortex', 'prefrontal.cortex'))

dat.no.match <- SPHM(svg=svg.hum.sub, bulk=se.fil.hum)
shm(data=dat.no.match, ID=c('OLFM4'), lay.shm='con', ncol=1, legend.r=0.8, legend.nrow=2, h=0.6)

SHMs before re-matching.

dat.rematch <- SPHM(svg=svg.hum.sub, bulk=se.fil.hum, match=remat.hum)
shm(data=dat.rematch, ID=c('OLFM4'), lay.shm='con', ncol=1, legend.r=0.8, legend.nrow=2, h=0.6)

SHMs after re-matching. The spatial feature frontal.cortex in assay data is re-matched to aSVG features frontal.cortex and prefrontal.cortex.

Re-matching has various applications, such as mapping gene expression profiles from one species to a closely related species in the same anatomical region. It enables leveraging existing data when direct measurements are challenging. Additionally, re-matching can be used to map gene assay profiles of a sub-tissue to the whole tissue, when assaying the entire tissue is difficult. It allows extrapolating and analyzing gene expression patterns of the entire tissue based on data from the sub-tissues.

SHMs of Time Series

This section generates a SHM plot for chicken data from the Expression Atlas. The code components are very similar to the Human Brain example. For brevity, the corresponding text explaining the code has been reduced to a minimum.

all.chk <- read_cache(cache.pa, 'all.chk') # Retrieve data from cache.
if (is.null(all.chk)) { # Save downloaded data to cache if it is not cached.
  all.chk <- searchAtlasExperiments(properties="heart", species="gallus")
  save_cache(dir=cache.pa, overwrite=TRUE, all.chk)
}

all.chk[3, ]
rse.chk <- read_cache(cache.pa, 'rse.chk') # Read data from cache.
if (is.null(rse.chk)) { # Save downloaded data to cache if it is not cached.
  rse.chk <- getAtlasData('E-MTAB-6769')[[1]][[1]]
  save_cache(dir=cache.pa, overwrite=TRUE, rse.chk)
}

colData(rse.chk)[1:3, ]

## DataFrame with 3 rows and 8 columns
##            AtlasAssayGroup      organism         strain           genotype
##                <character>   <character>    <character>        <character>
## ERR2576379              g1 Gallus gallus Red Junglefowl wild type genotype
## ERR2576380              g1 Gallus gallus Red Junglefowl wild type genotype
## ERR2576381              g2 Gallus gallus Red Junglefowl wild type genotype
##            developmental_stage         age         sex organism_part
##                    <character> <character> <character>   <character>
## ERR2576379              embryo      10 day      female         brain
## ERR2576380              embryo      10 day      female         brain
## ERR2576381              embryo      10 day      female    cerebellum

tmp.dir <- file.path(tempdir(check=TRUE), 'shm')
# Query aSVGs.
feature.df <- return_feature(feature=c('heart', 'kidney'), species=c('gallus'), dir=tmp.dir, remote=remote)

feature.df <- return_feature(feature=c('heart', 'kidney'), species=c('gallus'), dir=svg.dir, remote=NULL)
feature.df[1:2, c('feature', 'stroke', 'SVG')] # A slice of the search result.

svg.chk.pa <- system.file("extdata/shinyApp/data", "gallus_gallus.svg", package="spatialHeatmap")
svg.chk <- read_svg(svg.chk.pa)

chk.tar <- system.file('extdata/shinyApp/data/target_chicken.txt', package='spatialHeatmap')
target.chk <- read.table(chk.tar, header=TRUE, row.names=1, sep='\t') # Importing
colData(rse.chk) <- DataFrame(target.chk) # Loading
target.chk[1:3, ]
##            AtlasAssayGroup      organism         strain           genotype
## ERR2576379              g1 Gallus gallus Red Junglefowl wild type genotype
## ERR2576380              g1 Gallus gallus Red Junglefowl wild type genotype
## ERR2576381              g2 Gallus gallus Red Junglefowl wild type genotype
##            developmental_stage   age    sex organism_part
## ERR2576379              embryo day10 female         brain
## ERR2576380              embryo day10 female         brain
## ERR2576381              embryo day10 female    cerebellum

se.nor.chk <- norm_data(data=rse.chk, norm.fun='ESF', log2.trans=TRUE) # Normalization
se.aggr.chk <- aggr_rep(data=se.nor.chk, sam.factor='organism_part', con.factor='age', aggr='mean') # Replicate agggregation using mean 
se.fil.chk <- filter_data(data=se.aggr.chk, sam.factor='organism_part', con.factor='age', pOA=c(0.01, 5), CV=c(0.6, 100)) # Filtering of genes with low counts and varince

SHM: Time Series

The expression profile for gene ENSGALG00000006346 is plotted across nine time points in four organs in form of a composite SHM with 9 panels. Their layout in three columns is controlled with the argument setting ncol=3. In the legend plot, spatial features are labeled by label=TRUE.

dat.chk <- SPHM(svg=svg.chk, bulk=se.fil.chk)
shm(data=dat.chk, ID='ENSGALG00000006346', width=0.9, legend.width=0.9, legend.r=1.5, sub.title.size=9, ncol=3, legend.nrow=3, label=TRUE, verbose=FALSE)

Time course of chicken organs. The SHM shows the expression profile of a single gene across nine time points and four organs.