Workflow for Multi-omics Analysis with MultiAssayExperiment and curatedTCGAData

Instructor names and contact information

Workshop Description

This workshop demonstrates data management and analyses of multiple assays associated with a single set of biological specimens, using the MultiAssayExperiment data class and methods. It introduces the facilities made available by companion packages, curatedTCGAData and TCGAutils, to work with TCGA data. Built using the MultiAssayExperiment class, these packages make the management of multiple assays easier and more efficient. The workshop also covers relevant data classes such as SummarizedExperiment and RangedSummarizedExperiment, which provide efficient and powerful operations for representation of copy number, mutation, variant, and expression data that are represented by different genomic ranges for each specimen.

Pre-requisites

  • Basic knowledge of R syntax
  • Familiarity with the GRanges and SummarizedExperiment classes
  • Familiarity with ’omics data types including copy number and gene expression

Workshop Participation

Participants will have a chance to build a MultiAssayExperiment object from scratch, and will also work with more complex objects provided by the curatedTCGAData package.

R/Bioconductor packages used

First, please make sure you have the latest Bioconductor release version 3.9. See Bioconductor installation instructions:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

## Bioconductor 3.9 Stable Release
BiocManager::install(version = "3.9")
library(MultiAssayExperiment)
library(GenomicRanges)
library(RaggedExperiment)
library(curatedTCGAData)
library(GenomicDataCommons)
library(SummarizedExperiment)
library(SingleCellExperiment)
library(TCGAutils)
library(UpSetR)
library(mirbase.db)
library(AnnotationFilter)
library(EnsDb.Hsapiens.v86)
library(survival)
library(survminer)
library(pheatmap)
library(readxl)
library(kableExtra)

Time outline

1h 40m total

Activity Time
Overview of key data classes and packages 30m
Building a MultiAssayExperiment from scratch 10m
curatedTCGAData multi-assay dataset 20 min
TCGAutils functionality 10m
Subsetting and reshaping multi-assay data 20m
Plotting, correlation, analysis and other examples 10m

Workshop goals and objectives

Learning goals

  • overview common data classes available in Bioconductor
  • identify appropriate data structures for different ’omics data types
  • identify advantages of providing integrative data in an analysis-ready platform
  • gain familiarity with available functionality in TCGAutils for the management and coordination of TGCA data

Learning objectives

  • use curatedTCGAData to create custom on-the-fly TCGA MultiAssayExperiment objects
  • create a MultiAssayExperiment for TCGA or other multi’omics data
  • perform subsetting, reshaping, growing, and extraction of a MultiAssayExperiment
  • link MultiAssayExperiment data with packages for differential expression, machine learning, and plotting

Overview of key packages and data classes

Packages

curatedTCGAData

Many tools exist for accessing and downloading The Cancer Genome Atlas (TCGA) data. These include but are not limited to RTCGAToolbox, GenomicDataCommons (package and website), TCGAbiolinks, cBioPortal, and Broad GDAC Firehose. These tools encompass a spectrum of strengths in ease-of-use, integration, and completeness of data. Few tools provide an integrative and user-friendly representation of TCGA data in a widely used analysis platform such as Bioconductor.

The curatedTCGAData experiment data package provides on-the-fly construction of TCGA datasets for 33 different cancer types from the Broad GDAC Firehose pipeline. Generally, it provides data using build hg19. curatedTCGAData facilitates access and integration of TCGA data by providing multi-’omics data objects using the MultiAssayExperiment data class. Where other platforms provide fragmented datasets, curatedTCGAData ensures that all data provided is matched and accounted for within the phenotypic metadata.

A list of available cancer types can be obtained from TCGAutils::diseaseCodes:

data("diseaseCodes")
knitr::kable(diseaseCodes, align = "l", escape = FALSE,
    caption = "List of available cancer types from curatedTCGAData") %>%
    kableExtra::kable_styling(
        bootstrap_options = c("hover", "striped", "responsive"),
        full_width = FALSE
    )
Table 1: List of available cancer types from curatedTCGAData
Study.Abbreviation Available SubtypeData Study.Name
ACC Yes Yes Adrenocortical carcinoma
BLCA Yes Yes Bladder Urothelial Carcinoma
BRCA Yes Yes Breast invasive carcinoma
CESC Yes No Cervical squamous cell carcinoma and endocervical adenocarcinoma
CHOL Yes No Cholangiocarcinoma
CNTL No No Controls
COAD Yes Yes Colon adenocarcinoma
DLBC Yes No Lymphoid Neoplasm Diffuse Large B-cell Lymphoma
ESCA Yes No Esophageal carcinoma
FPPP No No FFPE Pilot Phase II
GBM Yes Yes Glioblastoma multiforme
HNSC Yes Yes Head and Neck squamous cell carcinoma
KICH Yes Yes Kidney Chromophobe
KIRC Yes Yes Kidney renal clear cell carcinoma
KIRP Yes Yes Kidney renal papillary cell carcinoma
LAML Yes Yes Acute Myeloid Leukemia
LCML No No Chronic Myelogenous Leukemia
LGG Yes Yes Brain Lower Grade Glioma
LIHC Yes No Liver hepatocellular carcinoma
LUAD Yes Yes Lung adenocarcinoma
LUSC Yes Yes Lung squamous cell carcinoma
MESO Yes No Mesothelioma
MISC No No Miscellaneous
OV Yes Yes Ovarian serous cystadenocarcinoma
PAAD Yes No Pancreatic adenocarcinoma
PCPG Yes No Pheochromocytoma and Paraganglioma
PRAD Yes Yes Prostate adenocarcinoma
READ Yes No Rectum adenocarcinoma
SARC Yes No Sarcoma
SKCM Yes Yes Skin Cutaneous Melanoma
STAD Yes Yes Stomach adenocarcinoma
TGCT Yes No Testicular Germ Cell Tumors
THCA Yes Yes Thyroid carcinoma
THYM Yes No Thymoma
UCEC Yes Yes Uterine Corpus Endometrial Carcinoma
UCS Yes No Uterine Carcinosarcoma
UVM Yes No Uveal Melanoma

A descriptive table of available ’omics types in curatedTCGAData (thanks to Ludwig G. @lgeistlinger):

## from dataTypesHTML.R
knitr::kable(dataTypes, align = "l", escape = FALSE,
    caption = "Descriptions of data types available by Bioconductor data class") %>%
    kableExtra::kable_styling(bootstrap_options = c("hover", "striped", "responsive"),
        full_width = FALSE) %>% group_rows(index = setNames(iwidth, groups)) %>%
    footnote(symbol = footnote)
