In waldronlab/MultiAssayExperiment: Software for the integration of multi-omics experiments in Bioconductor

library(MultiAssayExperiment)
library(S4Vectors)

This quick-start guide shows key features of MultiAssayExperiment using a subset of the TCGA adrenocortical carcinoma (ACC) dataset. This dataset provides five assays on 92 patients, although all five assays were not performed for every patient:

1. RNASeq2GeneNorm: gene mRNA abundance by RNA-seq
2. gistict: GISTIC genomic copy number by gene
3. RPPAArray: protein abundance by Reverse Phase Protein Array
4. Mutations: non-silent somatic mutations by gene
5. miRNASeqGene: microRNA abundance by microRNA-seq.

data(miniACC)
miniACC

Component slots

colData - information on biological units

A DataFrame describing the characteristics of biological units, for example clinical data for patients. In the prepared datasets from The Cancer Genome Atlas, each row is one patient and each column is a clinical, pathological, subtype, or other variable. The $ function provides a shortcut for accessing or setting colData columns.

colData(miniACC)[1:4, 1:4]
table(miniACC$race)

Key points: One row per patient Each row maps to zero or more observations in each experiment in the ExperimentList, below.

ExperimentList - experiment data

A base list or ExperimentList object containing the experimental datasets for the set of samples collected. This gets converted into a class ExperimentList during construction.

experiments(miniACC)

