---
title: "Introduction to TreeSummarizedExperiment"
author:
- name: Ruizhu HUANG
affiliation:
- Institute of Molecular Life Sciences, University of Zurich.
- SIB Swiss Institute of Bioinformatics.
- name: Charlotte Soneson
affiliation:
- Institute of Molecular Life Sciences, University of Zurich.
- SIB Swiss Institute of Bioinformatics.
- name: Mark Robinson
affiliation:
- Institute of Molecular Life Sciences, University of Zurich.
- SIB Swiss Institute of Bioinformatics.
package: TreeSummarizedExperiment
output:
BiocStyle::html_document
vignette: >
%\VignetteIndexEntry{Introduction to TreeSE}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{knitr::rmarkdown}
editor_options:
chunk_output_type: console
bibliography: TreeSE_vignette.bib
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
```
# Introduction
The `TreeSummarizedExperiment` class is an extension of the
`SingleCellExperiment` class [@LunA2020]. It's used to store rectangular data
of experimental results as in a `SingleCellExperiment`, and also supports the
storage of a hierarchical structure and its link information to the rectangular
data.
# TreeSummarizedExperiment {#tse-class}
## Anatomy of TreeSummarizedExperiment
```{r strTSE, echo=FALSE, fig.cap= "The structure of the TreeSummarizedExperiment class."}
knitr::include_graphics("tse.png")
```
Compared to the `SingleCellExperiment` objects, `TreeSummarizedExperiment` has
four additional slots:
* `rowTree`: the hierarchical structure on the rows of the `assays`.
* `rowLinks`: the link information between rows of the `assays` and the `rowTree`.
* `colTree`: the hierarchical structure on the columns of the `assays`.
* `colLinks`: the link information between columns of the `assays` and the
`colTree`.
The `rowTree` and$/$or `colTree` can be left empty (`NULL`) if no trees are available; in this case, the `rowLinks` and$/$or `colLinks` are also set to `NULL`. All other `TreeSummarizedExperiment` slots are inherited from `SingleCellExperiment`.
The `rowTree` and `colTree` slots require the tree to be an object of the `phylo`
class. If a tree is available in an alternative format, it can often be converted to a `phylo` object using dedicated R packages (e.g., `r Biocpkg("treeio")` [@Wang2019]).
Functions in the `r Biocpkg("TreeSummarizedExperiment")` package fall in two main categories: operations on the `TreeSummarizedExperiment` object or operations on the tree (`phylo`) objects. The former includes constructors and accessors, and the latter serves as "pieces" to be assembled as accessors or functions that manipulate the `TreeSummarizedExperiment` object. Given that more than 200 R packages make use of the `phylo` class, there are many resources (e.g., `r CRANpkg("ape")`) for users to manipulate the small "pieces" in addition to those provided in `r Biocpkg("TreeSummarizedExperiment")`.
## Toy data
We generate a toy dataset that has observations of 6 entities collected from 4 samples as an example to show how to construct a `TreeSummarizedExperiment` object.
```{r}
library(TreeSummarizedExperiment)
# assays data (typically, representing observed data from an experiment)
assay_data <- rbind(rep(0, 4), matrix(1:20, nrow = 5))
colnames(assay_data) <- paste0("sample", 1:4)
rownames(assay_data) <- paste("entity", seq_len(6), sep = "")
assay_data
```
The information of entities and samples are given in the **row_data** and **col_data**, respectively.
```{r}
# row data (feature annotations)
row_data <- data.frame(Kingdom = "A",
Phylum = rep(c("B1", "B2"), c(2, 4)),
Class = rep(c("C1", "C2", "C3"), each = 2),
OTU = paste0("D", 1:6),
row.names = rownames(assay_data),
stringsAsFactors = FALSE)
row_data
# column data (sample annotations)
col_data <- data.frame(gg = c(1, 2, 3, 3),
group = rep(LETTERS[1:2], each = 2),
row.names = colnames(assay_data),
stringsAsFactors = FALSE)
col_data
```
The hierarchical structure of the `r nrow(assay_data)` entities and
`r ncol(assay_data)` samples are denoted as **row_tree** and **col_tree**,
respectively. The two trees are `phylo` objects randomly created with `rtree`
from the package `r CRANpkg("ape")`. Note that the row tree has 5 rather than 6
leaves; this is used later to show that multiple rows in the `assays` are allowed to map to a single node in the tree.
```{r}
library(ape)
# The first toy tree
set.seed(12)
row_tree <- rtree(5)
# The second toy tree
set.seed(12)
col_tree <- rtree(4)
# change node labels
col_tree$tip.label <- colnames(assay_data)
col_tree$node.label <- c("All", "GroupA", "GroupB")
```
We visualize the tree using the package `r Biocpkg("ggtree")` [@Yu2017]. Here, the
internal nodes of the **row_tree** have no labels as shown in Figure \@ref(fig:plot-rtree).
```{r plot-rtree, fig.cap="\\label{plot-rtree} The structure of the row tree. The node labels and the node numbers are in orange and blue text, respectively.", out.width="90%"}
library(ggtree)
library(ggplot2)
# Visualize the row tree
ggtree(row_tree, size = 2, branch.length = "none") +
geom_text2(aes(label = node), color = "darkblue",
hjust = -0.5, vjust = 0.7, size = 5.5) +
geom_text2(aes(label = label), color = "darkorange",
hjust = -0.1, vjust = -0.7, size = 5.5)
```
The **col_tree** has labels for internal nodes.
```{r plot-ctree, fig.cap="\\label{plot-ctree} The structure of the column tree. The node labels and the node numbers are in orange and blue text, respectively.", out.width="90%"}
# Visualize the column tree
ggtree(col_tree, size = 2, branch.length = "none") +
geom_text2(aes(label = node), color = "darkblue",
hjust = -0.5, vjust = 0.7, size = 5.5) +
geom_text2(aes(label = label), color = "darkorange",
hjust = -0.1, vjust = -0.7, size = 5.5)+
ylim(c(0.8, 4.5)) +
xlim(c(0, 2.2))
```
## The construction of `TreeSummarizedExperiment`
The `TreeSummarizedExperiment` class is used to store the toy data generated in the previous section:
**assay_data**, **row_data**, **col_data**, **col_tree** and **row_tree**. To
correctly store data, the link information between the rows (or columns) of
**assay_data** and the nodes of the **row_tree** (or **col_tree**) can be provided via a character vector `rowNodeLab` (or `colNodeLab`), with length equal to the number of rows (or columns) of the `assays`; otherwise the row (or column) names are used. Those columns
or rows with labels that are not present among the node labels of the tree are removed with
warnings. The link data between the `assays` tables and the tree data is
automatically generated in the construction.
The row and column trees can be included simultaneously in the construction.
Here, the column names of **assay_data** can be found in the node labels of the
column tree, which enables the link to be created between the column dimension
of **assay_data** and the column tree **col_tree**. If the row names of
**assay_data** are not in the node labels of **row_tree**, we would need to
provide their corresponding node labels (**row_lab**) to `rowNodeLab` in the
construction of the object. It is allowed to have multiple rows or/and columns
mapped to a node, for example, the same leaf label is used for the first two
rows in **row_lab**.
```{r}
# all column names could be found in the node labels of the column tree
all(colnames(assay_data) %in% c(col_tree$tip.label, col_tree$node.label))
# provide the node labels in rowNodeLab
tip_lab <- row_tree$tip.label
row_lab <- tip_lab[c(1, 1:5)]
row_lab
both_tse <- TreeSummarizedExperiment(assays = list(Count = assay_data),
rowData = row_data,
colData = col_data,
rowTree = row_tree,
rowNodeLab = row_lab,
colTree = col_tree)
```
```{r}
both_tse
```
When printing out **both_tse**, we see a similar message as
`SingleCellExperiment` with four additional lines for `rowLinks`, `rowTree`,
`colLinks` and `colTree`.
## The accessor functions
### Assays, rowData, colData, and metadata
For slots inherited from the `SingleCellExperiment` class, the accessors are
exactly the same as shown in `r Biocpkg("SingleCellExperiment")`.
```{r}
# to get the first table in the assays
(count <- assays(both_tse)[[1]])
```
```{r}
# to get row data
rowData(both_tse)
```
```{r}
# to get column data
colData(both_tse)
```
```{r}
# to get metadata: it's empty here
metadata(both_tse)
```
### rowLinks, colLinks {#linkData}
For new slots, we provide `rowTree` (and `colTree`) accessors to retrieve the row (column) trees, and `rowLinks` (and `colLinks`) to retrieve the link information between `assays` and nodes of the row (column) tree. If the tree is not available, the corresponding link data is `NULL`.
```{r}
# access trees
rowTree(both_tse)
colTree(both_tse)
```
```{r}
# access the link data
(r_link <- rowLinks(both_tse))
(c_link <- colLinks(both_tse))
```
The link data objects are of the `LinkDataFrame` class, which extends the `DataFrame` class with the restriction that it has at least four columns:
* `nodeLab`: the labels of nodes on the tree
* `nodeLab_alias`: the alias labels of nodes on the tree
* `nodeNum`: the numbers of nodes on the tree
* `isLeaf`: whether the node is a leaf node
More details about
the `DataFrame` class could be found in the `r Biocpkg("S4Vectors")` R/Bioconductor package.
```{r}
class(r_link)
showClass("LinkDataFrame")
```
The link data is automatically generated when constructing the
`TreeSummarizedExperiment` object. We highly recommend users not to modify it
manually; otherwise the link might be broken. For R packages developers, we show
in the Section \@ref(modifyLink) about how to update the link.
## The subseting function
A `TreeSummarizedExperiment` object can be subset in two different ways: `[` to
subset by rows or columns, and `subsetByNode` to subset by nodes of a tree. To
keep track of the original data, the `rowTree` and `colTree` stay the same after subsetting, while `rowLinks` and `rowData` are updated accordingly.
```{r}
sub_tse <- both_tse[1:2, 1]
sub_tse
```
```{r}
# the row data
rowData(sub_tse)
# the row link data
rowLinks(sub_tse)
```
```{r}
# The first four columns are from colLinks data and the others from colData
cbind(colLinks(sub_tse), colData(sub_tse))
```
To subset by nodes, we specify the node by its node label or node number. Here, *entity1* and *entity2* are both mapped to the same node `t3`, so both of them are retained.
```{r}
node_tse <- subsetByNode(x = both_tse, rowNode = "t3")
rowLinks(node_tse)
```
Subsetting simultaneously in both dimensions is also allowed.
```{r}
node_tse <- subsetByNode(x = both_tse, rowNode = "t3",
colNode = c("sample1", "sample2"))
assays(node_tse)[[1]]
```
## Changing the tree
The current tree can be replaced by a new one using `changeTree`. If the hierarchical information is available as a `data.frame` with each column representing a taxonomic level (e.g., *row_data*), we provide `toTree` to convert it into a `phylo` object.
```{r plot-taxa2phylo, fig.cap="\\label{plot-taxa2phylo} The structure of the taxonomic tree that is generated from the taxonomic table.", out.width="90%"}
# The toy taxonomic table
(taxa <- rowData(both_tse))
# convert it to a phylo tree
taxa_tree <- toTree(data = taxa)
# Viz the new tree
ggtree(taxa_tree)+
geom_text2(aes(label = node), color = "darkblue",
hjust = -0.5, vjust = 0.7, size = 5.5) +
geom_text2(aes(label = label), color = "darkorange",
hjust = -0.1, vjust = -0.7, size = 5.5) +
geom_point2()
```
A mapping to match nodes of the two trees is required if nodes are labeled differently.
```{r}
taxa_tse <- changeTree(x = both_tse, rowTree = taxa_tree,
rowNodeLab = taxa[["OTU"]])
taxa_tse
rowLinks(taxa_tse)
```
## Aggregation
Since it may be of interest to report or analyze observed data on multiple resolutions based on the provided tree, the `TreeSummarizedExperiment` package offers functionionality to flexibly aggregate data to different levels of a tree.
### The column dimension {#aggCol}
Here, we show the aggregation along the column dimension. The desired
aggregation level is given in the `colLevel` argument, which can be specified via the
node label (orange text in Figure \@ref(fig:plot-ctree)) or the node number (blue text in Figure \@ref(fig:plot-ctree)). We could further specify how to aggregate
via the argument `FUN`.
```{r}
# use node labels to specify colLevel
agg_col <- aggValue(x = taxa_tse,
colLevel = c("GroupA", "GroupB"),
FUN = sum)
# or use node numbers to specify colLevel
agg_col <- aggValue(x = taxa_tse, colLevel = c(6, 7), FUN = sum)
```
```{r}
assays(agg_col)[[1]]
```
The `rowData` does not change, but the `colData` adjusts with the change of the
`assays` table. For example, the column **group** has the `A` value for
`GroupA` because the descendant nodes of `GroupA` all have the value `A`; the
column **gg** has the `NA` value for `GroupA` because the descendant nodes of
`GroupA` have different values, (1 and 2).
```{r}
# before aggregation
colData(taxa_tse)
# after aggregation
colData(agg_col)
```
The `colLinks` is updated to link the new rows of `assays` tables and the column
tree.
```{r}
# the link data is updated
colLinks(agg_col)
```
From Figure \@ref(fig:plot-ctree), nodes 6 and 7 are
labeled with `GroupA` and `GroupB`, respectively, which agrees with the
column link data.
### The row dimension {#aggRow}
Similarly, we could aggregate the data to the phylum level by providing the names of the internal nodes that represent the phylum level (see `taxa_one` below).
```{r}
# the phylum level
taxa <- c(taxa_tree$tip.label, taxa_tree$node.label)
(taxa_one <- taxa[startsWith(taxa, "Phylum:")])
# aggregation
agg_taxa <- aggValue(x = taxa_tse, rowLevel = taxa_one, FUN = sum)
assays(agg_taxa)[[1]]
```
The user is nonetheless free to choose nodes from different taxonomic ranks. Note that not all rows in the original table are included in one of the aggregated rows. Similarly, it is possible for a row to contribute to multiple aggregated rows
```{r}
# A mixed level
taxa_mix <- c("Class:C3", "Phylum:B1")
agg_any <- aggValue(x = taxa_tse, rowLevel = taxa_mix, FUN = sum)
rowData(agg_any)
```
### Both dimensions
The aggregation on both dimensions could be performed in one step using the same
function specified via `FUN`. If different functions are required for different
dimensions, the aggregation should be performed in two steps because the aggregation order matters.
```{r}
agg_both <- aggValue(x = both_tse, colLevel = c(6, 7),
rowLevel = 7:9, FUN = sum)
```
As expected, we obtain a table with 3 rows (`rowLevel = 7:9`) and 2 columns
(`colLevel = c(6, 7)`).
```{r}
assays(agg_both)[[1]]
```
## Functions operating on the `phylo` object.
Next, we highlight some functions to manipulate and/or to extract information
from the `phylo` object. Further operations can be found in other packages, such
as `r CRANpkg("ape")` [@ape2019], `r CRANpkg("tidytree")`[@R-tidytree]. These functions are
useful when users want to customize functions for the `TreeSummarizedExperiment`
class.
To show these functions, we use the tree shown in Figure \@ref(fig:plot-exTree).
```{r plot-exTree, fig.cap= "\\label{plot-exTree} An example tree with node labels and numbers in black and orange texts, respectively.", out.width="90%"}
data("tinyTree")
ggtree(tinyTree, branch.length = "none") +
geom_text2(aes(label = label), hjust = -0.1, size = 5.5) +
geom_text2(aes(label = node), vjust = -0.8,
hjust = -0.2, color = 'orange', size = 5.5)
```
### Conversion of the node label and the node number
The translation between the node labels and node numbers can be achieved by the function `convertNode`.
```{r}
convertNode(tree = tinyTree, node = c(12, 1, 4))
```
```{r}
convertNode(tree = tinyTree, node = c("t4", "Node_18"))
```
### Find the descendants
To get descendants that are at the leaf level, we could set the argument
`only.leaf = TRUE` for the function `findDescendant`.
```{r}
# only the leaf nodes
findDescendant(tree = tinyTree, node = 17, only.leaf = TRUE)
```
When `only.leaf = FALSE`, all descendants are returned.
```{r}
# all descendant nodes
findDescendant(tree = tinyTree, node = 17, only.leaf = FALSE)
```
### More functions
We list some functions that might also be useful in Table \@ref(tab:phyloFun).
More are available in the package, and we encourage users to contribute their
functions that might be helpful for others.
| Functions | Goal |
| ----------- | ------------------------------------------------------------ |
| printNode | print out the information of nodes |
| countNode | count the number of nodes |
| distNode | give the distance between a pair of nodes |
| matTree | list paths of a tree |
| findAncestor| find ancestor nodes |
| findChild | find child nodes |
| findSibling | find sibling nodes |
| shareNode | find the first node shared in the paths of nodes to the root |
| unionLeaf | find the union of descendant leaves |
| trackNode | track nodes by adding alias labels to a phylo object |
| isLeaf | test whether a node is a leaf node |
: (\#tab:phyloFun) A table lists some functions operating on the `phylo` object that are available in the `TreeSummarizedExperiment`.
## Custom functions for the `TreeSummarizedExperiment` class {#modifyLink}
Here, we show examples on how to write custom functions for `TreeSummarizedExperiment`
objects. To extract data of specific leaves, we created a function
`subsetByLeaf` by combining functions working on the `phylo` class (e.g.,
`convertNode`, `keep.tip`, `trackNode`, `isLeaf`) with the accessor function `subsetByNode`.
Here, `convertNode`, `trackNode` and `isLeaf` are available in `TreeSummarizedExperiment`,
and `keep.tip` is from the `r CRANpkg("ape")` package. Since the node number of a node is changed after pruning a tree with `keep.tip`, `trackNode` is provided to track the node and further update the link between the data and the new tree.
```{r}
# tse: a TreeSummarizedExperiment object
# rowLeaf: specific leaves
subsetByLeaf <- function(tse, rowLeaf) {
# if rowLeaf is provided as node labels, convert them to node numbers
if (is.character(rowLeaf)) {
rowLeaf <- convertNode(tree = rowTree(tse), node = rowLeaf)
}
# subset data by leaves
sse <- subsetByNode(tse, rowNode = rowLeaf)
# update the row tree
## -------------- new tree: drop leaves ----------
oldTree <- rowTree(sse)
newTree <- ape::keep.tip(phy = oldTree, tip = rowLeaf)
## -------------- update the row link ----------
# track the tree
track <- trackNode(oldTree)
track <- ape::keep.tip(phy = track, tip = rowLeaf)
# row links
rowL <- rowLinks(sse)
rowL <- DataFrame(rowL)
# update the row links:
# 1. use the alias label to track and updates the nodeNum
# 2. the nodeLab should be updated based on the new tree using the new
# nodeNum
# 3. lastly, update the nodeLab_alias
rowL$nodeNum <- convertNode(tree = track, node = rowL$nodeLab_alias,
message = FALSE)
rowL$nodeLab <- convertNode(tree = newTree, node = rowL$nodeNum,
use.alias = FALSE, message = FALSE)
rowL$nodeLab_alias <- convertNode(tree = newTree, node = rowL$nodeNum,
use.alias = TRUE, message = FALSE)
rowL$isLeaf <- isLeaf(tree = newTree, node = rowL$nodeNum)
rowNL <- new("LinkDataFrame", rowL)
## update the row tree and links
BiocGenerics:::replaceSlots(sse,
rowLinks = rowNL,
rowTree = list(phylo = newTree))
}
```
The row tree is updated after the subsetting. It now has only two leaves, `t2` and `t3`.
```{r}
(both_sse <- subsetByLeaf(tse = both_tse, rowLeaf = c("t2", "t3")))
rowLinks(both_sse)
```
# Session Info
```{r}
sessionInfo()
```
# Reference