Table 2: Descriptions of data types available by Bioconductor data class
ExperimentList data types Description
SummarizedExperiment*
RNASeqGene RSEM TPM gene expression values
RNASeq2GeneNorm Upper quartile normalized RSEM TPM gene expression values
miRNAArray Probe-level miRNA expression values
miRNASeqGene Gene-level log2 RPM miRNA expression values
mRNAArray Unified gene-level mRNA expression values
mRNAArray_huex Gene-level mRNA expression values from Affymetrix Human Exon Array
mRNAArray_TX_g4502a Gene-level mRNA expression values from Agilent 244K Array
mRNAArray_TX_ht_hg_u133a Gene-level mRNA expression values from Affymetrix Human Genome U133 Array
GISTIC_AllByGene Gene-level GISTIC2 copy number values
GISTIC_ThresholdedByGene Gene-level GISTIC2 thresholded discrete copy number values
RPPAArray Reverse Phase Protein Array normalized protein expression values
RangedSummarizedExperiment
GISTIC_Peaks GISTIC2 thresholded discrete copy number values in recurrent peak regions
SummarizedExperiment with HDF5Array DelayedMatrix
Methylation_methyl27 Probe-level methylation beta values from Illumina HumanMethylation 27K BeadChip
Methylation_methyl450 Probe-level methylation beta values from Infinium HumanMethylation 450K BeadChip
RaggedExperiment
CNASNP Segmented somatic Copy Number Alteration calls from SNP array
CNVSNP Segmented germline Copy Number Variant calls from SNP Array
CNASeq Segmented somatic Copy Number Alteration calls from low pass DNA Sequencing
Mutation* Somatic mutations calls
CNACGH_CGH_hg_244a Segmented somatic Copy Number Alteration calls from CGH Agilent Microarray 244A
CNACGH_CGH_hg_415k_g4124a Segmented somatic Copy Number Alteration calls from CGH Agilent Microarray 415K
* All can be converted to RangedSummarizedExperiment (except RPPAArray) with TCGAutils

TCGAutils

TCGAutils is a companion package that enhances curatedTCGAData by allowing additional exploration and manipulation of samples and metadata in rows and columns.

Available operations in TCGAutils and MultiAssayExperiment enable user-friendly operations for subsetting, separating, converting, and reshaping of sample and feature TCGA data.

TCGAutils was developed specifically for TCGA data and for curatedTCGAData products. It provides convenience / helper functions in three major areas:

  1. conversion / summarization of row annotations to genomic ranges
  2. identification and separation of samples
  3. translation and interpretation of TCGA identifiers

along with additional reference data sets to explore TCGA data.

## from TCGAutils_cheatHTML.R script
knitr::kable(tabInfo, align = "l", escape = FALSE,
    caption = "Summary of available functionality in TCGAutils") %>%
    kable_styling(bootstrap_options = c("hover", "striped", "responsive"),
        full_width = FALSE) %>% group_rows(index = setNames(iwidth, groups))
Table 3: Summary of available functionality in TCGAutils
Category and Function Description
MultiAssayExperiment helpers
generateMap Automatically generate a sampleMap structure from assays
imputeAssay Impute values in numerical assays based on KNN
mergeColData Add additional data to the colData slot of MAE
qreduceTCGA Convert RaggedExperiment assays to RSE based on heuristics
symbolsToRanges Convert gene symbols to genomic ranges using org.db
mirToRanges Convert microRNA sequences to genomic ranges with mirbase.db
simplifyTCGA Use qreduceTCGA, symbolsToRanges, and mirToRanges in succession
curatedTCGAData helpers
getSubtypeMap Obtain the available subtype information
getClinicalNames Get a list of clinical variable names for all cancer types
splitAssays Separate assays based on sample data found in barcodes
sampleTables Get a list of samples in each assay in a MAE
TCGA Identifiers
TCGAbarcode Chop TCGA barcode into sections
TCGAbiospec Get a table of information extracted from a vector of barcodes
TCGAsampleSelect Indicate which barcodes belong to a specific sample type
UUIDtoBarcode Translate universal identifiers to TCGA barcodes
UUIDtoUUID Translate between case and file universal identifiers
barcodeToUUID Translate TCGA barcodes to universal identifiers
filenameToBarcode Obtain TCGA barcodes from a vector of TCGA file names
Flat to Bioconductor classes
makeGRangesListFromCopyNumber Create a GRangesList from a copy number data.frame
makeGRangesListFromExonFiles Obtain a GRangesList object from a list of individual exon files
makeSummarizedExperimentFromGISTIC Create a SummarizedExperiment object from a Firehose GISTIC RTCGAToolbox object
Genome Builds
translateBuild Translate build version name between UCSC and NCBI
extractBuild Find build in string pattern such as a file name
uniformBuilds Homogenize a vector of builds based on a threshold for the alternative build name
Reference data
diseaseCodes Get a table of TCGA cancer codes and subtype availability
sampleTypes Get a table of sample type codes and their definition
clinicalNames Obtain a CharacterList of common variable names for each TCGA disease code
getFilename Obtain a file name string for the relevant data query
Miscellaneous
findGRangesCols Find the minimum necessary variable names for conversion to GRanges

To better understand how it all fits together, this schematic shows the relationship among all as part of the curatedTCGAData pipeline.

Schematic of curatedTCGAData Pipeline

Figure 1: Schematic of curatedTCGAData Pipeline

Classes

This section summarizes three fundamental data classes for the representation of multi-omics experiments.

(Ranged)SummarizedExperiment

A matrix-like container where rows represent features of interest and columns represent samples. The objects contain one or more assays, each represented by a matrix-like object of numeric or other mode.

Figure 2: A matrix-like container where rows represent features of interest and columns represent samples. The objects contain one or more assays, each represented by a matrix-like object of numeric or other mode.

SummarizedExperiment is the most important Bioconductor class for matrix-like experimental data, including from RNA sequencing and microarray experiments. It can store multiple experimental data matrices of identical dimensions, with associated metadata on the rows/genes/transcripts/other measurements (rowData), column/sample phenotype or clinical data (colData), and the overall experiment (metadata). The derivative class RangedSummarizedExperiment associates a GRanges or GRangesList vector with the rows. These classes supersede the use of ExpressionSet. Note that many other classes for experimental data are actually derived from SummarizedExperiment; for example, the SingleCellExperiment class for single-cell RNA sequencing experiments extends RangedSummarizedExperiment, which in turn extends SummarizedExperiment:

