library(BiocStyle)
The r Rpackage("SEtools")
package is a set of convenience functions for the Bioconductor class r Biocpkg("SummarizedExperiment")
. It facilitates merging, melting, and plotting SummarizedExperiment
objects.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("SEtools")
Or, to install the latest development version:
BiocManager::install("plger/SEtools")
To showcase the main functions, we will use an example object which contains (a subset of) whole-hippocampus RNAseq of mice after different stressors:
suppressPackageStartupMessages({ library(SummarizedExperiment) library(SEtools) }) data("SE", package="SEtools") SE
This is taken from Floriou-Servou et al., Biol Psychiatry 2018.
There are two main wrappers for plotting heatmaps from SummarizedExperiment
objects:
sehm
function uses the pheatmap enginesechm
function uses the ComplexHeatmap engineBoth functions were made to function very similarly, but the sechm
function is especially useful to combine heatmaps (for instance, from different SummarizedExperiment
objects). We'll showcase sehm
(the main functionalities being replicable with sechm
), and will then provide examples of multiple heatmaps.
The sehm
function simplifies the generation of heatmaps from SummarizedExperiment
. It uses r CRANpkg("pheatmap")
, so any argument supported by it can in principle be passed:
g <- c("Egr1", "Nr4a1", "Fos", "Egr2", "Sgk1", "Arc", "Dusp1", "Fosb", "Sik1") sehm(SE, genes=g) sehm(SE, assayName="logcpm", genes=g, do.scale=TRUE)
When scaling data, the function will automatically center the colour scale around zero, and handle the extreme values (0.5\% percentile on each side) in a non-linear fashion to retain a useful visualization.
This behavior can be manually controlled via the breaks
parameter (either setting it to FALSE, to a percentile until which the scale should be linear, of manually inputting breaks).
Annotation from the object's rowData
and colData
can be plotted simply by specifying the column name (some will be shown by default if found):
sehm(SE, assayName="logcpm", genes=g, do.scale=TRUE, anno_rows="meanTPM")
These can also be used to create gaps:
sehm(SE, genes=g, do.scale=TRUE, anno_rows="meanTPM", gaps_at="Condition")
The specific assay to use for plotting can be specified with the assayName
argument.
By default, rows are sorted not with hierarchical clustering, but from the angle on a MDS plot, which tends to give nicer results than bottom-up hierarchical clustering. This can be disabled using sortRowsOn=NULL
or cluster_rows=TRUE
(to avoid any row reordering and use the order given, use sortRowsOn=NULL, cluster_rows=FALSE
). Column clustering is disabled by default, but this can be changed with cluster_cols=TRUE
.
It is common to cluster features into groups, and such a clustering can be used simultaneously with row sorting using the toporder
argument. For instance:
lfcs <- assays(SE)$logcpm-rowMeans(assays(SE)$logcpm[,which(SE$Condition=="Homecage")]) rowData(SE)$cluster <- as.character(kmeans(lfcs,4)$cluster) sehm(SE, genes=g, do.scale=TRUE, anno_rows="cluster", toporder="cluster", gaps_at="Condition")
For some arguments (for instance colors), if they are not specified in the function call, SEtools
will try to see whether the object itself contains it, or whether the corresponding global options have been set, before using default colors. This means that if, in the context of a given project, the same colors are repeatedly being used, they can be specified a single time, and all subsequent plots will be affected.
Storing colors in the object:
metadata(SE)$hmcols <- c("purple","white","gold") ancols <- list( Condition=c( Homecage="#DB918B", Handling="#B86FD3", Restraint="#A9CED5", Swim="#B5DF7C" ) ) metadata(SE)$anno_colors <- ancols sehm(SE, g, do.scale = TRUE)
Using the global options:
options("SEtools_def_hmcols"=c("white","grey","black")) options("SEtools_def_anno_colors"=ancols) sehm(SE, g, do.scale = TRUE)
At the moment, the following arguments can be set as global options:
assayName
, hmcols
, anno_columns
, anno_rows
, anno_colors
, gaps_at
, breaks
.
Options must be set with the prefix SEtools_def_
, followed by the name of the argument.
To remove the predefined colors:
resetAllSEtoolsOptions() metadata(SE)$hmcols <- NULL metadata(SE)$anno_colors <- NULL
In order of priority, the arguments in the function call trump the object's metadata, which trumps the global options.
The sechm
function works like the sehm
function, but the fact that it outputs a Heatmap
object from ComplexHeatmap means that these can be easily combined:
sechm(SE, g, do.scale = TRUE) + sechm(SE, g, do.scale = FALSE)
However, doing so involves manual work to ensure that the labels and colors are nice and coherent, and that the rows names match. As a convenience, we provide the crossHm
function to handle these issues. crossHm
works with a list of SummarizedExperiment
objects:
# we build another SE object: SE2 <- SE assays(SE2)$logcpm <- jitter(assays(SE2)$logcpm, factor=1000) crossHm(list(SE1=SE, SE2=SE2), g, do.scale = TRUE)
A unique color scale can be enforced:
crossHm(list(SE1=SE, SE2=SE2), g, do.scale = TRUE, uniqueScale = TRUE)
se1 <- SE[,1:10] se2 <- SE[,11:20] se3 <- mergeSEs( list(se1=se1, se2=se2) ) se3
All assays were merged, along with rowData and colData slots.
By default, row z-scores are calculated for each object when merging. This can be prevented with:
se3 <- mergeSEs( list(se1=se1, se2=se2), do.scale=FALSE)
If more than one assay is present, one can specify a different scaling behavior for each assay:
se3 <- mergeSEs( list(se1=se1, se2=se2), use.assays=c("counts", "logcpm"), do.scale=c(FALSE, TRUE))
It is also possible to merge by rowData columns, which are specified through the mergeBy
argument.
In this case, one can have one-to-many and many-to-many mappings, in which case two behaviors are possible:
aggFun
, the features of each object will by aggregated by mergeBy
using this function before merging.rowData(se1)$metafeature <- sample(LETTERS,nrow(se1),replace = TRUE) rowData(se2)$metafeature <- sample(LETTERS,nrow(se2),replace = TRUE) se3 <- mergeSEs( list(se1=se1, se2=se2), do.scale=FALSE, mergeBy="metafeature", aggFun=median) sehm(se3, genes=row.names(se3))
A single SE can also be aggregated by using the aggSE
function:
se1b <- aggSE(se1, by = "metafeature") se1b
If the aggregation function(s) are not specified, aggSE
will try to guess decent aggregation functions from the assay names.
To facilitate plotting features with r CRANpkg("ggplot2")
, the meltSE
function combines assay values along with row/column data:
d <- meltSE(SE, genes=g[1:4]) head(d) suppressPackageStartupMessages(library(ggplot2)) ggplot(d, aes(Condition, counts, fill=Condition)) + geom_violin() + facet_wrap(~feature, scale="free")
Calculate an assay of log-foldchanges to the controls:
SE <- log2FC(SE, fromAssay="logcpm", controls=SE$Condition=="Homecage")
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.