spatialHeatmap: Visualizing Spatial Assays in Anatomical Images and Network Graphs
Last update: 27 October, 2020
Package
spatialHeatmap 1.0.0
1 Introduction
1.1 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 a
numeric color key. The color scheme used to represent the assay values can be
customized by the user. This core functionality of the package is called a
spatial heatmap (SHM) plot. It is enhanced with visualization tools for groups
of measured items (e.g. gene modules) sharing related abundance profiles, including
matrix heatmaps combined with hierarchical clustering dendrograms and network representations.
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 also part of this package. 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 other users less
familiar with R. Moreover, the Shiny App can be used on both local computers as
well as centralized server-based deployments (e.g. cloud-based or custom
servers) that can be accessed remotely as a public web service for using
spatialHeatmap’s functionalities with community and/or private data. The
functionalities of the spatialHeatmap package are illustrated in Figure
1.
Figure 1: Overview of spatialHeatmap
(A) The saptialHeatmap package plots numeric assay data onto spatially annotated images. A wide range of omics technologies is supported including genomic, transcriptomic, proteomic and metabolomic profiling data. The assay data can be provided as numeric vectors, tabular data, or SummarizedExperiment objects. The latter is a widely used data container for organizing both assay data as well as associated annotation and experimental design data. (B) Anatomical and other spatial images need to be provided as annotated SVG (aSVG) files where the spatial features and the corresponding data components of the assay data have matching labels (e.g. tissue labels). (C) The assay data are used to color the matching spatial features in one or more aSVG images according to a color key. The result is called a spatial heatmap (SHM) plot. Multiple measurements can be visualized in the same plot, such as several factors (e.g. genes, proteins, metabolites), treatment conditions, growth stages and more. (D) Data mining graphics, such as matrix heatmaps and network graphs, are integrated to facilitate the identification of factors with similar assay profiles. The functionalities of spatialHeatmap can be accessed from local computers via the R console or a graphical user interface based on Shiny. In addtion, the latter can be deployed as a web service on custom servers or cloud-based systems.
As anatomical images the package supports both tissue maps from public repositories and custom images provided by the user. In general any type of image can be used as long as it can be provided in SVG (Scalable Vector Graphics) format, where the corresponding spatial features have been defined (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 a 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 (Lekschas et al. 2015; Papatheodorou et al. 2018; Winter et al. 2007; Waese et al. 2017) or local tools (Maag 2018; 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.
1.2 Design
The core feature of spatialHeatmap is to map assay values (e.g.
gene expression data) of one or many items (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 3). 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 4). Moreover, one can plot multiple distinct
aSVGs in a single SHM plot as shown in Figure 7. 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 dynamic animations
embedded in interactive HTML files or generate videos for them.
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 1A). 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
(Morgan et al. 2018). 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, while the colData slot contains sample data including replicate
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.
1.3 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. 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. Correct assignment of image features and assay results is assured by using for both the same feature identifiers. 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.
1.4 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.
1.5 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 toy data set, and then more complex mRNA profiling data from the
Expression Atlas and GEO databases. First, SHM plots are generated for both the toy
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 co-expression network plots
(Network). Third, 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., n.d.; Chang and Borges Ribeiro 2018). Fourth, more advanced features for plotting customized
SHMs are covered using the Human Brain data set as an example.
2 Getting Started
2.1 Installation
The spatialHeatmap package should be installed from an R (version \(\ge\) 3.6)
session with the BiocManager::install command.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("spatialHeatmap")2.2 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)The following lists the vignette(s) of this package in an HTML browser. Clicking the corresponding name will open this vignette.
browseVignettes('spatialHeatmap')3 Spatial Heatmaps (SHMs)
3.1 Toy Example
SHMs are plotted with the spatial_hm function. To provide a quick
and intuitive overview how these plots are generated, the following uses a
generalized toy 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.
The path to this image on a user's system, where spatialHeatmap is
installed, can be obtained with the system.file function.
3.1.1 aSVG Image
The following commands obtain the directory of the aSVG collection and the full path to the chosen target aSVG image on a user’s system, respectively.
svg.dir <- system.file("extdata/shinyApp/example", package="spatialHeatmap")
svg.hum <- system.file("extdata/shinyApp/example", 'homo_sapiens.brain.svg', package="spatialHeatmap")To identify feature labels of interest in annotated aSVG images, the return_feature
function can be used. The following searches the aSVG images stored in dir
for the query terms ‘lobe’ and ‘homo sapiens’ under the feature and species
fields, respectively. The identified matches are returned as a data.frame.
feature.df <- return_feature(feature=c('lobe'), species=c('homo sapiens'), remote=FALSE, dir=svg.dir)## Accessing features...
## arabidopsis_thaliana.organ_shm.svg, arabidopsis_thaliana.organ_shm1.svg, arabidopsis_thaliana.organ_shm2.svg, arabidopsis_thaliana.root.cross_shm.svg, arabidopsis_thaliana.root.ebi_shm.svg, arabidopsis_thaliana.root.roottip_shm.svg, arabidopsis_thaliana.shoot.root_shm.svg, arabidopsis_thaliana.shoot_shm.svg, gallus_gallus.svg, homo_sapiens.brain.svg, mus_musculus.male.svg, us_map_shm.svg,feature.df## feature id SVG parent index index1
## 1 occipital lobe UBERON_0002021 homo_sapiens.brain.svg LAYER_EFO 9 7
## 2 parietal lobe UBERON_0001872 homo_sapiens.brain.svg LAYER_EFO 10 8
## 3 temporal lobe UBERON_0001871 homo_sapiens.brain.svg LAYER_EFO 26 24fnames <- feature.df[, 1]3.1.2 Numeric Data
The following example generates a small numeric toy vector, where the data slot
contains four numbers and its name slot is populated with the three feature
names obtained from the above aSVG image. In addition, a non-matching entry
(here ‘notMapped’) is included for demonstration purposes. Note, the numbers
are mapped to features via matching names among the numeric vector and the aSVG,
respectively. Accordingly, only numbers and features with matching name
counterparts can be colored in the aSVG image. Entries without name matches
are indicated by a message printed to the R console, here “notMapped”. This
behavior can be turned off with verbose=FALSE in the corresponding function
call. In addition, a summary of the numeric assay to feature mappings is stored
in the result data.frame returned by the spatial_hm function (see below).
my_vec <- sample(1:100, length(unique(fnames))+1)
names(my_vec) <- c(unique(fnames), 'notMapped')
my_vec## occipital lobe parietal lobe temporal lobe notMapped
## 69 83 88 893.1.3 Plot SHM
Next, the SHM is plotted with the spatial_hm function (Figure
2). Internally, the numbers in my_vec are translated into
colors based on the color key assigned to the col.com argument, and then
painted onto the corresponding features in the aSVG, where the path to the image
file is defined by svg.path=svg.hum. The remaining arguments used here include:
ID for defining the title of the plot; ncol for setting the column-wise layout
of the plot excluding the feature legend plot on the right; and height for defining
the height of the SHM relative to its width. In the given example
(Figure 2) only three features in my_vec (‘occipital lobe’,
‘parietal lobe’, and ‘temporal lobe’) have matching entries in the corresponding
aSVG.
shm.df <- spatial_hm(svg.path=svg.hum, data=my_vec, ID='toy', ncol=1, height=0.9, width=0.8, sub.title.size=20, legend.nrow=2)## Syntactically valid column names are made!
## Coordinates: homo_sapiens.brain.svg ...
## Enrties not mapped: notMapped
## Grobs: homo_sapiens.brain.svg ...
## SHMs and legend...
Figure 2: SHM of human brain with toy data
The plots from left to right represent: color key, SHM and legend. The colors in the first two plots depict the user provided numeric values, whereas in the legend plot they are used to map the feature labels to the corresponding spatial regions in the image.
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.
# The SHM and mapped features are stored in a list
names(shm.df)## [1] "spatial_heatmap" "mapped_feature"# Mapped features
shm.df[['mapped_feature']]## rowID featureSVG value SVG
## 1 toy occipital.lobe 69 homo_sapiens.brain.svg
## 2 toy parietal.lobe 83 homo_sapiens.brain.svg
## 3 toy temporal.lobe 88 homo_sapiens.brain.svg3.2 Human Brain
This subsection introduces how to find cell- and tissue-specific assay data in
the Expression Atlas database. After choosing a gene expression experiment, the
data is downloaded directly into a user's R session. Subsequently, the
expression values for selected genes can be plotted onto a chosen aSVG image with
or without prior preprocessing steps (e.g. normalization). For querying and
downloading expression data from the Expression Atlas database, functions from
the ExpressionAtlas package are used (Keays 2019).
3.2.1 Gene Expression Data
The following example searches the Expression Atlas for expression data derived from specific tissues and species of interest, here ‘cerebellum’ and ‘Homo sapiens’, respectively.
all.hum <- searchAtlasExperiments(properties="cerebellum", species="Homo sapiens")The search result is stored in a DFrame containing 13
accessions matching the above query. For the following sample code, the
accession
‘E-GEOD-67196’
from Prudencio et al. (2015) has been chosen, which corresponds
to an RNA-Seq profiling experiment of ‘cerebellum’ and ‘frontal cortex’ brain
tissue from patients with amyotrophic lateral sclerosis (ALS). Details about the
corresponding record can be returned as follows.
all.hum[2, ]## DataFrame with 1 row and 4 columns
## Accession Species Type Title
## <character> <character> <character> <character>
## 1 E-MTAB-3358 Homo sapiens RNA-seq of coding RNA RNA-Seq CAGE (Cap An..The getAtlasData function allows to download the chosen RNA-Seq experiment
from the Expression Atlas and import it into a RangedSummarizedExperiment
object of a user's R session.
rse.hum <- getAtlasData('E-GEOD-67196')[[1]][[1]]The design of the downloaded RNA-Seq experiment is described in the colData slot of
rse.hum. The following returns only its first five rows and columns.
colData(rse.hum)[1:5, 1:5]## DataFrame with 5 rows and 5 columns
## AtlasAssayGroup organism individual organism_part
## <character> <character> <character> <character>
## SRR1927019 g1 Homo sapiens individual1 cerebellum
## SRR1927020 g2 Homo sapiens individual1 frontal cortex
## SRR1927021 g1 Homo sapiens individual2 cerebellum
## SRR1927022 g2 Homo sapiens individual2 frontal cortex
## SRR1927023 g1 Homo sapiens individual34 cerebellum
## disease
## <character>
## SRR1927019 amyotrophic lateral ..
## SRR1927020 amyotrophic lateral ..
## SRR1927021 amyotrophic lateral ..
## SRR1927022 amyotrophic lateral ..
## SRR1927023 amyotrophic lateral ..3.2.2 aSVG Image
The following example shows how to download from the above described SVG
repositories an aSVG image that matches the tissues and species
assayed in the gene expression data set downloaded in the previous subsection.
The return_feature function queries the repository for feature- and
species-related keywords, here c('frontal cortex', 'cerebellum') and c('homo sapiens', 'brain'), respectively. To return matching aSVGs, the argument
keywords.any is set to TRUE by default. When return.all=FALSE, only aSVGs
matching the query keywords are returned and saved under dir. Otherwise, all
aSVGs are returned regardless of the keywords. To avoid overwriting of existing
SVG files, it is recommended to start with an empty target directory, here
tmp.dir. To search a local directory for matching aSVG images, the argument
setting remote=FALSE needs to be used, while specifying the path of the
corresponding directory under dir. All or only matching features are returned
if match.only is set to FALSE or TRUE, respectively.
tmp.dir <- paste0(normalizePath(tempdir(check=TRUE), winslash="/", mustWork=FALSE), '/shm') # Create empty directory
feature.df <- return_feature(feature=c('frontal cortex', 'cerebellum'), species=c('homo sapiens', 'brain'), keywords.any=TRUE, return.all=FALSE, dir=tmp.dir, remote=TRUE, match.only=TRUE, desc=FALSE) # Query aSVGs
feature.df[1:8, ] # Return first 8 rows for checking
unique(feature.df$SVG) # Return all matching aSVGsTo build this vignettes according to the R/Bioconductor package requirements, the
following code section uses the aSVG file instance included in the
spatialHeatmap package rather than the downloaded instance from the previous
example.
feature.df <- return_feature(feature=c('frontal cortex', 'cerebellum'), species=c('homo sapiens', 'brain'), keywords.any=TRUE, return.all=FALSE, dir=svg.dir, remote=FALSE)## Accessing features...
## arabidopsis_thaliana.organ_shm.svg, arabidopsis_thaliana.organ_shm1.svg, arabidopsis_thaliana.organ_shm2.svg, arabidopsis_thaliana.root.cross_shm.svg, arabidopsis_thaliana.root.ebi_shm.svg, arabidopsis_thaliana.root.roottip_shm.svg, arabidopsis_thaliana.shoot.root_shm.svg, arabidopsis_thaliana.shoot_shm.svg, gallus_gallus.svg, homo_sapiens.brain.svg, mus_musculus.male.svg, us_map_shm.svg,Note, the target tissues frontal cortex and cerebellum are included in both
the experimental design slot of the downloaded expression data as well as the
annotations of the aSVG. This way these features can be colored in the downstream
SHM plots. If necessary users can also change from within R the feature identifiers
and names in an aSVG. Details on this utility are provided in the Supplementary Section.
feature.df## feature id SVG parent index
## 1 middle frontal gyrus UBERON_0002702 homo_sapiens.brain.svg LAYER_EFO 8
## 2 cingulate cortex UBERON_0003027 homo_sapiens.brain.svg LAYER_EFO 21
## 3 prefrontal cortex UBERON_0000451 homo_sapiens.brain.svg LAYER_EFO 23
## 4 frontal cortex UBERON_0001870 homo_sapiens.brain.svg LAYER_EFO 24
## 5 cerebral cortex UBERON_0000956 homo_sapiens.brain.svg LAYER_EFO 25
## 6 cerebellum UBERON_0002037 homo_sapiens.brain.svg LAYER_EFO 27
## index1
## 1 6
## 2 19
## 3 21
## 4 22
## 5 23
## 6 25Since the Expression Atlas supports the cross-species anatomy
ontology, the corresponding UBERON identifiers are
included in the id column of the data.frame returned by the above function
call of return_feature (Mungall et al. 2012). This ontology is also supported
by the rols Bioconductor package (Gatto 2019).
3.2.3 Experimental Design
For organizing experimental designs and downstream plotting purposes, it can be
desirable to customize the text in certain columns of colData. This way one can
use the source data for displaying ‘pretty’ sample names in columns and legends
of all downstream tables and plots, respectively, in a consistent and automated
manner. To achieve this, the following example imports a ‘targets’ file that
can be generated and edited by the user in a text or spreadsheet program. In
the following example the target file content is used to replace the text in the
colData slot of the RangedSummarizedExperiment object with a version containing
shorter sample names that are more suitable for plotting purposes.
The following imports a custom target file containing simplified sample labels and experimental design information.
hum.tar <- system.file('extdata/shinyApp/example/target_human.txt', package='spatialHeatmap')
target.hum <- read.table(hum.tar, header=TRUE, row.names=1, sep='\t')Load custom target data into colData slot.
colData(rse.hum) <- DataFrame(target.hum)A slice of the simplified colData object is shown below, where the disease
column contains now shorter labels than in the original data set. Additional
details for generating and using target files in spatialHeatmap are provided
in the Supplementary Section of this vignette.
colData(rse.hum)[c(1:3, 41:42), 4:5]## DataFrame with 5 rows and 2 columns
## organism_part disease
## <character> <character>
## SRR1927019 cerebellum ALS
## SRR1927020 frontal cortex ALS
## SRR1927021 cerebellum ALS
## SRR1927059 cerebellum normal
## SRR1927060 frontal cortex normal3.2.4 Preprocess Assay Data
The actual gene expression data of the downloaded RNA-Seq experiment is stored
in the assay slot of rse.hum. Since it contains raw count data, it can be
desirable to apply basic preprocessing routines prior to plotting spatial
heatmaps. The following shows how to normalize the count data, aggregate
replicates and then remove genes with unreliable expression responses. These
preprocessing steps are optional and can be skipped if needed. For this,
the expression data can be provided to the spatial_hm function directly, where
it is important to assign to the sam.factor and con.factor arguments
the corresponding sample and condition column names (Table 2).
For normalizing raw count data from RNA-Seq experiments, the norm_data
function can be used. It supports the following pre-existing functions from
widely used packages for analyzing count data in the next generation sequencing
(NGS) field: calcNormFactors (CNF) from edgeR (Robinson, McCarthy, and Smyth 2010); as well as
estimateSizeFactors (ESF), varianceStabilizingTransformation (VST), and
rlog from DESeq2 (Love, Huber, and Anders 2014). The argument norm.fun specifies one of the
four internal normalizing methods: CNF, ESF, VST, and rlog. If
norm.fun='none', no normalization is applied. The arguments for each
normalizing function are provided via a parameter.list, which is a list
with named slots. For example, norm.fun='ESF' and
parameter.list=list(type='ratio') is equivalent to
estimateSizeFactors(object, type='ratio'). If paramter.list=NULL, the
default arguments are used by the normalizing function assigned to norm.fun.
For additional details, users want to consult the help file of the norm_data
function by typing ?norm_data in the R console.
The following example uses the ESF normalization option. This method has been
chosen mainly due to its good time performance.
se.nor.hum <- norm_data(data=rse.hum, norm.fun='ESF', data.trans='log2')## Normalising: ESF
## type
## "ratio"Replicates are aggregated with the aggr_rep function, where the summary
statistics can be chosen under the aggr argument (e.g. aggr='mean'). The
columns specifying replicates can be assigned to the sam.factor and
con.factor arguments corresponding to samples and conditions, respectively.
For tracking, the corresponding sample/condition labels are used as column
titles in the aggregated assay instance, where they are concatenated with a
double underscore as separator. In addition, the corresponding rows in the
colData slot are collapsed accordingly.
se.aggr.hum <- aggr_rep(data=se.nor.hum, sam.factor='organism_part', con.factor='disease', aggr='mean')## Syntactically valid column names are made!assay(se.aggr.hum)[1:3, ]## cerebellum__ALS frontal.cortex__ALS cerebellum__normal
## ENSG00000000003 7.024054 7.091484 6.406157
## ENSG00000000005 0.000000 1.540214 0.000000
## ENSG00000000419 7.866582 8.002549 8.073264
## frontal.cortex__normal
## ENSG00000000003 7.004446
## ENSG00000000005 1.403110
## ENSG00000000419 7.955709To remove unreliable expression measures, filtering can be applied.
The following example retains genes with expression values
larger than 5 (log2 space) in at least 1% of all samples (pOA=c(0.01, 5)), and
a coefficient of variance (CV) between 0.30 and 100 (CV=c(0.30, 100)).
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), dir=NULL)## Syntactically valid column names are made!To inspect the results, the following returns three selected rows of the fully preprocessed data matrix (Table 1).
assay(se.fil.hum)[c(5, 733:734), ]| cerebellum__ALS | frontal.cortex__ALS | cerebellum__normal | frontal.cortex__normal | |
|---|---|---|---|---|
| ENSG00000006047 | 1.134172 | 5.2629629 | 0.5377534 | 5.3588310 |
| ENSG00000268433 | 5.324064 | 0.3419665 | 3.4780744 | 0.1340332 |
| ENSG00000268555 | 5.954572 | 2.6148548 | 4.9349736 | 2.0351776 |
3.2.5 SHM: Single Gene
The preprocessed expression values for any gene in the assay slot of
se.fil.hum can be plotted as an SHM. The following uses gene
ENSG00000268433 as an example. The chosen aSVG is a depiction of the human
brain where the assayed featured are colored by the corresponding expression
values in se.fil.hum.
shm.df <- spatial_hm(svg.path=svg.hum, data=se.fil.hum, ID=c('ENSG00000268433'), height=0.7, legend.r=1.5, legend.key.size=0.02, legend.text.size=12, legend.nrow=2)## Coordinates: homo_sapiens.brain.svg ...
## Grobs: homo_sapiens.brain.svg ...
## SHMs and legend...The plotting instructions of the SHM along with the corresponding
mapped features are stored as a list, here named shm.df. Its components
can be accessed as follows.
names(shm.df)## [1] "spatial_heatmap" "mapped_feature"# Mapped features
shm.df[['mapped_feature']]## rowID featureSVG condition value SVG
## 1 ENSG00000268433 cerebellum ALS 5.3240638 homo_sapiens.brain.svg
## 2 ENSG00000268433 frontal.cortex ALS 0.3419665 homo_sapiens.brain.svg
## 3 ENSG00000268433 cerebellum normal 3.4780744 homo_sapiens.brain.svg
## 4 ENSG00000268433 frontal.cortex normal 0.1340332 homo_sapiens.brain.svgIn the above example, the normalized expression values of gene ENSG00000268433
are colored in the frontal cortex and cerebellum, where the different conditions,
here normal and ALS, are given in separate SHMs plotted next to
each other. The color and feature mappings are defined
by the corresponding color key and legend plot on the left and right, respectively.
3.2.6 SHM: Multiple Genes
SHMs for multiple genes can be plotted by providing the
corresponding gene IDs under the ID argument as a character vector. The
spatial_hm function will then sequentially arrange the SHMs for
each gene in a single composite plot. To facilitate comparisons among expression
values across genes and/or conditions, the lay.shm parameter can be assigned
'gene' or 'con', respectively. For instance, in Figure 3 the
SHMs of the genes ENSG00000268433 and ENSG00000006047 are organized
by condition in a horizontal view. This functionality is particularly useful when
comparing gene families. Users can also customize the order of the SHM subplots, by
assigning lay.shm='none'. With this setting the SHM subplots are organized according
to the gene and condition ordering under ID and data, respectively.
spatial_hm(svg.path=svg.hum, data=se.fil.hum, ID=c('ENSG00000268433', 'ENSG00000006047'), lay.shm='con', width=0.8, height=1, legend.r=1.5, legend.nrow=2)## Coordinates: homo_sapiens.brain.svg ...
## Grobs: homo_sapiens.brain.svg ...
## SHMs and legend...
Figure 3: SHMs of two genes
The subplots are organized by “condition” with the lay.shm='con' setting.
3.2.7 SHM: HTML and Video
SHMs can be saved to interactive HTML files as well as video files. To trigger
this export behavior, the argument out.dir needs to be assinged a directory
path where the HTML and video files will be stored. Each HTML file
contains an interactive SHM with zoom in and out functionality. Hovering over
graphics features will display data, gene, condition and other information. The
video will play the SHM subplots in the order specified under the lay.shm
argument.
The following example saves the interactive HTML and video files under
a directory named tmp.dir.
tmp.dir <- paste0(normalizePath(tempdir(check=TRUE), winslash="/"), '/shm')
spatial_hm(svg.path=svg.hum, data=se.fil.hum, ID=c('ENSG00000268433', 'ENSG00000006047'), lay.shm='con', width=0.8, height=1, legend.r=1.5, legend.nrow=2, out.dir=tmp.dir)3.2.8 SHM: Customization
To provide a high level of flexibility, the spatial_hm contains many arguments.
An overview of important arguments and their utility is provided in Table 2.
| argument | description |
|---|---|
| svg.path | Path of aSVG |
| data | Input data of SummarizedExperiment (SE), data frame, or vector |
| sam.factor | Applies to SE. Column name of sample replicates in colData slot. Default is NULL |
| con.factor | Applies to SE. Column name of condition replicates in colData slot. Default is NULL |
| ID | A character vector of row items for plotting spatial heatmaps |
| col.com | A character vector of color components for building colour scale. Default is c(‘yellow’, ‘orange’,‘red’) |
| col.bar | ‘selected’ or ‘all’, the former means use values of ID to build the colour scale while the latter use all values in data. Default is ‘selected’. |
| bar.width | A numeric of colour bar width. Default is 0.7 |
| trans.scale | One of ‘log2’, ‘exp2’, ‘row’, ‘column’, or NULL, which means transform the data by ‘log2’ or ‘2-base expoent’, scale by ‘row’ or ‘column’, or no manipuation respectively. |
| tis.trans | A vector of aSVG features to be transparent. Default is NULL. |
| width, height | Two numerics of overall width and height of all subplots, between 0 and 1 repsectively. The default is 1, 1. |
| legend.r | A numeric to adjust the dimension of the legend plot. Default is 1. The larger, the higher ratio of width to height. |
| sub.title.size | The title size of each spatial heatmap subplot. Default is 11. |
| lay.shm | ‘gen’ or ‘con’, applies to multiple genes or conditions respectively. ‘gen’ means spatial heatmaps are organised by genes while ‘con’ organised by conditions. Default is ‘gen’ |
| ncol | The total column number of spatial heatmaps, not including legend plot. Default is 2. |
| sam.legend | ‘identical’, ‘all’, or a vector of samples/features in aSVG to show in legend plot. ‘identical’ only shows matching features while ‘all’ shows all features. |
| legend.ncol, legend.nrow | Two numbers of columns and rows of legend keys respectively. Default is NULL, NULL, since they are automatically set. |
| legend.position | the position of legend keys (‘none’, ‘left’, ‘right’,‘bottom’, ‘top’), or two-element numeric vector. Default is ‘bottom’. |
| legend.key.size, legend.text.size | The size of legend keys and labels respectively. Default is 0.5 and 8 respectively. |
| line.size, line.color | The size and colour of all plogyon outlines respectively. Default is 0.2 and ‘grey70’ respectively. |
| verbose | TRUE or FALSE. Default is TRUE and the aSVG features not mapped are printed to R console. |
| out.dir | The directory to save HTML and video files of spatial heatmaps. Default is NULL. |
3.3 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.
3.3.1 Gene Expression Data
The chosen mouse RNA-Seq data compares tissue level gene expression across mammalian species (Merkin et al. 2012). The following searches the Expression Atlas for expression data from ‘heart’ and ‘Mus musculus’.
all.mus <- searchAtlasExperiments(properties="heart", species="Mus musculus")Among the many matching entries, accession ‘E-MTAB-2801’ will be downloaded.
all.mus[7, ]## DataFrame with 1 row and 4 columns
## Accession Species Type Title
## <character> <character> <character> <character>
## 1 E-MTAB-2801 Mus musculus RNA-seq of coding RNA Strand-specific RNA-..rse.mus <- getAtlasData('E-MTAB-2801')[[1]][[1]]The design of the downloaded RNA-Seq experiment is described in the colData slot of
rse.mus. The following returns only its first three rows.
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/2J3.3.2 aSVG Image
The following example shows how to download from the above described SVG
repositories an aSVG image that matches the tissues and species
assayed in the gene expression data set downloaded in the previous subsection.
As before the image is saved to a directory named tmp.dir.
tmp.dir <- paste0(normalizePath(tempdir(check=TRUE), winslash="/", mustWork=FALSE), '/shm')
feature.df <- return_feature(feature=c('heart', 'kidney'), species=c('Mus musculus'), keywords.any=TRUE, return.all=FALSE, dir=tmp.dir, remote=TRUE, match.only=FALSE)To build this vignettes according to the R/Bioconductor package requirements, the
following code section uses the aSVG file instance included in the
spatialHeatmap package rather than the downloaded instance from the example in
the previous step.
feature.df <- return_feature(feature=c('heart', 'kidney'), species=NULL, keywords.any=TRUE, return.all=FALSE, dir=svg.dir, remote=FALSE, match.only=FALSE) ## Accessing features...
## arabidopsis_thaliana.organ_shm.svg, arabidopsis_thaliana.organ_shm1.svg, arabidopsis_thaliana.organ_shm2.svg, arabidopsis_thaliana.root.cross_shm.svg, arabidopsis_thaliana.root.ebi_shm.svg, arabidopsis_thaliana.root.roottip_shm.svg, arabidopsis_thaliana.shoot.root_shm.svg, arabidopsis_thaliana.shoot_shm.svg, gallus_gallus.svg, homo_sapiens.brain.svg, mus_musculus.male.svg, us_map_shm.svg,Return the names of the matching aSVG files.
unique(feature.df$SVG)## [1] "gallus_gallus.svg" "mus_musculus.male.svg"The following first selects mus_musculus.male.svg as target aSVG, then
returns the first three rows of the resulting feature.df, and finally prints
the unique set of all aSVG features.
feature.df <- subset(feature.df, SVG=='mus_musculus.male.svg')
feature.df[1:3, ]## feature id SVG parent index index1
## 10 kidney UBERON_0002113 mus_musculus.male.svg LAYER_EFO 14 12
## 11 heart UBERON_0000948 mus_musculus.male.svg LAYER_EFO 51 49
## 12 path4204 path4204 mus_musculus.male.svg LAYER_OUTLINE 1 1unique(feature.df[, 1])## [1] "kidney" "heart"
## [3] "path4204" "aorta"
## [5] "circulatory system" "blood vessel"
## [7] "brown adipose tissue" "white adipose tissue"
## [9] "skin" "stomach"
## [11] "duodenum" "pancreas"
## [13] "spleen" "adrenal gland"
## [15] "colon" "small intestine"
## [17] "caecum" "jejunum"
## [19] "ileum" "esophagus"
## [21] "gall bladder" "parotid gland"
## [23] "submandibular gland" "lymph node"
## [25] "parathyroid gland" "tongue"
## [27] "Peyer’s patch" "prostate gland"
## [29] "vas deferens" "epididymis"
## [31] "testis" "seminal vesicle"
## [33] "penis" "urinary bladder"
## [35] "thymus" "femur"
## [37] "bone marrow" "cartilage"
## [39] "quadriceps femoris" "spinal cord"
## [41] "lung" "diaphragm"
## [43] "peripheral nervous system" "trachea"
## [45] "hindlimb" "trigeminal nerve"
## [47] "eye" "sciatic nerve"
## [49] "intestinal mucosa" "liver"
## [51] "brain" "skeletal muscle"Obtain path of target aSVG on user system.
svg.mus <- system.file("extdata/shinyApp/example", "mus_musculus.male.svg", package="spatialHeatmap")3.3.3 Experimental Design
The following imports a sample target file that is included in this package. To inspect its content, the first three rows of the target file are printed to the screen.
mus.tar <- system.file('extdata/shinyApp/example/target_mouse.txt', package='spatialHeatmap')
target.mus <- read.table(mus.tar, header=TRUE, row.names=1, sep='\t')
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.2Junique(target.mus[, 3])## [1] "brain" "colon" "heart" "kidney"
## [5] "liver" "lung" "skeletal muscle" "spleen"
## [9] "testis"Load custom target data into colData slot.
colData(rse.mus) <- DataFrame(target.mus)3.3.4 Preprocess Assay Data
The raw RNA-Seq count are preprocessed with the following steps: (1) normalization, (2) aggregation of replicates, and (3) filtering of reliable expression data. The details of these steps are explained in the sub-section above using data from human.
se.nor.mus <- norm_data(data=rse.mus, norm.fun='ESF', data.trans='log2') # Normalization## Normalising: ESF
## type
## "ratio"se.aggr.mus <- aggr_rep(data=se.nor.mus, sam.factor='organism_part', con.factor='strain', aggr='mean') # Aggregation of replicates## Syntactically valid column names are made!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), dir=NULL) # Filtering of genes with low counts and variance ## Syntactically valid column names are made!3.3.5 SHM: Transparency
The pre-processed expression data for gene ‘ENSMUSG00000000263’ is plotted in form of an SHM. In this case the plot includes expression data for 8 tissues across 3 mouse strains.
spatial_hm(svg.path=svg.mus, data=se.fil.mus, ID=c('ENSMUSG00000000263'), height=0.7, legend.width=0.7, legend.text.size=10, sub.title.size=9, ncol=3, tis.trans=c('skeletal muscle'), legend.nrow=4, line.size=0.2, line.color='grey70')## Coordinates: mus_musculus.male.svg ...
## Grobs: mus_musculus.male.svg ...
## SHMs and legend...
Figure 4: 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 4 and below demonstrate
the usage of the transparency feature via the tis.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 4, the skelectal
muscle tissue is set transparent with tis.trans=c('skeletal muscle'). To view
in the same aSVG the skeletal muscle tissue instead, tis.trans is assigned
NULL as shown below.
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.
spatial_hm(svg.path=svg.mus, data=se.fil.mus, ID=c('ENSMUSG00000000263'), height=0.6, legend.text.size=10, sub.title.size=9, ncol=3, tis.trans=NULL, legend.ncol=2, line.size=0.2, line.color='grey70')## Coordinates: mus_musculus.male.svg ...
## Grobs: mus_musculus.male.svg ...
## SHMs and legend...3.4 Chicken Organs
This section generates an 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.
3.4.1 Gene Expression Data
The chosen chicken RNA-Seq experiment compares the developmental changes across nine time points of seven organs (Cardoso-Moreira et al. 2019).
The following searches the Expression Atlas for expression data from ‘heart’ and ‘gallus’.
all.chk <- searchAtlasExperiments(properties="heart", species="gallus")Among the matching entries, accession ‘E-MTAB-6769’ will be downloaded.
all.chk[3, ]## DataFrame with 1 row and 4 columns
## Accession Species Type Title
## <character> <character> <character> <character>
## 1 E-MTAB-6769 Gallus gallus RNA-seq of coding RNA Chicken RNA-seq time..rse.chk <- getAtlasData('E-MTAB-6769')[[1]][[1]]The design of the downloaded RNA-Seq experiment is described in the colData
slot of rse.chk. The following returns only its first three rows.
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 cerebellum3.4.2 aSVG Image
The following example shows how to download from the above introduced SVG
repositories an aSVG image that matches the tissues and species
assayed in the gene expression data set downloaded in the previous subsection.
As before the image is saved to a directory named tmp.dir.
tmp.dir <- paste0(normalizePath(tempdir(check=TRUE), winslash="/", mustWork=FALSE), '/shm')
# Query aSVGs.
feature.df <- return_feature(feature=c('heart', 'kidney'), species=c('gallus'), keywords.any=TRUE, return.all=FALSE, dir=tmp.dir, remote=TRUE, match.only=FALSE)To build this vignettes according to the R/Bioconductor package requirements, the
following code section uses the aSVG file instance included in the
spatialHeatmap package rather than the downloaded instance from the previous
step.
feature.df <- return_feature(feature=c('heart', 'kidney'), species=c('gallus'), keywords.any=TRUE, return.all=FALSE, dir=svg.dir, remote=FALSE, match.only=FALSE)## Accessing features...
## arabidopsis_thaliana.organ_shm.svg, arabidopsis_thaliana.organ_shm1.svg, arabidopsis_thaliana.organ_shm2.svg, arabidopsis_thaliana.root.cross_shm.svg, arabidopsis_thaliana.root.ebi_shm.svg, arabidopsis_thaliana.root.roottip_shm.svg, arabidopsis_thaliana.shoot.root_shm.svg, arabidopsis_thaliana.shoot_shm.svg, gallus_gallus.svg, homo_sapiens.brain.svg, mus_musculus.male.svg, us_map_shm.svg,feature.df## feature id SVG parent index
## 1 heart UBERON_0000948 gallus_gallus.svg LAYER_EFO 4
## 2 kidney UBERON_0002113 gallus_gallus.svg LAYER_EFO 5
## 3 chicken_outline chicken_outline gallus_gallus.svg LAYER_OUTLINE 1
## 4 brain UBERON_0000955 gallus_gallus.svg LAYER_EFO 3
## 5 liver UBERON_0002107 gallus_gallus.svg LAYER_EFO 6
## 6 skeletal muscle organ UBERON_0014892 gallus_gallus.svg LAYER_EFO 7
## 7 colon UBERON_0001155 gallus_gallus.svg LAYER_EFO 8
## 8 spleen UBERON_0002106 gallus_gallus.svg LAYER_EFO 9
## 9 lung UBERON_0002048 gallus_gallus.svg LAYER_EFO 10
## index1
## 1 2
## 2 3
## 3 1
## 4 1
## 5 4
## 6 5
## 7 6
## 8 7
## 9 8Obtain path of target aSVG on user system.
svg.chk <- system.file("extdata/shinyApp/example", "gallus_gallus.svg", package="spatialHeatmap")3.4.3 Experimental Design
The following imports a sample target file that is included in this package. To inspect its content, the first three rows of the target file are printed to the screen.
chk.tar <- system.file('extdata/shinyApp/example/target_chicken.txt', package='spatialHeatmap')
target.chk <- read.table(chk.tar, header=TRUE, row.names=1, sep='\t')
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 cerebellumLoad custom target data into colData slot.
colData(rse.chk) <- DataFrame(target.chk)Return samples used for plotting SHMs.
unique(colData(rse.chk)[, 'organism_part'])## [1] "brain" "cerebellum" "heart" "kidney" "ovary"
## [6] "testis" "liver"Return conditions considered for plotting downstream SHM.
unique(colData(rse.chk)[, 'age'])## [1] "day10" "day12" "day14" "day17" "day0" "day155" "day35" "day7"
## [9] "day70"3.4.4 Preprocess Assay Data
The raw RNA-Seq count are preprocessed with the following steps: (1) normalization, (2) aggregation of replicates, and (3) filtering of reliable expression data. The details of these steps are explained in the above sub-section on human data.
se.nor.chk <- norm_data(data=rse.chk, norm.fun='ESF', data.trans='log2') # Normalization## Normalising: ESF
## type
## "ratio"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), dir=NULL) # Filtering of genes with low counts and varince3.4.5 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.
spatial_hm(svg.path=svg.chk, data=se.fil.chk, ID='ENSGALG00000006346', width=0.9, legend.width=0.9, legend.r=1.5, sub.title.size=9, ncol=3, legend.nrow=2, label=TRUE)## Coordinates: gallus_gallus.svg ...
## Enrties not mapped: cerebellum, ovary, testis
## Grobs: gallus_gallus.svg ...
## SHMs and legend...
Figure 5: Time course of chicken organs
The SHM shows the expression profile of a single gene across nine time points and four organs.
3.5 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.
3.5.1 Gene Expression Data
The GSE14502 data set will be downloaded with the getGEO function from the
GEOquery package. Intermediately, the expression data is stored in an
ExpressionSet container (Huber et al. 2015), and then converted to a
SummarizedExperiment object.
gset <- getGEO("GSE14502", GSEMatrix=TRUE, getGPL=TRUE)[[1]]
se.sh <- as(gset, "SummarizedExperiment")The gene symbol identifiers are extracted from the rowData component to be used
as row names. Similarly, one can work with AGI identifiers by providing below AGI
under Gene.Symbol.
rownames(se.sh) <- make.names(rowData(se.sh)[, 'Gene.Symbol'])The following returns a slice of the experimental design stored in the
colData slot. Both the samples and conditions are contained in the title column.
The samples include promoters (pGL2, pCO2, pSCR, pWOL, p35S), tissues
and organs (root atrichoblast epidermis, root cortex meristematic zone, root
endodermis, root vasculature, root_total and shoot_total); and the conditions
are control and hypoxia.
colData(se.sh)[60:63, 1:4]## DataFrame with 4 rows and 4 columns
## title geo_accession status
## <character> <character> <character>
## GSM362227 shoot_hypoxia_pGL2_r.. GSM362227 Public on Oct 12 2009
## GSM362228 shoot_hypoxia_pGL2_r.. GSM362228 Public on Oct 12 2009
## GSM362229 shoot_control_pRBCS_.. GSM362229 Public on Oct 12 2009
## GSM362230 shoot_control_pRBCS_.. GSM362230 Public on Oct 12 2009
## submission_date
## <character>
## GSM362227 Jan 21 2009
## GSM362228 Jan 21 2009
## GSM362229 Jan 21 2009
## GSM362230 Jan 21 20093.5.2 aSVG Image
In this example, the aSVG image has been generated in Inkscape from
the corresponding figure in Mustroph et al. (2009). The resulting custom figure
has been included as a sample aSVG file in the spatialHeatmap package. Detailed
instructions for generating custom aSVG images in Inkscape are provided in the
SVG tutorial.
The annotations in the corresponding aSVG file located under svg.dir can be
queried with the return_features function.
feature.df <- return_feature(feature=c('pGL2', 'pRBCS'), species=c('shoot'), keywords.any=TRUE, return.all=FALSE, dir=svg.dir, remote=FALSE, match.only=FALSE)## Accessing features...
## arabidopsis_thaliana.organ_shm.svg, arabidopsis_thaliana.organ_shm1.svg, arabidopsis_thaliana.organ_shm2.svg, arabidopsis_thaliana.root.cross_shm.svg, arabidopsis_thaliana.root.ebi_shm.svg, arabidopsis_thaliana.root.roottip_shm.svg, arabidopsis_thaliana.shoot.root_shm.svg, arabidopsis_thaliana.shoot_shm.svg, gallus_gallus.svg, homo_sapiens.brain.svg, mus_musculus.male.svg, us_map_shm.svg,The unique set of the matching aSVG files can be returned as follows.
unique(feature.df$SVG)## [1] "arabidopsis_thaliana.shoot.root_shm.svg"
## [2] "arabidopsis_thaliana.shoot_shm.svg"The aSVG file arabidopsis_thaliana.shoot_shm.svg is chosen to generate the SHM in this section.
feature.df <- subset(feature.df, SVG=='arabidopsis_thaliana.shoot_shm.svg')
feature.df[1:3, ]## feature id SVG parent index
## 17 shoot_pGL2 shoot_pGL2 arabidopsis_thaliana.shoot_shm.svg container 2
## 18 shoot_pRBCS shoot_pRBCS arabidopsis_thaliana.shoot_shm.svg container 3
## 19 g258 g258 arabidopsis_thaliana.shoot_shm.svg container 1
## index1
## 17 2
## 18 3
## 19 1Obtain full path of target aSVG on user system.
svg.sh <- system.file("extdata/shinyApp/example", "arabidopsis_thaliana.shoot_shm.svg", package="spatialHeatmap")3.5.3 Experimental Design
The following imports a sample target file that is included in this package. To inspect its content, four selected rows of this target file are printed to the screen.
sh.tar <- system.file('extdata/shinyApp/example/target_arab.txt', package='spatialHeatmap')
target.sh <- read.table(sh.tar, header=TRUE, row.names=1, sep='\t')
target.sh[60:63, ]## col.name sample condition
## 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 controlReturn all samples present in target file.
unique(target.sh[, 'sample'])## [1] "root_total" "root_p35S" "root_pSCR" "root_pSHR"
## [5] "root_pWOL" "root_pGL2" "root_pSUC2" "root_pSultr2.2"
## [9] "root_pCO2" "root_pPEP" "root_pRPL11C" "shoot_total"
## [13] "shoot_p35S" "shoot_pGL2" "shoot_pRBCS" "shoot_pSUC2"
## [17] "shoot_pSultr2.2" "shoot_pCER5" "shoot_pKAT1"Return all conditions present in target file.
unique(target.sh[, 'condition'])## [1] "control" "hypoxia"Load custom target data into colData slot.
colData(se.sh) <- DataFrame(target.sh)3.5.4 Preprocess Assay Data
The downloaded GSE14502 data set has already been normalized with the RMA algorithm (Gautier et al. 2004). Thus, the pre-processing steps can be restricted to the aggregation of replicates and filtering of reliably expressed genes. For the latter, the following code will retain genes with expression values larger than 6 (log2 space) in at least 3% of all samples (pOA=c(0.03, 6)), and with a coefficient of variance (CV) between 0.30 and 100 (CV=c(0.30, 100)).
se.aggr.sh <- aggr_rep(data=se.sh, sam.factor='sample', con.factor='condition', aggr='mean') # Replicate agggregation using mean
se.fil.arab <- filter_data(data=se.aggr.sh, sam.factor='sample', con.factor='condition', pOA=c(0.03, 6), CV=c(0.30, 100), dir=NULL) # Filtering of genes with low intensities and variance3.5.5 SHM: Microarray
The expression profile for the HRE2 gene is plotted for the control and the hypoxia treatment across six cell types (Figure 6).
spatial_hm(svg.path=svg.sh, data=se.fil.arab, ID=c("HRE2"), height=0.7, legend.nrow=3, legend.text.size=11)## Coordinates: arabidopsis_thaliana.shoot_shm.svg ...
## Enrties not mapped: root_total, root_p35S, root_pSCR, root_pSHR, root_pWOL, root_pGL2, root_pSUC2, root_pSultr2.2, root_pCO2, root_pPEP, root_pRPL11C, shoot_total, shoot_p35S
## Grobs: arabidopsis_thaliana.shoot_shm.svg ...
## SHMs and legend...
Figure 6: SHM of Arabidopsis shoots
The expression profile of the HRE2 gene is plotted for control and hypoxia treatment across six cell types.
3.5.6 SHM: Development
The spatial_hm function can arrange multiple aSVGs in a single SHM plot, such
as aSVGs from different development stages. To organize the subplots, the names
of the aSVG files are expected to include the following suffixes: *_shm1.svg,
*_shm2.svg, etc. The paths to the aSVG files are provided under the
svg.path argument. By default, every aSVG image will have a legend plot on
the right. The legend argument provides fine control over which legend plots
to display.
As a simple toy example, the following stores random numbers in a
data.frame.
df.random <- data.frame(matrix(sample(x=1:100, size=50, replace=TRUE), nrow=10))
colnames(df.random) <- c('shoot_totalA__condition1', 'shoot_totalA__condition2', 'shoot_totalB__condition1', 'shoot_totalB__condition2', 'notMapped') # Assign column names
rownames(df.random) <- paste0('gene', 1:10) # Assign row names
df.random[1:3, ]## shoot_totalA__condition1 shoot_totalA__condition2
## gene1 85 13
## gene2 69 71
## gene3 33 75
## shoot_totalB__condition1 shoot_totalB__condition2 notMapped
## gene1 6 55 91
## gene2 8 59 46
## gene3 48 12 9
Next the paths to the aSVG files are obtained, here for younger and older plants using *_shm1
and *_shm1, respectively, which are made from Mustroph et al. (2009).
svg.sh1 <- system.file("extdata/shinyApp/example", "arabidopsis_thaliana.organ_shm1.svg", package="spatialHeatmap")
svg.sh2 <- system.file("extdata/shinyApp/example", "arabidopsis_thaliana.organ_shm2.svg", package="spatialHeatmap")The following generates the corresponding SHMs plot for gene1. The orginal
image dimensions can be preserved by assigning TRUE to the preserve.scale argument.
spatial_hm(svg.path=c(svg.sh1, svg.sh2), data=df.random, ID=c('gene1'), width=0.7, legend.r=0.9, legend.width=1, preserve.scale=TRUE) ## Coordinates: arabidopsis_thaliana.organ_shm1.svg ...
## Coordinates: arabidopsis_thaliana.organ_shm2.svg ...
## Enrties not mapped: shoot_totalB, notMapped
## Enrties not mapped: shoot_totalA, notMapped
## Grobs: arabidopsis_thaliana.organ_shm1.svg ...
## Grobs: arabidopsis_thaliana.organ_shm2.svg ...
## SHMs and legend...
Figure 7: SHMs of Arabidopsis at two growth stages
The expression profile of gene1 under condition1 and condition2 is plotted for two growth stages (top and bottom row).
4 Matrix Heatmaps
SHMs are suitable for comparing assay profiles among small number of items
(e.g. few genes or proteins) across cell types and conditions. To also
support analysis routines of larger number of items, spatialHeatmap integrates
functionalities for identifying groups of items with similar and/or dissimilar
assay profiles, and subsequently analyzing the results with data mining methods
methods that scale to larger numbers of items than SHMs, such as hierarchical
clustering and network analysis methods. The following introduces both
approaches using as sample data the gene expression data from Arabidopsis
shoots from the previous section.
To identify similar items, the submatrix function can be used. The latter
identifies items sharing similar profiles with one or more query items of
interest. The given example uses correlation coefficients as similarity metric.
Three filtering parameters are provided to control the similarity and number of
items in the result. First, the p argument allows to restrict the number of
similar items to return based on a percentage cutoff relative to the number of
items in the assay data set (e.g. 1% of the top most similar items). Second,
a fixed number n of the most similar items can be returned. Third, all items
above a user definable correlation coefficient cutoff value can be returned,
such as a Pearson correlation coefficient (PCC) of at least 0.6. If several
query items are provided then the function returns the similar genes for each
query, while assuring uniqueness among items in the result.
The following example uses as query the two genes: ‘RCA’ and ‘HRE2’. The ann
argument corresponds to the column name in the rowData slot containing gene
annotation information. The latter is only relevant if users want to see these
annotations when mousing over a node in the interactive network
plot generated in the next section.
sub.mat <- submatrix(data=se.fil.arab, ann='Target.Description', ID=c('RCA', 'HRE2'), p=0.1)The result from the previous step is the assay matrix subsetted to the genes sharing similar assay profiles with the query genes ‘RCA’ and ‘HRE2’.
sub.mat[c('RCA', 'HRE2'), c(1:3, 37)] # Subsetted assay matrix## root_total__control root_total__hypoxia root_p35S__control
## RCA 6.569305 6.416811 7.443822
## HRE2 5.486920 11.370161 5.578123
## Target.Description
## RCA hypothetical protein ;supported by full-length cDNA: Ceres:7114.
## HRE2 putative AP2 domain transcription factor
Subsequently, hierarchical clustering is applied to the subsetted assay matrix
containing only the genes that share profile similarities with the query genes
‘RCA’ and ‘HRE2’. The clustering result is displayed as a matrix heatmap where
the rows and columns are sorted by the corresponding hierarchical clustering
dendrograms (Figure 8). The position of the query genes (‘RCA’
and ‘HRE2’) 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.
matrix_hm(ID=c('RCA', 'HRE2'), data=sub.mat, angleCol=80, angleRow=35, cexRow=0.8, cexCol=0.8, margin=c(10, 6), static=TRUE, arg.lis1=list(offsetRow=0.01, offsetCol=0.01))
Figure 8: Matrix Heatmap
Rows are genes and columns are samples. The input genes are tagged by black lines.
5 Network Graphs
5.1 Module Identification
Network analysis is performed with the WGCNA algorithm (Langfelder and Horvath 2008; Ravasz et al. 2002) using as input the subsetted assay matrix generated in the
previous section. The objective is to identify network modules that can be
visualized in the following step in form of network graphs. Applied to the gene
expression sample data used here, these network modules represent groups of
genes sharing highly similar expression profiles. Internally, the network
module identification includes five major steps. First, a correlation matrix is
computed using distance or correlation-based similarity metrics. Second, the
obtained matrix is transformed into an adjacency matrix defining the
connections among items. Third, the adjacency matrix is used to calculate a
topological overlap matrix (TOM) where shared neighborhood information among
items is used to preserve robust connections, while removing spurious
connections. Fourth, the distance transformed TOM is used for hierarchical
clustering. To maximize time performance, the hierarchical clustering is
performed with the flashClust package (Langfelder and Horvath 2012). Fifth, network modules
are identified with the dynamicTreeCut package (Langfelder, Zhang, and Steve Horvath 2016). Its ds
(deepSplit) argument can be assigned integer values from 0 to 3, where
higher values increase the stringency of the module identification process. To
simplify the network module idenification process with WGCNA, the individual
steps can be executed with a single function called adj_mod. The result is a
list containing the adjacency matrix and the final module assignments stored in
a data.frame. Since the interactive network feature used in the
visualization step below performs best on smaller modules, only modules are
returned that were obtained with stringent ds settings (here ds=2 and ds=3).
adj.mod <- adj_mod(data=sub.mat)A slice of the adjacency matrix is shown below.
adj.mod[['adj']][1:3, 1:3]## CA1 PSAH.1 AT2G26500
## CA1 1.0000000 0.9514016 0.9636366
## PSAH.1 0.9514016 1.0000000 0.9611725
## AT2G26500 0.9636366 0.9611725 1.0000000The module assignments are stored in a data frame. Its columns contain the results
for the ds=2 and ds=3 settings. Integer values \(>0\) are the module labels, while \(0\)
indicates unassigned items. The following returns the first three rows of the module
assignment result.
adj.mod[['mod']][1:3, ] ## 2 3
## CA1 1 1
## PSAH.1 1 1
## AT2G26500 1 15.2 Module Visualization
Network modules can be visualized with the network function. To plot a module
containing an item (gene) of interest, its ID (e.g. ‘HRE2’) needs to be
provided under the corresponding argument. Typically, this could be one of the
items chosen for the above SHM plots. There are two modes to visualize the
selected module: static or interactive. Figure 9 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 item in the network. The adjacency and connectivity levels
are also indicated by colors. For example, in Figure 9
connectivity increases from “turquoise” to “violet”, while the adjacency
increases from “yellow” to “blue”. If a gene of interest has been assigned
under ID, then its ID is labeled in the plot with the suffix tag: ’*_target’.
network(ID="HRE2", data=sub.mat, adj.mod=adj.mod, adj.min=0.90, vertex.label.cex=1.2, vertex.cex=2, static=TRUE)
Figure 9: 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 lables 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 subsetted expression matrix contains gene/protein annotation information in the last column, then it will be shown when moving the cursor over a network node.
network(ID="HRE2", data=sub.mat, adj.mod=adj.mod, static=FALSE)6 Shiny App
In additon to running spatialHeatmap from R, the package includes a Shiny
App that provides access to the same
functionalities from an intuitive-to-use web browser interface. Apart from
being very user-friendly, this App conveniently organizes the results of the
entire visualization workflow in a single browser window with options to adjust
the parameters of the individual components interactively. For instance, genes
can be selected and replotted in the SHM simply by clicking the corresponding
rows in the expression table included in the same window.
This representation is very efficient in guiding the interpretation of the results
in a visual and user-friendly manner. For testing purposes, the spatialHeatmap
Shiny App also includes ready-to-use sample expression data and aSVG images
along with embedded user instructions.
6.1 Local System
The Shiny App of spatialHeatmap can be launched from an R session with the following
function call.
shiny_all()The dashboard panels of the Shiny App are organized as follows:
- Left Side Panel: menu for instructions, and selecting input files and parameters
- Data Matrix: scrollable assay matrix uploaded by users
- Spatial Heatmap: spatially colored images based on aSVG images selected and/or uploaded by user
- Matrix Heatmap: hierarchical clustering results of user selected items along with nearest neighbors showing most similar abundance profiles
- Network Graph: interactive network graph for selected module
- Parameter Control Menus: postioned on top of each visualization panel
A screenshot is shown below depicting an SHM plot generated with the Shiny App of spatialHeatmap
(Figure 10).
Figure 10: Screenshot of spatialHeatmap’s Shiny App
After launching, the Shiny App displays by default one of the included data sets.
The assay data and aSVG images are uploaded to the Shiny App as tabular files
(e.g. in CSV or TSV format) and SVG files, respectively. To also allow users
to upload gene expression data stored in SummarizedExperiment objects, one
can export it from R to a tabular file with the filter_data function, where
the user specifies the directory path under the dir argument. This will create
in the target directory a tablular file named customData.txt in TSV format.
The column names in this file preserve the experimental design information from the
colData slot by concatenating the corresponding sample and condition
information separated by double underscores. An example of this format is shown
in Table 1.
se.fil.arab <- filter_data(data=se.aggr.sh, ann="Target.Description", sam.factor='sample', con.factor='condition', pOA=c(0.03, 6), CV=c(0.30, 100), dir='./')To interactively access gene-, transcript- or protein-level annotations in the plots and
tables of the Shiny App, such as viewing functional descriptions by moving the
cursor over network nodes, the corresponding annotation column needs to be
present in the rowData slot and its column name assigned to the ann
argument. In the exported tabular file, the extra annotation column is appended
to the expression matrix.
6.2 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.
6.3 Custom Shiny App
The spatialHeatmap package also allows users to create customized Shiny App
instances using the custom_shiny function. This function provides options to include
custom example data and aSVGs, and define default values within most visualization
panels (e.g. color schemes, image dimensions). For details users want
to consult the help file of the custom_shiny function. To maximize
flexibilty, the default parameters are stored in a yaml file on the Shiny App.
This makes it easy to refine and optimize default parameters simply by changing
this yaml file.
7 Supplementary Section
To generate SHMs with custom data, proper formatting of the numeric (assay) data and aSVG images is required. A tabular target file describing the numeric data can also be provided but is optional. This section provides additional details on these three input types.
7.1 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 items,
such as many genes or proteins. In contrast to this, the
vector class is only suitable for data from single items. Due to its
simplicity this less complex container is often useful for testing or when
dealing with simple data sets.
7.1.1 Object Types
7.1.1.1 vector
When using numeric vectors as input to spatial_hm, then their name slot needs
to be populated with strings matching the feature names in the corresponding aSVG.
To also specify conditions, their labels need to be appended to the feature names
with double underscores as separator, i.e. ’feature__condition’.
The following example replots the toy example for two spatial features (‘occipital lobe’ and ‘parietal lobe’) and two conditions (‘1’ and ‘2’).
vec <- sample(x=1:100, size=5) # Random numeric values
names(vec) <- c('occipital lobe__condition1', 'occipital lobe__condition2', 'parietal lobe__condition1', 'parietal lobe__condition2', 'notMapped') # Assign unique names to random values
vec## occipital lobe__condition1 occipital lobe__condition2
## 8 70
## parietal lobe__condition1 parietal lobe__condition2
## 63 87
## notMapped
## 11With this configuration the resulting plot contains two spatial heatmap plots
for the human brain, one for ‘condition 1’ and another one for ‘contition 2’.
To keep the build time and storage size of this package to a minimum, the
spatial_hm function call in the code block below is not evaluated. Thus,
the corresponding SHM is not shown in this vignette.
spatial_hm(svg.path=svg.hum, data=vec, ID='toy', ncol=1, legend.r=1.2, sub.title.size=14)7.1.1.2 data.frame
Compared to the above vector input, data.frames are structured here like row-wise
appended vectors, where the name slot information in the vectors is stored in the
column names. Each row also contains a name that corresponds to the corresponding
item name, such as a gene ID. The naming of spatial features and conditions in the
column names follows the same conventions as the naming used for the name slots in
the above vector example.
The following illustrates this with an example where a numeric data.frame with
random numbers is generated containing 20 rows and 5 columns. To avoid name clashes,
the values in the rows and columns should be unique.
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, ]## occipital lobe__condition1 occipital lobe__condition2
## gene1 849 809
## gene2 777 370
## gene3 964 47
## parietal lobe__condition1 parietal lobe__condition2 notMapped
## gene1 538 880 257
## gene2 610 421 628
## gene3 992 654 765With the resulting data.frame one can generate the same SHM as in the previous
example, that used a named vector as input to the spatial_hm function. Additionally,
one can now select each row by its name (here gene ID) under the ID argument.
spatial_hm(svg.path=svg.hum, data=df.test, ID=c('gene1'), ncol=1, legend.r=1.2, sub.title.size=14)Additional information can be appended to the data.frame column-wise, such as
annotation data including gene descriptions. This information can then be displayed
interactively in the network plots of the Shiny App by placing the curser over
network nodes.
df.test$ann <- paste0('ann', 1:20)
df.test[1:3, ]## occipital lobe__condition1 occipital lobe__condition2
## gene1 849 809
## gene2 777 370
## gene3 964 47
## parietal lobe__condition1 parietal lobe__condition2 notMapped ann
## gene1 538 880 257 ann1
## gene2 610 421 628 ann2
## gene3 992 654 765 ann37.1.1.3 SummarizedExperiment
The SummarizedExperiment class is a much more extensible and flexible container
for providing metadata for both rows and columns of numeric data stored in tabular
format.
To import experimental design information from tabular files, users can provide
a target file that will be stored in the colData slot of the
SummarizedExperiment (SE, Morgan et al. (2018)) object. In other words, the
target file provides the metadata for the columns of the numeric assay data. Usually,
the target file contains at least two columns: one for the features/samples and
one for the conditions. Replicates are indicated by identical entries in these
columns. The actual numeric matrix representing the assay data is stored in
the assay slot, where the rows correspond to items, such as gene IDs.
Optionally, additional annotation information for the rows (e.g. gene
descriptions) can be stored in the rowData slot.
For constructing a valid SummarizedExperiment object, that can be used by
the spatial_hm function, the target file should meet the following requirements.
It should be imported with
read.tableorread.deliminto adata.frameor thedata.framecan be constructed in R on the fly (as shown below).It should contain at least two columns. One column represents the features/samples and the other one the conditions. The rows in the target file correspond to the columns of the numeric data stored in the
assayslot. If the condition column is empty, then the same condition is assumed under the corresponding features/samples entry.The feature/sample names must have matching entries in to corresponding aSVG to be considered in the final SHM. Note, the double underscore is a special string that is reserved for specific purposes in spatialHeatmap, and thus should be avoided for naming feature/samples and conditions.
The following example illustrates the design of a valid SummarizedExperiment
object for generating SHMs. In this example, the ‘occipital lobe’ tissue has 2
conditions and each condition has 2 replicates. Thus, there are 4 assays for
occipital lobe, and the same design applies to the parietal lobe tissue.
sample <- c(rep('occipital lobe', 4), rep('parietal lobe', 4))
condition <- rep(c('condition1', 'condition1', 'condition2', 'condition2'), 2)
target.test <- data.frame(sample=sample, condition=condition, row.names=paste0('assay', 1:8))
target.test## sample condition
## assay1 occipital lobe condition1
## assay2 occipital lobe condition1
## assay3 occipital lobe condition2
## assay4 occipital lobe condition2
## assay5 parietal lobe condition1
## assay6 parietal lobe condition1
## assay7 parietal lobe condition2
## assay8 parietal lobe condition2The assay slot is populated with a 8 x 20 data.frame containing random
numbers. Each column corresponds to an assay in the target file (here imported
into colData), while each row corresponds to a gene.
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 712 513 145 680 278 659 865 870
## gene2 817 827 252 284 31 988 718 141
## gene3 570 377 310 490 612 652 25 136Next, the final SummarizedExperiment object is constructed by providing the
numeric and target data under the assays and colData arguments,
respectively.
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): sample conditionIf needed row-wise annotation information (e.g. for genes) can be included in
the SummarizedExperiment object as well. This can be done during the
construction of the initial object, or as below by injecting the information
into an existing SummarizedExperiment object.
rowData(se) <- df.test['ann']In this simple example, possible normalization and filtering steps are skipped. Yet, the aggregation of replicates is performed as shown below.
se.aggr <- aggr_rep(data=se, sam.factor='sample', con.factor='condition', aggr='mean')## Syntactically valid column names are made!assay(se.aggr)[1:3, ]## occipital.lobe__condition1 occipital.lobe__condition2
## gene1 612.5 412.5
## gene2 822.0 268.0
## gene3 473.5 400.0
## parietal.lobe__condition1 parietal.lobe__condition2
## gene1 468.5 867.5
## gene2 509.5 429.5
## gene3 632.0 80.5With the fully configured SummarizedExperiment object, a similar SHM is plotted as
in the previous examples.
spatial_hm(svg.path=svg.hum, data=se.aggr, ID=c('gene1'), ncol=1, legend.r=1.2, sub.title.size=14)7.2 aSVG File
7.2.1 aSVG repository
A public aSVG repository, that can be used by spatialHeatmap directly, is the one
maintained by the EBI Gene Expression
Group.
It contains annatomical aSVG images from different species. The same aSVG
images are also used by the web service of the Expression Atlas project. In addition, the
spatialHeatmap has its own repository called spatialHeatmap aSVG
Repository,
that stores custom aSVG files developed for this project (e.g. Figure
6).
If users cannot find a suitable aSVG image in these two repositories, we have developed a step-by-step SVG tutorial for creating custom aSVG images. For example, the BAR eFP browser at University of Toronto contains many anatomical images. These images can be used as templates in the SVG tutorial to construct custom aSVGs.
Moreover, in the future we will add more aSVGs to our repository, and users are welcome to deposit their own aSVGs to this repository to share them with the community.
7.2.2 Update aSVG features
To create and edit existing feature identifiers in aSVGs, the update_feature function
can be used. The demonstration below, creates an empty folder tmp.dir1 and copies
into it the homo_sapiens.brain.svg image provided by the spatialHeatmap
package. Subsequently, selected feature annotations are updated in this file.
tmp.dir1 <- paste0(normalizePath(tempdir(check=TRUE), winslash="/", mustWork=FALSE), '/shm1')
if (!dir.exists(tmp.dir1)) dir.create(tmp.dir1)
svg.hum <- system.file("extdata/shinyApp/example", 'homo_sapiens.brain.svg', package="spatialHeatmap")
file.copy(from=svg.hum, to=tmp.dir1, overwrite=TRUE) # Copy "homo_sapiens.brain.svg" file into 'tmp.dir1'Query the above aSVG with feature and species keywords, and return the resulting matches
in a data.frame.
feature.df <- return_feature(feature=c('frontal cortex', 'prefrontal cortex'), species=c('homo sapiens', 'brain'), dir=tmp.dir1, remote=FALSE, keywords.any=FALSE)
feature.dfSubsequently, create a character vector of new feature identifiers corresponding to
each of the returned features. In the following examples spaces in strings will be
filled with dots. This character vector must be added to the first column of the
feature data.frame. The latter is used by the update_feature function to look
up new features.
Sample code that creates new feature names and stores them in a character vector.
f.new <- c('prefrontal.cortex', 'frontal.cortex')Next, new features are added to the first column of the feature data.frame.
feature.df.new <- cbind(featureNew=f.new, feature.df)
feature.df.newFinally, the features are updated in the aSVG file(s) located in the tmp.dir1 directory.
update_feature(feature=feature.df.new, dir=tmp.dir1)8 Version Informaion
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] GEOquery_2.58.0 ExpressionAtlas_1.18.0
## [3] xml2_1.3.2 limma_3.46.0
## [5] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [7] GenomicRanges_1.42.0 GenomeInfoDb_1.26.0
## [9] IRanges_2.24.0 S4Vectors_0.28.0
## [11] BiocGenerics_0.36.0 MatrixGenerics_1.2.0
## [13] matrixStats_0.57.0 spatialHeatmap_1.0.0
## [15] knitr_1.30 BiocStyle_2.18.0
##
## loaded via a namespace (and not attached):
## [1] backports_1.1.10 rols_2.18.0 Hmisc_4.4-1
## [4] av_0.5.1 igraph_1.2.6 lazyeval_0.2.2
## [7] shinydashboard_0.7.1 splines_4.0.3 BiocParallel_1.24.0
## [10] ggplot2_3.3.2 digest_0.6.27 foreach_1.5.1
## [13] htmltools_0.5.0 magick_2.5.0 GO.db_3.12.0
## [16] fansi_0.4.1 magrittr_1.5 checkmate_2.0.0
## [19] memoise_1.1.0 cluster_2.1.0 doParallel_1.0.16
## [22] readr_1.4.0 fastcluster_1.1.25 annotate_1.68.0
## [25] prettyunits_1.1.1 jpeg_0.1-8.1 colorspace_1.4-1
## [28] blob_1.2.1 xfun_0.18 dplyr_1.0.2
## [31] crayon_1.3.4 RCurl_1.98-1.2 jsonlite_1.7.1
## [34] genefilter_1.72.0 impute_1.64.0 survival_3.2-7
## [37] iterators_1.0.13 glue_1.4.2 gtable_0.3.0
## [40] zlibbioc_1.36.0 XVector_0.30.0 DelayedArray_0.16.0
## [43] scales_1.1.1 DBI_1.1.0 edgeR_3.32.0
## [46] Rcpp_1.0.5 viridisLite_0.3.0 xtable_1.8-4
## [49] progress_1.2.2 htmlTable_2.1.0 gridGraphics_0.5-0
## [52] flashClust_1.01-2 foreign_0.8-80 bit_4.0.4
## [55] preprocessCore_1.52.0 Formula_1.2-4 rsvg_2.1
## [58] htmlwidgets_1.5.2 httr_1.4.2 gplots_3.1.0
## [61] RColorBrewer_1.1-2 ellipsis_0.3.1 farver_2.0.3
## [64] pkgconfig_2.0.3 XML_3.99-0.5 nnet_7.3-14
## [67] locfit_1.5-9.4 dynamicTreeCut_1.63-1 labeling_0.4.2
## [70] ggplotify_0.0.5 tidyselect_1.1.0 rlang_0.4.8
## [73] later_1.1.0.1 AnnotationDbi_1.52.0 visNetwork_2.0.9
## [76] munsell_0.5.0 tools_4.0.3 cli_2.1.0
## [79] generics_0.0.2 RSQLite_2.2.1 fastmap_1.0.1
## [82] evaluate_0.14 stringr_1.4.0 ggdendro_0.1.22
## [85] yaml_2.2.1 bit64_4.0.5 caTools_1.18.0
## [88] purrr_0.3.4 mime_0.9 compiler_4.0.3
## [91] rstudioapi_0.11 curl_4.3 plotly_4.9.2.1
## [94] png_0.1-7 tibble_3.0.4 geneplotter_1.68.0
## [97] stringi_1.5.3 ps_1.4.0 highr_0.8
## [100] lattice_0.20-41 Matrix_1.2-18 vctrs_0.3.4
## [103] pillar_1.4.6 lifecycle_0.2.0 BiocManager_1.30.10
## [106] data.table_1.13.2 bitops_1.0-6 grImport_0.9-3
## [109] httpuv_1.5.4 R6_2.4.1 latticeExtra_0.6-29
## [112] bookdown_0.21 promises_1.1.1 KernSmooth_2.23-17
## [115] gridExtra_2.3 codetools_0.2-16 assertthat_0.2.1
## [118] MASS_7.3-53 gtools_3.8.2 DESeq2_1.30.0
## [121] GenomeInfoDbData_1.2.4 hms_0.5.3 grid_4.0.3
## [124] rpart_4.1-15 tidyr_1.1.2 rmarkdown_2.5
## [127] rvcheck_0.1.8 shiny_1.5.0 WGCNA_1.69
## [130] base64enc_0.1-39 Funding
This project has been funded by NSF awards: PGRP-1546879, PGRP-1810468, PGRP-1936492.
10 References
Cardoso-Moreira, Margarida, Jean Halbert, Delphine Valloton, Britta Velten, Chunyan Chen, Yi Shao, Angélica Liechti, et al. 2019. “Gene Expression Across Mammalian Organ Development.” Nature 571 (7766):505–9.
Chang, Winston, and Barbara Borges Ribeiro. 2018. Shinydashboard: Create Dashboards with ’Shiny’. https://CRAN.R-project.org/package=shinydashboard.
Chang, Winston, Joe Cheng, JJ Allaire, Yihui Xie, and Jonathan McPherson. n.d. Shiny: Web Application Framework for R. https://CRAN.R-project.org/package=shiny.
Davis, Sean, and Paul Meltzer. 2007. “GEOquery: A Bridge Between the Gene Expression Omnibus (Geo) and Bioconductor.” Bioinformatics 14:1846–7.
Gatto, Laurent. 2019. Rols: An R Interface to the Ontology Lookup Service. http://lgatto.github.com/rols/.
Gautier, Laurent, Leslie Cope, Benjamin M. Bolstad, and Rafael A. Irizarry. 2004. “Affy—Analysis of Affymetrix Genechip Data at the Probe Level.” Bioinformatics 20 (3). Oxford, UK: Oxford University Press:307–15. https://doi.org/10.1093/bioinformatics/btg405.
Huber, W., V. J. Carey, R. Gentleman, S. An ders, M. Carlson, B. S. Carvalho, H. C. Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis Wit H Bioconductor.” Nature Methods 12 (2):115–21. http://www.nature.com/nmeth/journal/v12/n2/full/nmeth.3252.html.
Keays, Maria. 2019. ExpressionAtlas: Download Datasets from Embl-Ebi Expression Atlas.
Langfelder, P., and S. Horvath. 2012. “Fast R Functions for Robust Correlations and Hierarchical Clustering.” J. Stat. Softw. 46 (11). https://www.ncbi.nlm.nih.gov/pubmed/23050260.
Langfelder, Peter, and Steve Horvath. 2008. “WGCNA: An R Package for Weighted Correlation Network Analysis.” BMC Bioinformatics 9 (December):559.
Langfelder, Peter, Bin Zhang, and with contributions from Steve Horvath. 2016. DynamicTreeCut: Methods for Detection of Clusters in Hierarchical Clustering Dendrograms. https://CRAN.R-project.org/package=dynamicTreeCut.
Lekschas, Fritz, Harald Stachelscheid, Stefanie Seltmann, and Andreas Kurtz. 2015. “Semantic Body Browser: Graphical Exploration of an Organism and Spatially Resolved Expression Data Visualization.” Bioinformatics 31 (5):794–96.
Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for Rna-Seq Data with Deseq2.” Genome Biology 15 (12):550. https://doi.org/10.1186/s13059-014-0550-8.
Maag, Jesper L V. 2018. “Gganatogram: An R Package for Modular Visualisation of Anatograms and Tissues Based on Ggplot2.” F1000Res. 7 (September):1576.
Merkin, Jason, Caitlin Russell, Ping Chen, and Christopher B Burge. 2012. “Evolutionary Dynamics of Gene and Isoform Regulation in Mammalian Tissues.” Science 338 (6114):1593–9.
Morgan, Martin, Valerie Obenchain, Jim Hester, and Hervé Pagès. 2018. SummarizedExperiment: SummarizedExperiment Container.
Mungall, Christopher J, Carlo Torniai, Georgios V Gkoutos, Suzanna E Lewis, and Melissa A Haendel. 2012. “Uberon, an integrative multi-species anatomy ontology.” Genome Biol. 13 (1):R5. https://doi.org/10.1186/gb-2012-13-1-r5.
Muschelli, John, Elizabeth Sweeney, and Ciprian Crainiceanu. 2014. “BrainR: Interactive 3 and 4D Images of High Resolution Neuroimage Data.” R J. 6 (1):41–48.
Mustroph, Angelika, M Eugenia Zanetti, Charles J H Jang, Hans E Holtan, Peter P Repetti, David W Galbraith, Thomas Girke, and Julia Bailey-Serres. 2009. “Profiling Translatomes of Discrete Cell Populations Resolves Altered Cellular Priorities During Hypoxia in Arabidopsis.” Proc Natl Acad Sci U S A 106 (44):18843–8.
Papatheodorou, Irene, Nuno A Fonseca, Maria Keays, Y Amy Tang, Elisabet Barrera, Wojciech Bazant, Melissa Burke, et al. 2018. “Expression Atlas: gene and protein expression across multiple studies and organisms.” Nucleic Acids Res. 46 (D1):D246–D251. https://doi.org/10.1093/nar/gkx1158.
Prudencio, Mercedes, Veronique V Belzil, Ranjan Batra, Christian A Ross, Tania F Gendron, Luc J Pregent, Melissa E Murray, et al. 2015. “Distinct Brain Transcriptome Profiles in C9orf72-Associated and Sporadic ALS.” Nat. Neurosci. 18 (8):1175–82.
Ravasz, E, A L Somera, D A Mongru, Z N Oltvai, and A L Barabási. 2002. “Hierarchical Organization of Modularity in Metabolic Networks.” Science 297 (5586):1551–5.
Robinson, Mark D, Davis J McCarthy, and Gordon K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1):139–40. https://doi.org/10.1093/bioinformatics/btp616.
Waese, Jamie, Jim Fan, Asher Pasha, Hans Yu, Geoffrey Fucile, Ruian Shi, Matthew Cumming, et al. 2017. “EPlant: Visualizing and Exploring Multiple Levels of Data for Hypothesis Generation in Plant Biology.” Plant Cell 29 (8):1806–21.
Winter, Debbie, Ben Vinegar, Hardeep Nahal, Ron Ammar, Greg V Wilson, and Nicholas J Provart. 2007. “An ‘Electronic Fluorescent Pictograph’ Browser for Exploring and Analyzing Large-Scale Biological Data Sets.” PLoS One 2 (8):e718.