library(SingleCellExperiment)
extends("SingleCellExperiment")
## [1] "SingleCellExperiment"       "RangedSummarizedExperiment"
## [3] "SummarizedExperiment"       "Vector"                    
## [5] "Annotated"                  "vector_OR_Vector"

Thus, although SingleCellExperiment provides additional methods over RangedSummarizedExperiment, it also inherits all the methods of SummarizedExperiment and RangedSummarizedExperiment, so everything you learn about SummarizedExperiment will be applicable to SingleCellExperiment.

RaggedExperiment

RaggedExperiment is a flexible data representation for segmented copy number, somatic mutations such as represented in .vcf files, and other ragged array schema for genomic location data. Like the GRangesList class from GenomicRanges, RaggedExperiment can be used to represent differing genomic ranges on each of a set of samples. In fact, RaggedExperiment contains a GRangesList:

showClass("RaggedExperiment")
## Class "RaggedExperiment" [package "RaggedExperiment"]
## 
## Slots:
##                                                       
## Name:       assays      rowidx      colidx    metadata
## Class: GRangesList     integer     integer        list
## 
## Extends: "Annotated"

However, RaggedExperiment provides a flexible set of Assay methods to support transformation of such data to matrix format.

RaggedExperiment object schematic. Rows and columns represent genomic ranges and samples, respectively. Assay operations can be performed with (from left to right) compactAssay, qreduceAssay, and sparseAssay.

Figure 3: RaggedExperiment object schematic. Rows and columns represent genomic ranges and samples, respectively. Assay operations can be performed with (from left to right) compactAssay, qreduceAssay, and sparseAssay.

MultiAssayExperiment

MultiAssayExperiment is an integrative container for coordinating multi-omics experiment data on a set of biological specimens. As much as possible, its methods adopt the same vocabulary as SummarizedExperiment. A MultiAssayExperiment can contain any number of assays with different representations. Assays may be ID-based, where measurements are indexed identifiers of genes, microRNA, proteins, microbes, etc. Alternatively, assays may be range-based, where measurements correspond to genomic ranges that can be represented as GRanges objects, such as gene expression or copy number. For ID-based assays, there is no requirement that the same IDs be present for different experiments. For range-based assays, there is also no requirement that the same ranges be present for different experiments; furthermore, it is possible for different samples within an experiment to be represented by different ranges. The following data classes have been tested to work as elements of a MultiAssayExperiment:

  1. matrix: the most basic class for ID-based datasets, could be used for example for gene expression summarized per-gene, microRNA, metabolomics, or microbiome data.
  2. SummarizedExperiment and derived methods: described above, could be used for miRNA, gene expression, proteomics, or any matrix-like data where measurements are represented by IDs.
  3. RangedSummarizedExperiment: described above, could be used for gene expression, methylation, or other data types referring to genomic positions.
  4. ExpressionSet: Another rich representation for ID-based datasets, supported only for legacy reasons
  5. RaggedExperiment: described above, for non-rectangular (ragged) ranged-based datasets such as segmented copy number, where segmentation of copy number alterations occurs and different genomic locations in each sample.
  6. RangedVcfStack: For VCF archives broken up by chromosome (see VcfStack class defined in the GenomicFiles package)
  7. DelayedMatrix: An on-disk representation of matrix-like objects for large datasets. It reduces memory usage and optimizes performance with delayed operations. This class is part of the DelayedArray package.

Note that any data class extending these classes, and in fact any data class supporting row and column names and subsetting can be used as an element of a MultiAssayExperiment.

MultiAssayExperiment object schematic. colData provides data about the patients, cell lines, or other biological units, with one row per unit and one column per variable. The experiments are a list of assay datasets of arbitrary class.  The sampleMap relates each column (observation) in ExperimentList to exactly one row (biological unit) in colData; however, one row of colData may map to zero, one, or more columns per assay, allowing for missing and replicate assays. sampleMap allows for per-assay sample naming conventions. Metadata can be used to store information in arbitrary format about the MultiAssayExperiment. Green stripes indicate a mapping of one subject to multiple observations across experiments.

Figure 4: MultiAssayExperiment object schematic. colData provides data about the patients, cell lines, or other biological units, with one row per unit and one column per variable. The experiments are a list of assay datasets of arbitrary class. The sampleMap relates each column (observation) in ExperimentList to exactly one row (biological unit) in colData; however, one row of colData may map to zero, one, or more columns per assay, allowing for missing and replicate assays. sampleMap allows for per-assay sample naming conventions. Metadata can be used to store information in arbitrary format about the MultiAssayExperiment. Green stripes indicate a mapping of one subject to multiple observations across experiments.

Working with MultiAssayExperiment

API cheat sheet

## New names:
## * `` -> ...2

The MultiAssayExperiment miniACC demo object

Get started by trying out MultiAssayExperiment using a subset of the TCGA adrenocortical carcinoma (ACC) dataset provided with the package. 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
## A MultiAssayExperiment object of 5 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 5: 
##  [1] RNASeq2GeneNorm: SummarizedExperiment with 198 rows and 79 columns 
##  [2] gistict: SummarizedExperiment with 198 rows and 90 columns 
##  [3] RPPAArray: SummarizedExperiment with 33 rows and 46 columns 
##  [4] Mutations: matrix with 97 rows and 90 columns 
##  [5] miRNASeqGene: SummarizedExperiment with 471 rows and 80 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

colData - information biological units

This slot is 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]
## DataFrame with 4 rows and 4 columns
##                 patientID years_to_birth vital_status days_to_death
##               <character>      <integer>    <integer>     <integer>
## TCGA-OR-A5J1 TCGA-OR-A5J1             58            1          1355
## TCGA-OR-A5J2 TCGA-OR-A5J2             44            1          1677
## TCGA-OR-A5J3 TCGA-OR-A5J3             23            0            NA
## TCGA-OR-A5J4 TCGA-OR-A5J4             23            1           423
table(miniACC$race)
## 
##                     asian black or african american 
##                         2                         1 
##                     white 
##                        78

