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:

• 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)

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

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:

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

## harmonizing input:
##   removing 311 sampleMap rows with 'colname' not in colnames of experiments
##   removing 74 colData rownames not in sampleMap 'primary'

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

miniACC[, , "RNASeq2GeneNorm"]

## harmonizing input:
##   removing 13 colData rownames not in sampleMap 'primary'

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"]]

## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 198 features, 79 samples 
##   element names: exprs 
## protocolData: none
## phenoData: none
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:

Patients with complete data

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)

## harmonizing input:
##   removing 170 sampleMap rows with 'colname' not in colnames of experiments
##   removing 49 colData rownames not in sampleMap 'primary'

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")])
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

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

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"))

## DataFrame with 518 rows and 7 columns
##               assay      primary     rowname                      colname
##               <Rle>        <Rle> <character>                        <Rle>
## 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"))

## DataFrame with 92 rows and 9 columns
##          primary vital_status days_to_death Mutations_CTNNB1
##         <factor>    <integer>     <integer>        <numeric>
## 1   TCGA-OR-A5J1            1          1355                0
## 2   TCGA-OR-A5J2            1          1677                1
## 3   TCGA-OR-A5J3            0            NA                0
## 4   TCGA-OR-A5J4            1           423                0
## 5   TCGA-OR-A5J5            1           365                0
## ...          ...          ...           ...              ...
## 88  TCGA-PK-A5H8            0            NA                0
## 89  TCGA-PK-A5H9            0            NA                0
## 90  TCGA-PK-A5HA            0            NA                0
## 91  TCGA-PK-A5HB            0            NA                0
## 92  TCGA-PK-A5HC            0            NA                0
##     Mutations_TP53 RNASeq2GeneNorm_CTNNB1 RNASeq2GeneNorm_TP53
##          <numeric>              <numeric>            <numeric>
## 1                0               5634.467             563.4006
## 2                1              62658.391             165.4811
## 3                0               6337.426             956.3028
## 4                0                     NA                   NA
## 5                0               5979.055            1169.6359
## ...            ...                    ...                  ...
## 88               0               3033.648             737.6640
## 89               0               5258.986             890.8663
## 90               0               8120.165             683.5722
## 91               0               5257.815             237.3697
## 92               0                     NA                   NA
##     gistict_CTNNB1 gistict_TP53
##          <numeric>    <numeric>
## 1                0            0
## 2                1            0
## 3                0            0
## 4                0            1
## 5                0            0
## ...            ...          ...
## 88              NA           NA
## 89               0            0
## 90               0           -1
## 91              -1           -1
## 92               1            1

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: ExpressionSet with 198 rows and 79 columns 
##  [2] gistict: SummarizedExperiment with 198 rows and 90 columns 
##  [3] RPPAArray: ExpressionSet with 33 rows and 46 columns 
##  [4] Mutations: matrix with 97 rows and 90 columns 
##  [5] miRNASeqGene: ExpressionSet 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 in .local(x, ...): 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: ExpressionSet with 198 rows and 79 columns 
##  [2] gistict: SummarizedExperiment with 198 rows and 90 columns 
##  [3] RPPAArray: ExpressionSet with 33 rows and 46 columns 
##  [4] Mutations: matrix with 97 rows and 90 columns 
##  [5] miRNASeqGene: ExpressionSet with 471 rows and 80 columns 
##  [6] log2rnaseq: matrix with 198 rows and 79 columns

back to top

How many samples have data for each combination of assays?

Solution

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

upsetSamples(miniACC)

## Loading required namespace: UpSetR

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.

suppressPackageStartupMessages({
  library(survival)
  library(survminer)
})
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), ]

## harmonizing input:
##   removing 248 sampleMap rows with 'colname' not in colnames of experiments
##   removing 58 colData rownames not in sampleMap 'primary'

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"))
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
##       <factor>    <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.1940            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.0372    0.9635   0.2821 -0.13 0.89499
## log2(RNASeq2GeneNorm_EZH2)  0.9773    2.6573   0.2811  3.48 0.00051
## pathology_N_stagen1         0.3775    1.4586   0.5699  0.66 0.50774
## 
## Likelihood ratio test=16.3  on 3 df, p=0.000994
## 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")]

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