Key points: One matrix-like dataset per list element (although they do not even need to be matrix-like, see for example the RaggedExperiment package) One matrix column per assayed specimen. Each matrix column must correspond to exactly one row of colData: in other words, you must know which patient or cell line the observation came from. However, multiple columns can come from the same patient, or there can be no data for that patient. Matrix rows correspond to variables, e.g. genes or genomic ranges ExperimentList elements can be genomic range-based (e.g. SummarizedExperiment::RangedSummarizedExperiment-class or RaggedExperiment::RaggedExperiment-class) or ID-based data (e.g. SummarizedExperiment::SummarizedExperiment-class, Biobase::eSet-class base::matrix-class, DelayedArray::DelayedArray-class, and derived classes) * Any data class can be included in the ExperimentList, as long as it supports: single-bracket subsetting ([), dimnames, and dim. Most data classes defined in Bioconductor meet these requirements.

sampleMap - relationship graph

sampleMap is a graph representation of the relationship between biological units and experimental results. In simple cases where the column names of ExperimentList data matrices match the row names of colData, the user won't need to specify or think about a sample map, it can be created automatically by the MultiAssayExperiment constructor. sampleMap is a simple three-column DataFrame:

1. assay column: the name of the assay, and found in the names of ExperimentList list names
2. primary column: identifiers of patients or biological units, and found in the row names of colData
3. colname column: identifiers of assay results, and found in the column names of ExperimentList elements Helper functions are available for creating a map from a list. See ?listToMap

sampleMap(miniACC)

Key points: relates experimental observations (colnames) to colData permits experiment-specific sample naming, missing, and replicate observations

back to top

metadata

Metadata can be used to keep additional information about patients, assays performed on individuals or on the entire cohort, or features such as genes, proteins, and genomic ranges. There are many options available for storing metadata. First, MultiAssayExperiment has its own metadata for describing the entire experiment:

metadata(miniACC)

Additionally, the DataFrame class used by sampleMap and colData, as well as the ExperimentList class, similarly support metadata. Finally, many experimental data objects that can be used in the ExperimentList support metadata. These provide flexible options to users and to developers of derived classes.

Subsetting

Single bracket [

In pseudo code below, the subsetting operations work on the rows of the following indices: 1. i experimental data rows 2. j the primary names or the column names (entered as a list or List) 3. k assay

multiassayexperiment[i = rownames, j = primary or colnames, k = assay]

Subsetting operations always return another MultiAssayExperiment. For example, the following will return any rows named "MAPK14" or "IGFBP2", and remove any assays where no rows match:

miniACC[c("MAPK14", "IGFBP2"), , ]

The following will keep only patients of pathological stage iv, and all their associated assays:

stg4 <- miniACC$pathologic_stage == "stage iv"
# remove NA values from vector
miniACC[, stg4 & !is.na(stg4), ]

And the following will keep only the RNA-seq dataset, and only patients for which this assay is available:

miniACC[, , "RNASeq2GeneNorm"]

Subsetting by genomic ranges

If any ExperimentList objects have features represented by genomic ranges (e.g. RangedSummarizedExperiment, RaggedExperiment), then a GRanges object in the first subsetting position will subset these objects as in GenomicRanges::findOverlaps().

Double bracket [[

The "double bracket" method ([[) is a convenience function for extracting a single element of the MultiAssayExperiment ExperimentList. It avoids the use of experiments(mae)[[1L]]. For example, both of the following extract the ExpressionSet object containing RNA-seq data:

miniACC[[1L]]  #or equivalently, miniACC[["RNASeq2GeneNorm"]]

Patients with complete data

complete.cases() shows which patients have complete data for all assays:

summary(complete.cases(miniACC))

The above logical vector could be used for patient subsetting. More simply, intersectColumns() will select complete cases and rearrange each ExperimentList element so its columns correspond exactly to rows of colData in the same order:

accmatched = intersectColumns(miniACC)

Note, the column names of the assays in accmatched are not the same because of assay-specific identifiers, but they have been automatically re-arranged to correspond to the same patients. In these TCGA assays, the first three - delimited positions correspond to patient, ie the first patient is TCGA-OR-A5J2:

colnames(accmatched)

Row names that are common across assays

intersectRows() keeps only rows that are common to each assay, and aligns them in identical order. For example, to keep only genes where data are available for RNA-seq, GISTIC copy number, and somatic mutations:

accmatched2 <- intersectRows(miniACC[, , c("RNASeq2GeneNorm",
                                           "gistict",
                                           "Mutations")])
rownames(accmatched2)

back to top

Extraction

assay and assays

The assay and assays methods follow SummarizedExperiment convention. The assay (singular) method will extract the first element of the ExperimentList and will return a matrix.

class(assay(miniACC))

The assays (plural) method will return a SimpleList of the data with each element being a matrix.

assays(miniACC)

Key point: * Whereas the [[ returned an assay as its original class, assay() and assays() convert the assay data to matrix form.

back to top

Summary of slots and accessors

Slot in the MultiAssayExperiment can be accessed or set using their accessor functions:

| Slot | Accessor | |------|----------| | ExperimentList | experiments()| | colData | colData() and $ * | | sampleMap | sampleMap() | | metadata | metadata() |

* The $ operator on a MultiAssayExperiment returns a single column of the colData.

Transformation / reshaping

The longFormat or wideFormat functions will "reshape" and combine experiments with each other and with colData into one DataFrame. These functions provide compatibility with most of the common R/Bioconductor functions for regression, machine learning, and visualization.

longFormat

In long format a single column provides all assay results, with additional optional colData columns whose values are repeated as necessary. Here assay is the name of the ExperimentList element, primary is the patient identifier (rowname of colData), rowname is the assay rowname (in this case genes), colname is the assay-specific identifier (column name), value is the numeric measurement (gene expression, copy number, presence of a non-silent mutation, etc), and following these are the vital_status and days_to_death colData columns that have been added:

longFormat(miniACC[c("TP53", "CTNNB1"), , ],
           colDataCols = c("vital_status", "days_to_death"))

wideFormat

In wide format, each feature from each assay goes in a separate column, with one row per primary identifier (patient). Here, each variable becomes a new column:

wideFormat(miniACC[c("TP53", "CTNNB1"), , ],
           colDataCols = c("vital_status", "days_to_death"))

MultiAssayExperiment class construction and concatenation

MultiAssayExperiment constructor function

The MultiAssayExperiment constructor function can take three arguments:

1. experiments - An ExperimentList or list of data
2. colData - A DataFrame describing the patients (or cell lines, or other biological units)
3. sampleMap - A DataFrame of assay, primary, and colname identifiers

The miniACC object can be reconstructed as follows:

MultiAssayExperiment(experiments=experiments(miniACC),
    colData=colData(miniACC),
    sampleMap=sampleMap(miniACC),
    metadata=metadata(miniACC))

prepMultiAssay - Constructor function helper

The prepMultiAssay function allows the user to diagnose typical problems when creating a MultiAssayExperiment object. See ?prepMultiAssay for more details.

c - concatenate to MultiAssayExperiment

The c function allows the user to concatenate an additional experiment to an existing MultiAssayExperiment. The optional sampleMap argument allows concatenating an assay whose column names do not match the row names of colData. For convenience, the mapFrom argument allows the user to map from a particular experiment provided that the order of the colnames is in the same. A warning will be issued to make the user aware of this assumption. For example, to concatenate a matrix of log2-transformed RNA-seq results:

miniACC2 <- c(miniACC,
    log2rnaseq = log2(assays(miniACC)$RNASeq2GeneNorm), mapFrom=1L)
assays(miniACC2)

back to top

Examples

UpsetR "Venn" diagram

We see that 43 samples have all 5 assays, 32 are missing reverse-phase protein (RPPAArray), 2 are missing Mutations, 1 is missing gistict, 12 have only mutations and gistict, etc:

library(UpSetR)
upsetSamples(miniACC)

Kaplan-meier plot stratified by a clinical variable

The colData can provide clinical data for things like a Kaplan-Meier plot for overall survival stratified by nodal stage. To simplify things, first add a "y" column to the colData, containing the Surv object for survival analysis:

Note: survfit method does not work well with DataFrame. To bypass the error, here we covert colData to a data.frame.

library(survival)
library(survminer)

coldat <- as.data.frame(colData(miniACC))
coldat$y <- Surv(miniACC$days_to_death, miniACC$vital_status)
colData(miniACC) <- DataFrame(coldat)

And remove any patients missing overall survival information:

miniACC <- miniACC[, complete.cases(coldat$y), ]
coldat <- as(colData(miniACC), "data.frame")

fit <- survfit(y ~ pathology_N_stage, data = coldat)
ggsurvplot(fit, data = coldat, risk.table = TRUE)

Multivariate Cox regression including RNA-seq, copy number, and pathology

Choose the EZH2 gene for demonstration. This subsetting will drop assays with no row named EZH2:

wideacc <- wideFormat(miniACC["EZH2", , ],
    colDataCols = c("vital_status", "days_to_death", "pathology_N_stage"))
wideacc$y <- Surv(wideacc$days_to_death, wideacc$vital_status)
head(wideacc)

Perform a multivariate Cox regression with EZH2 copy number (gistict), log2-transformed EZH2 expression (RNASeq2GeneNorm), and nodal status (pathology_N_stage) as predictors:

coxph(Surv(days_to_death, vital_status) ~ gistict_EZH2 +
          log2(RNASeq2GeneNorm_EZH2) + pathology_N_stage,  data=wideacc)

We see that EZH2 expression is significantly associated with overal survival (p < 0.001), but EZH2 copy number and nodal status are not. This analysis could easily be extended to the whole genome for discovery of prognostic features by repeated univariate regressions over columns, penalized multivariate regression, etc.

For further detail, see the main MultiAssayExperiment vignette.

Session info

sessionInfo()

back to top

waldronlab/MultiAssayExperiment documentation built on Nov. 4, 2024, 7:51 a.m.