Key points about the colData:

  • Each row maps to zero or more observations in each experiment in the ExperimentList, below.
  • One row per biological unit
    • MultiAssayExperiment supports both missing observations and replicate observations, ie one row of colData can map to 0, 1, or more columns of any of the experimental data matrices.
    • therefore you could treat replicate observations as one or multiple rows of colData, and this will result in different behaviors of functions you will learn later like subsetting, duplicated(), and wideFormat().
    • multiple time points, or distinct biological replicates, should probably be separate rows of the colData.

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)
## ExperimentList class object of length 5: 
##  [1] RNASeq2GeneNorm: SummarizedExperiment with 198 rows and 79 columns 
##  [2] gistict: SummarizedExperiment with 198 rows and 90 columns 
##  [3] RPPAArray: SummarizedExperiment with 33 rows and 46 columns 
##  [4] Mutations: matrix with 97 rows and 90 columns 
##  [5] miRNASeqGene: SummarizedExperiment with 471 rows and 80 columns

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)
## DataFrame with 385 rows and 3 columns
##               assay      primary                      colname
##            <factor>  <character>                  <character>
## 1   RNASeq2GeneNorm TCGA-OR-A5J1 TCGA-OR-A5J1-01A-11R-A29S-07
## 2   RNASeq2GeneNorm TCGA-OR-A5J2 TCGA-OR-A5J2-01A-11R-A29S-07
## 3   RNASeq2GeneNorm TCGA-OR-A5J3 TCGA-OR-A5J3-01A-11R-A29S-07
## 4   RNASeq2GeneNorm TCGA-OR-A5J5 TCGA-OR-A5J5-01A-11R-A29S-07
## 5   RNASeq2GeneNorm TCGA-OR-A5J6 TCGA-OR-A5J6-01A-31R-A29S-07
## ...             ...          ...                          ...
## 381    miRNASeqGene TCGA-PA-A5YG TCGA-PA-A5YG-01A-11R-A29W-13
## 382    miRNASeqGene TCGA-PK-A5H8 TCGA-PK-A5H8-01A-11R-A29W-13
## 383    miRNASeqGene TCGA-PK-A5H9 TCGA-PK-A5H9-01A-11R-A29W-13
## 384    miRNASeqGene TCGA-PK-A5HA TCGA-PK-A5HA-01A-11R-A29W-13
## 385    miRNASeqGene TCGA-PK-A5HB TCGA-PK-A5HB-01A-11R-A29W-13

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)
## $title
## [1] "Comprehensive Pan-Genomic Characterization of Adrenocortical Carcinoma"
## 
## $PMID
## [1] "27165744"
## 
## $sourceURL
## [1] "http://s3.amazonaws.com/multiassayexperiments/accMAEO.rds"
## 
## $RPPAfeatureDataURL
## [1] "http://genomeportal.stanford.edu/pan-tcga/show_target_selection_file?filename=Allprotein.txt"
## 
## $colDataExtrasURL
## [1] "http://www.cell.com/cms/attachment/2062093088/2063584534/mmc3.xlsx"

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.

MultiAssayExperiment 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"), , ]
## harmonizing input:
##   removing 170 sampleMap rows not in names(experiments)

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

miniACC[, miniACC$pathologic_stage == "stage iv", ]

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

miniACC[, , "RNASeq2GeneNorm"]
## harmonizing input:
##   removing 306 sampleMap rows not in names(experiments)
##   removing 13 colData rownames not in sampleMap 'primary'

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(). Any non-ranged ExperimentList element will be subset to zero rows.

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"]]
## class: SummarizedExperiment 
## dim: 198 79 
## metadata(3): experimentData annotation protocolData
## assays(1): exprs
## rownames(198): DIRAS3 MAPK14 ... SQSTM1 KCNJ13
## rowData names(0):
## colnames(79): TCGA-OR-A5J1-01A-11R-A29S-07
##   TCGA-OR-A5J2-01A-11R-A29S-07 ... TCGA-PK-A5HA-01A-11R-A29S-07
##   TCGA-PK-A5HB-01A-11R-A29S-07
## colData names(0):

Complete cases

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

summary(complete.cases(miniACC))
##    Mode   FALSE    TRUE 
## logical      49      43

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)
## CharacterList of length 5
## [["RNASeq2GeneNorm"]] TCGA-OR-A5J2-01A-11R-A29S-07 ...
## [["gistict"]] TCGA-OR-A5J2-01A-11D-A29H-01 ...
## [["RPPAArray"]] TCGA-OR-A5J2-01A-21-A39K-20 ...
## [["Mutations"]] TCGA-OR-A5J2-01A-11D-A29I-10 ...
## [["miRNASeqGene"]] TCGA-OR-A5J2-01A-11R-A29W-13 ...

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")])
## harmonizing input:
##   removing 126 sampleMap rows not in names(experiments)
rownames(accmatched2)
## CharacterList of length 3
## [["RNASeq2GeneNorm"]] DIRAS3 G6PD KDR ERBB3 ... RET CDKN2A MACC1 CHGA
## [["gistict"]] DIRAS3 G6PD KDR ERBB3 AKT1S1 ... PREX1 RET CDKN2A MACC1 CHGA
## [["Mutations"]] DIRAS3 G6PD KDR ERBB3 AKT1S1 ... RET CDKN2A MACC1 CHGA

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))
## [1] "matrix"

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

