---
title: "ExpressionSet and SummarizedExperiment"
author: "Mikhail Dozmorov"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 3
    number_sections: true
    theme: united
    highlight: tango
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, 
                      warning = FALSE, 
                      message = FALSE,
                      eval = TRUE)  # Set to TRUE when running
```

# Bioconductor basics

Bioconductor provides software for the analysis and comprehension of high-throughput genomic data. It extends R with specialized data containers such as the `ExpressionSet`.

To install Bioconductor:

```{r eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install()
```

To install additional Bioconductor packages, use `BiocManager::install("package_name")`.

# Bioconductor basics

Once the base Bioconductor packages have been installed, you can access the vignettes for a specific package as follows:

```{r eval=FALSE}
# BiocManager::install("Biobase")
library("Biobase")
openVignette()
```

```
Please select a vignette:  

 1: airway - Airway smooth muscle cells
 2: Biobase - An introduction to Biobase and ExpressionSets
 3: Biobase - esApply Introduction
 4: Biobase - Notes for eSet developers
```

Press "1" to read the first one - it is the foundation of genomics data formats used in R. Or, press "0" to quit.

# Most common Bioconductor objects

| Object                               | Origin / Era                        | Main Use                                             | Structure of Assay Data                                     | Sample Annotations                    | Feature Annotations                  | Genomic Coordinates             | Metadata                      | Notes                                                                               |
| ------------------------------------ | ----------------------------------- | ---------------------------------------------------- | ----------------------------------------------------------- | ------------------------------------- | ------------------------------------ | ------------------------------- | ----------------------------- | ----------------------------------------------------------------------------------- |
| **ExpressionSet**                    | Early Bioconductor (microarray era) | Microarray data, early RNA-seq                       | Matrix (features × samples)                                 | `phenoData` (`AnnotatedDataFrame`)    | `featureData` (`AnnotatedDataFrame`) | ❌ Not built-in                  | `experimentData` (MIAME)      | First standardized container; still widely used because many packages expect it.    |
| **SummarizedExperiment (SE)**        | Modern Bioconductor                 | RNA-seq, proteomics, general omics                   | One or more assays (matrix-like, e.g., counts, TPM, logCPM) | `colData` (`DataFrame`)               | `rowData` (`DataFrame`)              | ❌ Not required                  | `metadata` list               | Designed to be flexible beyond microarrays, supports multiple assays.               |
| **RangedSummarizedExperiment (RSE)** | Subclass of SE                      | Genomics data (RNA-seq, ChIP-seq, ATAC-seq)          | Same as SE                                                  | Same as SE                            | `rowRanges` (`GRanges`)              | ✅ Yes, built-in via `rowRanges` | Same as SE                    | Extends SE with proper genomic coordinates per feature (chrom, start, end, strand). |
| **SingleCellExperiment (SCE)**       | Subclass of SE                      | Single-cell RNA-seq, multimodal                      | Multiple assays (counts, logcounts, etc.)                   | `colData` for per-cell annotations    | `rowData` for gene annotations       | Optional (`rowRanges`)          | `reducedDims`, `altExps`      | Adds dimensionality reduction slots, alternative experiments, useful for scRNA-seq. |
| **MultiAssayExperiment (MAE)**       | Higher-level container              | Multi-omics (RNA-seq + ATAC-seq + methylation, etc.) | List of multiple SE/RSE objects                             | Unified sample map across experiments | Experiment-specific                  | Depends on components           | `experimentData` & `metadata` | Organizes several SE-like objects under one umbrella.                               |

# Key Points

* **ExpressionSet**: Legacy, microarray-focused, standardized assays + sample + feature + experiment metadata.
* **SummarizedExperiment**: General-purpose, flexible, supports multiple assays, modern replacement.
* **RangedSummarizedExperiment**: Adds genomic ranges to rows → best for sequencing data with genomic positions.
* **SingleCellExperiment**: Tailored for single-cell; adds PCA/UMAP/t-SNE slots and multimodal support.
* **MultiAssayExperiment**: Higher-level container that combines multiple SEs/RSEs across omics types.

# ExpressionSet

Recall that objects in R can be either a `vector`, `factor`, `matrix`, `array`, `data.frame`, `list`, or `ts`.  
The **Biobase** package of the Bioconductor project is fundamental, and established new objects that can be used to store genomic data.  
An `ExpressionSet` is a container that organizes:

- `assayData` — expression matrix (features × samples)  
- `phenoData` — sample-level annotations (phenotypes, treatments, conditions)  
- `featureData` — feature-level annotations (e.g., genes, probes)  
- `protocolData` — information about the processing protocol  
- `experimentData` — experiment-level metadata  

Let’s look at an `ExpressionSet` object, built from the **airway** dataset.

```{r}
library(Biobase)

data("sample.ExpressionSet")
sample.ExpressionSet
```

For adventurous, let's peek under the hood to see the slots of the ExpressionSet object. Access them as `sample.ExpressionSet@experimentData`

## assayData (gene expression)

First, the most important part of the high-throughput genomic experiment is the matrix of expression values. The underlying structure of an expression matrix in Bioconductor is that the probes (i.e., genes) are in rows while the samples are in columns. Let’s read in an example expression matrix and then store it as an ExpressionSet. Once created, `exprs` is the extractor function that is used to access the expression values.

Let's read in a gene expression matrix.

```{r}
expression <- exprs(sample.ExpressionSet) # This is how to get to the expression matrix !

dim(expression)
class(expression)
rownames(expression)[1:4]
colnames(expression)[1:4]
expression[1:5, 1:3]
```

Having just expression values, we can construct minimal expression set.

```{r}
minimalSet <- ExpressionSet(assayData = as.matrix(expression))
minimalSet
exprs(minimalSet)[1:5, 1:3] 
featureNames(minimalSet)[1:4]
```

## phenoData (sample annotations)

Phenotypic data provides information about the samples, such as normal/abnormal, age, gender, etc. The phenotypic data is represented such that samples appear in rows while the variables appear in columns. Notice that when including phenotypic data in an ExpressionSet, the `row.names` in the `phenoData` must match the sample names in the expression matrix.

```{r}
characteristics <- pData(sample.ExpressionSet) # read.csv("lab/data/phenodata.csv", row.names = 1)
head(characteristics)
summary(characteristics)
all.equal(rownames(characteristics), colnames(expression))
```

You will get a warning if there is a mismatch. Before including the `phenoData` into the `ExpressionSet`, we may add some documentation describing information about each covariate (what does the variable name represent, what units the covariates are measure in, etc). This is done by creating a metadata table.

```{r}
# metadata must describe all three columns
metadata <- data.frame(
  labelDescription = c("Patient sex (Male or Female)",
                       "Sample type (Control or Case)",
                       "Phenotype score (numeric)"),
  row.names = c("sex", "type", "score")
)
metadata

# build AnnotatedDataFrame
phenoChar <- new("AnnotatedDataFrame", 
                 data = characteristics, 
                 varMetadata = metadata)

phenoChar
pData(phenoChar)[1:5, ]
pData(phenoChar)$sex[1:5]
```

Once a phenoData set is created, it can be accessed using the pData accessor function. Adding phenoData to samples from your ExpressionSet but ensure the phenotypic characteristics stored with it are properly aligned.

```{r}
anotherSet <- ExpressionSet(assayData = as.matrix(expression), phenoData = phenoChar)
anotherSet
males <- anotherSet[, pData(anotherSet)$Gender == "Male"]
pData(males)$Gender
```

The following code shows what happens when the phenotypic and expression data do not include matching sample names (output suppressed).

```{r eval=FALSE}
phony.pheno <- characteristics
rownames(phony.pheno)[1] <- "wrong.sample.name"
phenoPhony <- new("AnnotatedDataFrame", data = phony.pheno, varMetadata = metadata)
phony.pheno[1:3, ]
pData(phenoPhony)[1:3, ]
errorSet <- ExpressionSet(assayData = as.matrix(expression), phenoData = phenoPhony)
```

## Annotation (featureData, annotation)

After an analysis, one is usually left with cryptic manufacturer labels of the probes that were significant in your data analysis. To provide meaning to these probes, annotations represent meta data about the probes. The annotation package provides some basic tools for annotation packages.

```{r}
library(hgu95av2.db)
library(annotate)

# sample.ExpressionSet already includes exprs + phenoData
# just update its annotation slot
withannoSet <- sample.ExpressionSet
annotation(withannoSet) <- "hgu95av2.db"
withannoSet

# Get all probe IDs
all_probes <- featureNames(withannoSet)

# Remove spike-ins (those starting with "AFFX-")
gene_probes <- all_probes[!grepl("^AFFX-", all_probes)]

library(annotate)
# Map symbols/Entrez IDs for actual genes
gene_symbols <- getSYMBOL(gene_probes, "hgu95av2.db")
gene_entrez  <- getEG(gene_probes, "hgu95av2.db")

# Inspect first few
head(gene_symbols, n = 20)
head(gene_entrez)
```

Create `featureData` from scratch.

```{r}
# Assume `expression` is your expression matrix
probes <- rownames(expression)

# Create a data frame with some basic annotations
fdata_df <- data.frame(
  ProbeID = probes,
  Symbol  = getSYMBOL(probes, "hgu95av2.db"),
  Entrez  = getEG(probes, "hgu95av2.db"),
  stringsAsFactors = FALSE
)

# Convert to AnnotatedDataFrame
featureData_adf <- AnnotatedDataFrame(fdata_df)
featureData_adf
```



## experimentData

Data about the experiment can be stored in the `experimentData` slot.

```{r}
library(Biobase)

experimentData <- new("MIAME", 
                      name      = "Bioconductor Example Author", 
                      lab       = "Bioconductor Labs", 
                      contact   = "bioc-support@bioconductor.org", 
                      title     = "Example Microarray Study using sample.ExpressionSet", 
                      abstract  = "Toy dataset included in Biobase demonstrating structure of ExpressionSet with microarray expression data.", 
                      url       = "https://bioconductor.org/packages/Biobase", 
                      pubMedIds = "NA", 
                      other     = list(notes = "This dataset contains 500 features (probes) across 26 samples, including spike-in controls."))

experimentData
abstract(experimentData)
notes(experimentData)
```

### Putting it all together

```{r}
# Construct the ExpressionSet with featureData
withexpSet <- ExpressionSet(
  assayData     = as.matrix(expression),
  phenoData     = phenoChar,
  featureData   = featureData_adf,
  annotation    = "hgu95av2.db",
  experimentData= experimentData
)

withexpSet
experimentData(withexpSet)
abstract(experimentData(withexpSet))

exprs(withexpSet)[1:5, 1:5]
pData(withexpSet)
featureData(withexpSet)
```

# SummarizedExperiment

The next generation of an object that can hold annotated 'omics' data is `SummarizedExperiment`. It is not limited to genes, but instead holds information about genomic regions of interest.

```{r}
library(SummarizedExperiment)
?SummarizedExperiment
```

We'll look at an example of the `SummarizedExperiment` object in the **airway** package. The loaded data is a `SummarizedExperiment`, which contains RNA-seq read counts for human airway smooth muscle cells treated with dexamethasone or control. The `SummarizedExperiment` object has **63,677 rows**, which are genes (features), and **8 columns**, which are samples. The count data are stored in the assay called `counts`. Each row has metadata in `rowData(se)`, including Ensembl gene IDs, and each column has sample-level metadata in `colData(se)`, including cell line, treatment, and other phenotypic information. This object structure allows easy access to expression values, feature annotations, and sample annotations in a consistent way.

```{r}
### SummarizedExperiment demo with airway
library(airway)

# Load airway dataset (RNA-seq counts)
data("airway")
airway
se <- airway  # Treat as SummarizedExperiment
```

`assay` function can be used get access to the counts of RNA sequencing reads. `colData` function , the column data, is equivalent to the `pData` on the `ExpressionSet`. Each row in this data frame corresponds to a column in the `SummarizedExperiment`. We can see that there are indeed 27 rows here, which give information about the columns. Each sample in this case is treated with two treatments or control and we can see the number of replicates for each, using the as.numeric function again.

```{r}
# Dimension of the SummarizedExperiment
dim(se)
# Get access to the counts of RNA sequencing reads, using assay function.
assay(se)[1:3,1:3]
# Dimensions of this assay is a matrix, which has the same dimensions as the SummarizedExperiment.
dim(assay(se))
# Get information about samples
colData(se)[1:3,1:6]
# dimension of column data
dim(colData(se))
# characteristics of the samples
names(colData(se))
# Get access to treatment column of sample characteristics
colData(se)$dex
```

See [https://bioconductor.org/packages/devel/bioc/vignettes/SummarizedExperiment/inst/doc/SummarizedExperiment.html](https://bioconductor.org/packages/devel/bioc/vignettes/SummarizedExperiment/inst/doc/SummarizedExperiment.html) for the full description.

# Diagnostics

```{r session_info}
library(dplyr)
library(pander)
diagnostics <- devtools::session_info()
platform <- data.frame(diagnostics$platform %>% unlist, stringsAsFactors = FALSE)
colnames(platform) <- c("description")
pander(platform)

packages <- as.data.frame(diagnostics$packages)
pander(packages[ packages$`*` == "*", ])
```