subacc <- intersectColumns(subacc)

## harmonizing input:
##   removing 15 sampleMap rows with 'colname' not in colnames of experiments
##   removing 15 colData rownames not in sampleMap 'primary'

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
##       <factor>             <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)

## harmonizing input:
##   removing 170 sampleMap rows with 'colname' not in colnames of experiments
##   removing 49 colData rownames not in sampleMap 'primary'

miniACC2 <- c(miniACC2, rnaseqPCA=getLoadings(assays(miniACC2)[[1]], dolog=TRUE), mapFrom=1L)

## Warning in .local(x, ...): 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 in .local(x, ...): 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 in .local(x, ...): 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 in .local(x, ...): 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 in .local(x, ...): 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]
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

Convert all the ExperimentList elements in miniACC to SummarizedExperiment 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

First, make a new object and convert its experiments to SummarizedExperiment objects:

suppressPackageStartupMessages(library(SummarizedExperiment))
seACC <- miniACC
experiments(seACC)

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

seACC[[1]] <- SummarizedExperiment(exprs(seACC[[1]]))
seACC[[3]] <- SummarizedExperiment(exprs(seACC[[3]]))
seACC[[4]] <- SummarizedExperiment(seACC[[4]])
seACC[[5]] <- SummarizedExperiment(exprs(seACC[[5]]))
experiments(seACC)

## 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: SummarizedExperiment with 97 rows and 90 columns 
##  [5] miRNASeqGene: SummarizedExperiment with 471 rows and 80 columns

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="SymbolFilter"){
  suppressPackageStartupMessages({
    library(AnnotationFilter)
    library(EnsDb.Hsapiens.v86)
  })
    FUN <- get(EnsDbFilterFunc)
    edb <- EnsDb.Hsapiens.v86
    afl <- AnnotationFilterList(FUN(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(seACC)[[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
## 
## $MAPK14 
## GRanges object with 1 range and 7 metadata columns:
##                   seqnames             ranges strand |         gene_id
##   ENSG00000116285        1 [8004404, 8026308]      - | ENSG00000116285
##                   gene_name   gene_biotype seq_coord_system symbol
##   ENSG00000116285    ERRFI1 protein_coding       chromosome ERRFI1
##                   entrezid     tx_biotype
##   ENSG00000116285    54206 protein_coding
## 
## $YAP1 
## GRanges object with 1 range and 7 metadata columns:
##                   seqnames               ranges strand |         gene_id
##   ENSG00000198793        1 [11106535, 11262507]      - | ENSG00000198793
##                   gene_name   gene_biotype seq_coord_system symbol
##   ENSG00000198793      MTOR protein_coding       chromosome   MTOR
##                   entrezid     tx_biotype
##   ENSG00000198793     2475 protein_coding
## 
## ...
## <195 more elements>
## -------
## seqinfo: 22 sequences from GRCh38 genome

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

for (i in 1:4){
  rowRanges(seACC[[i]]) <- getrr(rownames(seACC[[i]]))
}

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

experiments(seACC)

## 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: RangedSummarizedExperiment with 97 rows and 90 columns 
##  [5] miRNASeqGene: SummarizedExperiment 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:

seACC[GRanges(seqnames="1:1-1e9"), , 1:4]

## A MultiAssayExperiment object of 4 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 4: 
##  [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 
##  [4] Mutations: RangedSummarizedExperiment with 11 rows and 90 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

Shiny app

The maeView.R function defined in this workshop opens a Shiny app for similar TCGA objects, to identify and visualize relationships between RNA-seq expression, GISTIC copy number peaks, and microRNA expression. For a specified gene, you can view a boxplot of expression vs. copy number, and use limma to identify microRNA correlated to expression of that gene.

MultiAssayExperimentWorkshop::maeView(miniACC)

back to top