assays(miniACC)
## List of length 5
## names(5): RNASeq2GeneNorm gistict RPPAArray Mutations miRNASeqGene

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"))
## harmonizing input:
##   removing 126 sampleMap rows not in names(experiments)
## DataFrame with 518 rows and 7 columns
##               assay      primary     rowname                      colname
##         <character>  <character> <character>                  <character>
## 1   RNASeq2GeneNorm TCGA-OR-A5J1        TP53 TCGA-OR-A5J1-01A-11R-A29S-07
## 2   RNASeq2GeneNorm TCGA-OR-A5J1      CTNNB1 TCGA-OR-A5J1-01A-11R-A29S-07
## 3   RNASeq2GeneNorm TCGA-OR-A5J2        TP53 TCGA-OR-A5J2-01A-11R-A29S-07
## 4   RNASeq2GeneNorm TCGA-OR-A5J2      CTNNB1 TCGA-OR-A5J2-01A-11R-A29S-07
## 5   RNASeq2GeneNorm TCGA-OR-A5J3        TP53 TCGA-OR-A5J3-01A-11R-A29S-07
## ...             ...          ...         ...                          ...
## 514       Mutations TCGA-PK-A5HA      CTNNB1 TCGA-PK-A5HA-01A-11D-A29I-10
## 515       Mutations TCGA-PK-A5HB        TP53 TCGA-PK-A5HB-01A-11D-A29I-10
## 516       Mutations TCGA-PK-A5HB      CTNNB1 TCGA-PK-A5HB-01A-11D-A29I-10
## 517       Mutations TCGA-PK-A5HC        TP53 TCGA-PK-A5HC-01A-11D-A30A-10
## 518       Mutations TCGA-PK-A5HC      CTNNB1 TCGA-PK-A5HC-01A-11D-A30A-10
##          value vital_status days_to_death
##      <numeric>    <integer>     <integer>
## 1     563.4006            1          1355
## 2    5634.4669            1          1355
## 3     165.4811            1          1677
## 4   62658.3913            1          1677
## 5     956.3028            0            NA
## ...        ...          ...           ...
## 514          0            0            NA
## 515          0            0            NA
## 516          0            0            NA
## 517          0            0            NA
## 518          0            0            NA

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"))
## harmonizing input:
##   removing 126 sampleMap rows not in names(experiments)
## DataFrame with 92 rows and 9 columns
##          primary vital_status days_to_death RNASeq2GeneNorm_CTNNB1
##      <character>    <integer>     <integer>              <numeric>
## 1   TCGA-OR-A5J1            1          1355              5634.4669
## 2   TCGA-OR-A5J2            1          1677             62658.3913
## 3   TCGA-OR-A5J3            0            NA              6337.4256
## 4   TCGA-OR-A5J4            1           423                     NA
## 5   TCGA-OR-A5J5            1           365               5979.055
## ...          ...          ...           ...                    ...
## 88  TCGA-PK-A5H9            0            NA              5258.9863
## 89  TCGA-PK-A5HA            0            NA              8120.1654
## 90  TCGA-PK-A5HB            0            NA              5257.8148
## 91  TCGA-PK-A5HC            0            NA                     NA
## 92  TCGA-P6-A5OG            1           383              6390.0997
##     RNASeq2GeneNorm_TP53 gistict_CTNNB1 gistict_TP53 Mutations_CTNNB1
##                <numeric>      <numeric>    <numeric>        <numeric>
## 1               563.4006              0            0                0
## 2               165.4811              1            0                1
## 3               956.3028              0            0                0
## 4                     NA              0            1                0
## 5              1169.6359              0            0                0
## ...                  ...            ...          ...              ...
## 88              890.8663              0            0                0
## 89              683.5722              0           -1                0
## 90              237.3697             -1           -1                0
## 91                    NA              1            1                0
## 92              815.3446              1           -1               NA
##     Mutations_TP53
##          <numeric>
## 1                0
## 2                1
## 3                0
## 4                0
## 5                0
## ...            ...
## 88               0
## 89               0
## 90               0
## 91               0
## 92              NA

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))
## A MultiAssayExperiment object of 5 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 5: 
##  [1] RNASeq2GeneNorm: SummarizedExperiment with 198 rows and 79 columns 
##  [2] gistict: SummarizedExperiment with 198 rows and 90 columns 
##  [3] RPPAArray: SummarizedExperiment with 33 rows and 46 columns 
##  [4] Mutations: matrix with 97 rows and 90 columns 
##  [5] miRNASeqGene: SummarizedExperiment with 471 rows and 80 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

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)
## Warning: Assuming column order in the data provided 
##  matches the order in 'mapFrom' experiment(s) colnames
experiments(miniACC2)
## ExperimentList class object of length 6: 
##  [1] RNASeq2GeneNorm: SummarizedExperiment with 198 rows and 79 columns 
##  [2] gistict: SummarizedExperiment with 198 rows and 90 columns 
##  [3] RPPAArray: SummarizedExperiment with 33 rows and 46 columns 
##  [4] Mutations: matrix with 97 rows and 90 columns 
##  [5] miRNASeqGene: SummarizedExperiment with 471 rows and 80 columns 
##  [6] log2rnaseq: matrix with 198 rows and 79 columns

back to top

Building a MultiAssayExperiment from scratch

To start from scratch building your own MultiAssayExperiment, see the package Coordinating Analysis of Multi-Assay Experiments vignette. The package cheat sheet is also helpful.

If anything is unclear, please ask a question at https://support.bioconductor.org/ or create an issue on the MultiAssayExperiment issue tracker.

The Cancer Genome Atlas (TCGA) as MultiAssayExperiment objects

Most unrestricted TCGA data are available as MultiAssayExperiment objects from the curatedTCGAData package. This represents a lot of harmonization!

library(curatedTCGAData)
curatedTCGAData("ACC")
##                                    Title DispatchClass
## 1                    ACC_CNASNP-20160128           Rda
## 2                    ACC_CNVSNP-20160128           Rda
## 4          ACC_GISTIC_AllByGene-20160128           Rda
## 5              ACC_GISTIC_Peaks-20160128           Rda
## 6  ACC_GISTIC_ThresholdedByGene-20160128           Rda
## 8        ACC_Methylation-20160128_assays        H5File
## 9            ACC_Methylation-20160128_se           Rds
## 10             ACC_miRNASeqGene-20160128           Rda
## 11                 ACC_Mutation-20160128           Rda
## 12          ACC_RNASeq2GeneNorm-20160128           Rda
## 13                ACC_RPPAArray-20160128           Rda
suppressMessages({
    acc <- curatedTCGAData("ACC",
        assays = c("miRNASeqGene", "RPPAArray", "Mutation", "RNASeq2GeneNorm", "CNVSNP"),
        dry.run = FALSE)
})
acc
## A MultiAssayExperiment object of 5 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 5: 
##  [1] ACC_CNVSNP-20160128: RaggedExperiment with 21052 rows and 180 columns 
##  [2] ACC_miRNASeqGene-20160128: SummarizedExperiment with 1046 rows and 80 columns 
##  [3] ACC_Mutation-20160128: RaggedExperiment with 20166 rows and 90 columns 
##  [4] ACC_RNASeq2GeneNorm-20160128: SummarizedExperiment with 20501 rows and 79 columns 
##  [5] ACC_RPPAArray-20160128: SummarizedExperiment with 192 rows and 46 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

These objects contain most unrestricted TCGA assay and clinical / pathological data, as well as material curated from the supplements of published TCGA primary papers at the end of the colData columns:

dim(colData(acc))
## [1]  92 822
tail(colnames(colData(acc)), 10)
##  [1] "MethyLevel"       "miRNA.cluster"    "SCNA.cluster"    
##  [4] "protein.cluster"  "COC"              "OncoSign"        
##  [7] "purity"           "ploidy"           "genome_doublings"
## [10] "ADS"

The TCGAutils package provides additional helper functions, see below.

Utilities for TCGA

Aside from the available reshaping functions already included in the MultiAssayExperiment package, the TCGAutils package provides additional helper functions for working with TCGA data.

“Simplification” of curatedTCGAData objects

A number of helper functions are available for managing datasets from curatedTCGAData. These include:

  • Conversions of SummarizedExperiment to RangedSummarizedExperiment based on TxDb.Hsapiens.UCSC.hg19.knownGene for:
    • mirToRanges: microRNA
    • symbolsToRanges: gene symbols
  • qreduceTCGA: convert RaggedExperiment objects to RangedSummarizedExperiment with one row per gene symbol, for:
    • segmented copy number datasets (“CNVSNP” and “CNASNP”)
    • somatic mutation datasets (“Mutation”), with a value of 1 for any non-silent mutation and a value of 0 for no mutation or silent mutation

The simplifyTCGA function combines all of the above operations to create a more managable MultiAssayExperiment object and using RangedSummarizedExperiment assays where possible.

genome(acc[[3]]) <- vapply(genome(acc[[3]]), translateBuild, character(1L))
(simpa <- TCGAutils::simplifyTCGA(acc))
## 
## 'select()' returned 1:1 mapping between keys and columns
## harmonizing input:
##   removing 270 sampleMap rows not in names(experiments)
## harmonizing input:
##   removing 80 sampleMap rows not in names(experiments)
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
## harmonizing input:
##   removing 349 sampleMap rows not in names(experiments)
## A MultiAssayExperiment object of 7 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 7: 
##  [1] ACC_RPPAArray-20160128: SummarizedExperiment with 192 rows and 46 columns 
##  [2] ACC_miRNASeqGene-20160128_ranged: SummarizedExperiment with 1002 rows and 80 columns 
##  [3] ACC_miRNASeqGene-20160128_unranged: SummarizedExperiment with 44 rows and 80 columns 
##  [4] ACC_RNASeq2GeneNorm-20160128_ranged: SummarizedExperiment with 17449 rows and 79 columns 
##  [5] ACC_RNASeq2GeneNorm-20160128_unranged: SummarizedExperiment with 3052 rows and 79 columns 
##  [6] ACC_Mutation-20160128_simplified_ranged: SummarizedExperiment with 22938 rows and 90 columns 
##  [7] ACC_CNVSNP-20160128_simplified_ranged: SummarizedExperiment with 22938 rows and 180 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

What types of samples are in the data?

Solution

The sampleTables function gives you an overview of samples in each assay:

sampleTables(acc)
## $`ACC_CNVSNP-20160128`
## 
## 01 10 11 
## 90 85  5 
## 
## $`ACC_miRNASeqGene-20160128`
## 
## 01 
## 80 
## 
## $`ACC_Mutation-20160128`
## 
## 01 
## 90 
## 
## $`ACC_RNASeq2GeneNorm-20160128`
## 
## 01 
## 79 
## 
## $`ACC_RPPAArray-20160128`
## 
## 01 
## 46
head(sampleTypes)
##   Code                                      Definition Short.Letter.Code
## 1   01                             Primary Solid Tumor                TP
## 2   02                           Recurrent Solid Tumor                TR
## 3   03 Primary Blood Derived Cancer - Peripheral Blood                TB
## 4   04    Recurrent Blood Derived Cancer - Bone Marrow              TRBM
## 5   05                        Additional - New Primary               TAP
## 6   06                                      Metastatic                TM

Curated molecular subtypes

Is there subtype data available in the MultiAssayExperiment obtained from curatedTCGAData?

Solution

The getSubtypeMap function will show actual variable names found in colData that contain subtype information. This can only be obtained from MultiAssayExperiment objects provided by curatedTCGAData.

getSubtypeMap(acc)
##          ACC_annotations     ACC_subtype
## 1             Patient_ID       patientID
## 2  histological_subtypes       Histology
## 3          mrna_subtypes         C1A/C1B
## 4          mrna_subtypes         mRNA_K4
## 5                   cimp      MethyLevel
## 6      microrna_subtypes   miRNA cluster
## 7          scna_subtypes    SCNA cluster
## 8       protein_subtypes protein cluster
## 9   integrative_subtypes             COC
## 10     mutation_subtypes        OncoSign
head(colData(acc)$Histology)
## [1] "Usual Type" "Usual Type" "Usual Type" "Usual Type" "Usual Type"
## [6] "Usual Type"

Converting TCGA UUIDs to barcodes and back

TCGAutils provides a number of ID translation functions. These allow the user to translate from either file or case UUIDs to TCGA barcodes and back. These functions work by querying the Genomic Data Commons API via the GenomicDataCommons package (thanks to Sean Davis). These include:

  • UUIDtoBarcode()
  • barcodeToUUID()
  • UUIDtoUUID()
  • filenameToBarcode()

See the TCGAutils help pages for details.

Other TCGA data types

Helper functions to add TCGA exon files (legacy archive), copy number and GISTIC copy number calls to MultiAssayExperiment objects are also available in TCGAutils.

Plotting, correlation, and other analyses

These provide exercises and solutions using the miniACC dataset.

How many miniACC samples have data for each combination of assays?

Solution

The built-in upsetSamples creates an “upset” Venn diagram to answer this question:

upsetSamples(miniACC)

In this dataset only 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.

Kaplan-meier plot stratified by pathology_N_stage

Create a Kaplan-meier plot, using pathology_N_stage as a stratifying variable.

Solution

The colData provides clinical data for things like a Kaplan-Meier plot for overall survival stratified by nodal stage.

Surv(miniACC$days_to_death, miniACC$vital_status)
##  [1] 1355  1677    NA+  423   365    NA+  490   579    NA+  922   551 
## [12] 1750    NA+ 2105    NA+  541    NA+   NA+  490    NA+   NA+  562 
## [23]   NA+   NA+   NA+   NA+   NA+   NA+  289    NA+   NA+   NA+  552 
## [34]   NA+   NA+   NA+  994    NA+   NA+  498    NA+   NA+  344    NA+
## [45]   NA+   NA+   NA+   NA+   NA+   NA+   NA+   NA+   NA+  391   125 
## [56]   NA+ 1852    NA+   NA+   NA+   NA+   NA+   NA+   NA+ 1204   159 
## [67] 1197   662   445    NA+ 2385   436  1105    NA+ 1613    NA+   NA+
## [78] 2405    NA+   NA+   NA+   NA+   NA+  207     0    NA+   NA+   NA+
## [89]   NA+   NA+   NA+  383

And remove any patients missing overall survival information:

miniACCsurv <- miniACC[,
    complete.cases(miniACC$days_to_death, miniACC$vital_status), ]
fit <- survfit(Surv(days_to_death, vital_status) ~ pathology_N_stage,
    data = colData(miniACCsurv))
ggsurvplot(fit, data = colData(miniACCsurv), 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"))
## harmonizing input:
##   removing 216 sampleMap rows not in names(experiments)
wideacc$y = Surv(wideacc$days_to_death, wideacc$vital_status)
head(wideacc)
## DataFrame with 6 rows and 7 columns
##        primary vital_status days_to_death pathology_N_stage
##    <character>    <integer>     <integer>       <character>
## 1 TCGA-OR-A5J1            1          1355                n0
## 2 TCGA-OR-A5J2            1          1677                n0
## 3 TCGA-OR-A5J3            0            NA                n0
## 4 TCGA-OR-A5J4            1           423                n1
## 5 TCGA-OR-A5J5            1           365                n0
## 6 TCGA-OR-A5J6            0            NA                n0
##   RNASeq2GeneNorm_EZH2 gistict_EZH2      y
##              <numeric>    <numeric> <Surv>
## 1              75.8886            0 1355:1
## 2             326.5332            1 1677:1
## 3              190.194            1   NA:0
## 4                   NA           -2  423:1
## 5             366.3826            1  365:1
## 6              30.7371            1   NA:0

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)
## Call:
## coxph(formula = Surv(days_to_death, vital_status) ~ gistict_EZH2 + 
##     log2(RNASeq2GeneNorm_EZH2) + pathology_N_stage, data = wideacc)
## 
##                                coef exp(coef) se(coef)      z        p
## gistict_EZH2               -0.03723   0.96345  0.28205 -0.132 0.894986
## log2(RNASeq2GeneNorm_EZH2)  0.97731   2.65729  0.28105  3.477 0.000506
## pathology_N_stagen1         0.37749   1.45862  0.56992  0.662 0.507743
## 
## Likelihood ratio test=16.28  on 3 df, p=0.0009942
## n= 26, number of events= 26 
##    (66 observations deleted due to missingness)

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.

back to top

Correlation between RNA-seq and copy number

Part 1

For all genes where there is both recurrent copy number (gistict assay) and RNA-seq, calculate the correlation between log2(RNAseq + 1) and copy number. Create a histogram of these correlations. Compare this with the histogram of correlations between all unmatched gene - copy number pairs.

Solution

First, narrow down miniACC to only the assays needed:

subacc <- miniACC[, , c("RNASeq2GeneNorm", "gistict")]
## harmonizing input:
##   removing 216 sampleMap rows not in names(experiments)

Align the rows and columns, keeping only samples with both assays available:

subacc <- intersectColumns(subacc)
subacc <- intersectRows(subacc)

Create a list of numeric matrices:

subacc.list <- assays(subacc)

Log-transform the RNA-seq assay:

subacc.list[[1]] <- log2(subacc.list[[1]] + 1)

Transpose both, so genes are in columns:

subacc.list <- lapply(subacc.list, t)

Calculate the correlation between columns in the first matrix and columns in the second matrix:

corres <- cor(subacc.list[[1]], subacc.list[[2]])

And finally, create the histograms:

hist(diag(corres))

hist(corres[upper.tri(corres)])

Part 2

For the gene with highest correlation to copy number, make a box plot of log2 expression against copy number.

Solution

First, identify the gene with highest correlation between expression and copy number:

which.max(diag(corres))
## EIF4E 
##    91

You could now make the plot by taking the EIF4E columns from each element of the list subacc.list list that was extracted from the subacc MultiAssayExperiment, but let’s do it by subsetting and extracting from the MultiAssayExperiment:

df <- wideFormat(subacc["EIF4E", , ])
head(df)
## DataFrame with 6 rows and 3 columns
##        primary RNASeq2GeneNorm_EIF4E gistict_EIF4E
##    <character>             <numeric>     <numeric>
## 1 TCGA-OR-A5J1              460.6148             0
## 2 TCGA-OR-A5J2              371.2252             0
## 3 TCGA-OR-A5J3              516.0717             0
## 4 TCGA-OR-A5J5             1129.3571             1
## 5 TCGA-OR-A5J6              822.0782             0
## 6 TCGA-OR-A5J7              344.5648            -1
boxplot(RNASeq2GeneNorm_EIF4E ~ gistict_EIF4E,
        data=df, varwidth=TRUE,
        xlab="GISTIC Relative Copy Number Call",
        ylab="RNA-seq counts")

back to top

Identifying correlated principal components

Perform Principal Components Analysis of each of the five assays, using samples available on each assay, log-transforming RNA-seq data first. Using the first 10 components, calculate Pearson correlation between all scores and plot these correlations as a heatmap to identify correlated components across assays.

Solution

Here’s a function to simplify doing the PCAs:

getLoadings <- function(x, ncomp=10, dolog=FALSE, center=TRUE, scale.=TRUE){
  if(dolog){
    x <- log2(x + 1)
  }
  pc = prcomp(x, center=center, scale.=scale.)
  return(t(pc$rotation[, 1:10]))
}

Although it would be possible to do the following with a loop, the different data types do require different options for PCA (e.g. mutations are a 0/1 matrix with 1 meaning there is a somatic mutation, and gistict varies between -2 for homozygous loss and 2 for a genome doubling, so neither make sense to scale and center). So it is just as well to do the following one by one, concatenating each PCA results to the MultiAssayExperiment:

miniACC2 <- intersectColumns(miniACC)
miniACC2 <- c(miniACC2, rnaseqPCA=getLoadings(assays(miniACC2)[[1]], dolog=TRUE), mapFrom=1L)
## Warning: Assuming column order in the data provided 
##  matches the order in 'mapFrom' experiment(s) colnames
miniACC2 <- c(miniACC2, gistictPCA=getLoadings(assays(miniACC2)[[2]], center=FALSE, scale.=FALSE), mapFrom=2L)
## Warning: Assuming column order in the data provided 
##  matches the order in 'mapFrom' experiment(s) colnames
miniACC2 <- c(miniACC2, proteinPCA=getLoadings(assays(miniACC2)[[3]]), mapFrom=3L)
## Warning: Assuming column order in the data provided 
##  matches the order in 'mapFrom' experiment(s) colnames
miniACC2 <- c(miniACC2, mutationsPCA=getLoadings(assays(miniACC2)[[4]], center=FALSE, scale.=FALSE), mapFrom=4L)
## Warning: Assuming column order in the data provided 
##  matches the order in 'mapFrom' experiment(s) colnames
miniACC2 <- c(miniACC2, miRNAPCA=getLoadings(assays(miniACC2)[[5]]), mapFrom=5L)
## Warning: Assuming column order in the data provided 
##  matches the order in 'mapFrom' experiment(s) colnames

Now subset to keep only the PCA results:

miniACC2 <- miniACC2[, , 6:10]
## harmonizing input:
##   removing 215 sampleMap rows not in names(experiments)
experiments(miniACC2)
## ExperimentList class object of length 5: 
##  [1] rnaseqPCA: matrix with 10 rows and 43 columns 
##  [2] gistictPCA: matrix with 10 rows and 43 columns 
##  [3] proteinPCA: matrix with 10 rows and 43 columns 
##  [4] mutationsPCA: matrix with 10 rows and 43 columns 
##  [5] miRNAPCA: matrix with 10 rows and 43 columns

Note, it would be equally easy (and maybe better) to do PCA on all samples available for each assay, then do intersectColumns at this point instead.

Now, steps for calculating the correlations and plotting a heatmap: * Use wideFormat to paste these together, which has the nice property of adding assay names to the column names. * The first column always contains the sample identifier, so remove it. * Coerce to a matrix * Calculate the correlations, and take the absolute value (since signs of principal components are arbitrary) * Set the diagonals to NA (each variable has a correlation of 1 to itself).

df <- wideFormat(miniACC2)[, -1]
mycors <- cor(as.matrix(df))
mycors <- abs(mycors)
diag(mycors) <- NA

To simplify the heatmap, show only components that have at least one correlation greater than 0.5.

has.high.cor <- apply(mycors, 2, max, na.rm=TRUE) > 0.5
mycors <- mycors[has.high.cor, has.high.cor]
pheatmap::pheatmap(mycors)

The highest correlation present is between PC2 of the RNA-seq assay, and PC1 of the protein assay.

back to top

Annotating with ranges

This section doesn’t use the TCGAutils shortcuts, and is more generally applicable.

Convert all the ExperimentList elements in miniACC to RangedSummarizedExperiment objects. Then use rowRanges to annotate these objects with genomic ranges. For the microRNA assay, annotate instead with the genomic coordinates of predicted targets.

Solution

The following shortcut function takes a list of human gene symbols and uses AnnotationFilter and EnsDb.Hsapiens.v86 to look up the ranges, and return these as a GRangesList which can be used to replace the rowRanges of the SummarizedExperiment objects:

getrr <- function(identifiers, EnsDbFilterFunc=AnnotationFilter::SymbolFilter) {
    edb <- EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86
    afl <- AnnotationFilter::AnnotationFilterList(
        EnsDbFilterFunc(identifiers),
        SeqNameFilter(c(1:21, "X", "Y")),
        TxBiotypeFilter("protein_coding"))
    gr <- genes(edb, filter=afl)
    grl <- split(gr, factor(identifiers))
    grl <- grl[match(identifiers, names(grl))]
    stopifnot(identical(names(grl), identifiers))
    return(grl)
}

For example:

getrr(rownames(miniACC)[[1]])
## GRangesList object of length 198:
## $DIRAS3
## GRanges object with 1 range and 7 metadata columns:
##                   seqnames          ranges strand |         gene_id
##                      <Rle>       <IRanges>  <Rle> |     <character>
##   ENSG00000116288        1 7954291-7985505      + | ENSG00000116288
##                     gene_name   gene_biotype seq_coord_system      symbol
##                   <character>    <character>      <character> <character>
##   ENSG00000116288       PARK7 protein_coding       chromosome       PARK7
##                   entrezid     tx_biotype
##                     <list>    <character>
##   ENSG00000116288    11315 protein_coding
##   -------
##   seqinfo: 22 sequences from GRCh38 genome
## 
## $MAPK14
## GRanges object with 1 range and 7 metadata columns:
##                   seqnames          ranges strand |         gene_id
##                      <Rle>       <IRanges>  <Rle> |     <character>
##   ENSG00000116285        1 8004404-8026308      - | ENSG00000116285
##                     gene_name   gene_biotype seq_coord_system      symbol
##                   <character>    <character>      <character> <character>
##   ENSG00000116285      ERRFI1 protein_coding       chromosome      ERRFI1
##                   entrezid     tx_biotype
##                     <list>    <character>
##   ENSG00000116285    54206 protein_coding
##   -------
##   seqinfo: 22 sequences from GRCh38 genome
## 
## $YAP1
## GRanges object with 1 range and 7 metadata columns:
##                   seqnames            ranges strand |         gene_id
##                      <Rle>         <IRanges>  <Rle> |     <character>
##   ENSG00000198793        1 11106535-11262507      - | ENSG00000198793
##                     gene_name   gene_biotype seq_coord_system      symbol
##                   <character>    <character>      <character> <character>
##   ENSG00000198793        MTOR protein_coding       chromosome        MTOR
##                   entrezid     tx_biotype
##                     <list>    <character>
##   ENSG00000198793     2475 protein_coding
##   -------
##   seqinfo: 22 sequences from GRCh38 genome
## 
## ...
## <195 more elements>

Use this to set the rowRanges of experiments 1-4 with these GRangesList objects

rseACC <- miniACC
withRSE <- c(1:3, 5)
for (i in withRSE){
  rowRanges(rseACC[[i]]) <- getrr(rownames(rseACC[[i]]))
}

Note that the class of experiments 1-4 is now RangedSummarizedExperiment:

experiments(rseACC)
## ExperimentList class object of length 5: 
##  [1] RNASeq2GeneNorm: RangedSummarizedExperiment with 198 rows and 79 columns 
##  [2] gistict: RangedSummarizedExperiment with 198 rows and 90 columns 
##  [3] RPPAArray: RangedSummarizedExperiment with 33 rows and 46 columns 
##  [4] Mutations: matrix with 97 rows and 90 columns 
##  [5] miRNASeqGene: RangedSummarizedExperiment with 471 rows and 80 columns

With ranged objects in the MultiAssayExperiment, you can then do subsetting by ranges. For example, select all genes on chromosome 1 for the four rangedSummarizedExperiment objects:

rseACC[GRanges(seqnames="1:1-1e9"), , withRSE]
## harmonizing input:
##   removing 90 sampleMap rows not in names(experiments)
## harmonizing input:
##   removing 80 sampleMap rows not in names(experiments)
## A MultiAssayExperiment object of 3 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 3: 
##  [1] RNASeq2GeneNorm: RangedSummarizedExperiment with 22 rows and 79 columns 
##  [2] gistict: RangedSummarizedExperiment with 22 rows and 90 columns 
##  [3] RPPAArray: RangedSummarizedExperiment with 3 rows and 46 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

Note about microRNA: You can set ranges for the microRNA assay according to the genomic location of those microRNA, or the locations of their predicted targets, but we don’t do it here. Assigning genomic ranges of microRNA targets would be an easy way to subset them according to their targets.

back to top


  1. City University of New York, New York, NY, USA

  2. Roswell Park Comprehensive Cancer Center, Buffalo, NY

  3. City University of New York, New York, NY, USA

  4. City University of New York, New York, NY, USA