Generalized Canonical Correlation with missing individuals - Examples

Dolors Pelegrí-Sisó and Juan R. Gonzalez

Institute for Global Health (ISGlobal), Barcelona, Spain Bioinformatics Research Group in Epidemiolgy (BRGE) https://brge.isglobal.org

Package

mgcca 0.99.7

Contents

1 Introduction

This vignette shows how to integrate several data tables with mgcca, Generalized Canonical Correlation Analysis with missing individuals (whole missing rows; van de Velden and Bijmolt, 2006). A shared set of latent components is estimated across tables even when the tables do not share exactly the same individuals: an individual absent from a table simply does not contribute to that table.

The single entry point, mgcca(), has two backends chosen by whether you pass a filename. With a filename the whole numerical core runs in C++ on HDF5 files through the BigDataStatMeth API, so it scales from small tables to full-omics tables (hundreds of thousands of features) without ever loading a table completely into memory and without forming any variable-by-variable matrix; the results are written to the HDF5 file and can be pulled back into a compact R object for plotting with mgcca_results() (or mgcca(..., collect = TRUE)). With filename = NULL the same call fits in memory in pure R instead, for tables that fit in RAM – no HDF5 file is written. The two backends return the same "mgcca" object; the in-memory fit is not written to disk, so it is not reloadable with mgcca_load() and the reliability analyses (mgcca_sensitivity(), mgcca_stability()), which read the source blocks from HDF5, require a file-backed fit.

1.1 Related packages

Several R packages integrate multiple omics tables, and mgcca is best read beside them rather than on its own. RGCCA implements the regularized generalized canonical correlation framework this package builds on, with a wide choice of schemes, designs and sparsity options. mixOmics provides a large supervised and unsupervised multi-block toolbox, including sparse variants of PLS and CCA together with an extensive set of graphical outputs. MOFA2 takes a probabilistic latent-factor view, supporting several likelihoods per view and reporting the variance each factor explains in each of them. omicade4 and made4 come from the multiple co-inertia tradition and are built on the ade4 machinery, while mogsa combines multiple factor analysis with gene-set projection. What mgcca adds is narrow and specific: individuals missing from a table are a first-class case, so tables that do not cover the same participants are analysed on the union of individuals without imputing or dropping rows; the numerical core runs out-of-core on HDF5 through BigDataStatMeth, so a table with hundreds of thousands of features never has to fit in memory; and the package ships a reliability layer – availability audit, sensitivity and stability – for asking how much of a result rests on the design and on the sample. That is a difference in scope, not a claim of superiority: no comparison of accuracy or run time is made here, and which package suits an analysis depends on whether missing individuals, table size, or the modelling options the others offer are what it turns on.

2 Installation

mgcca is installed from Bioconductor:

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

BiocManager::install("mgcca")

library(mgcca)

3 A self-contained example

The package ships a small multi-table cardiovascular data set with three tables that do not share the same individuals – exactly the situation mgcca is designed for. Each table has individuals in rows (identified by their row names) and variables in columns.

The data set is simulated. It is generated from aggregate summary statistics of a real, de-identified cardiovascular methylation study, so it reproduces that study’s correlation and missingness structure, but no row corresponds to a real participant and the results below support no substantive conclusion. See ?cardiovascular for the generating script.

data(cardiovascular, package = "mgcca")

X <- list(methylation = as.matrix(X1),
          clinical    = as.matrix(X2),
          other       = as.matrix(X3))
X <- lapply(X, function(m) { storage.mode(m) <- "double"; m })

# Different number of individuals per table -> missing individuals
sapply(X, dim)

##      methylation clinical other
## [1,]         618      633   645
## [2,]          69        6     4

length(Reduce(union, lapply(X, rownames)))   # size of the individual union

## [1] 646

3.1 Running mgcca

mgcca() imports the tables into an HDF5 file, scales each variable (column-centre and standardise, done out-of-core inside the file), and computes the shared components. The tables here are small (few variables relative to individuals), so method = "solve" is well conditioned; for wide/omics tables use "penalized" (with a lambda per table) or "geninv".

Passing collect = TRUE returns an object of class "mgcca" (holding Y, the per-table corsY/pval.cor, scores and AVE) ready for the plotting helpers. The full results also remain in the HDF5 file.

h5file <- tempfile(fileext = ".h5")

fit <- mgcca(X, filename = h5file, nfac = 2,
             method = "solve", scores = TRUE, collect = TRUE)

fit                       # print method: a compact overview

## Generalized Canonical Correlation Analysis (mgcca)
##   tables      : 3  (methylation, clinical, other)
##   individuals : 646
##   components  : 2
##   eigenvalues : 0.7271, 0.6393
##   AVE (outer) : 0.1692, 0.13
##   spectral gap: 0.03843 (relative, external at L = 2)
##   scores      : available (per table)

summary() adds the variance explained per table and the variables most correlated with each shared component:

summary(fit, top = 4)

## Generalized Canonical Correlation Analysis (mgcca)
##   tables      : 3  (methylation, clinical, other)
##   individuals : 646
##   components  : 2
##   eigenvalues : 0.7271, 0.6393
##   AVE (outer) : 0.1692, 0.13
##   spectral gap: 0.03843 (relative, external at L = 2)
##   scores      : available (per table)
## 
## AVE per table (outer model):
##              comp1  comp2
## methylation 0.1625 0.1343
## clinical    0.1177 0.0059
## other       0.3633 0.2415
## 
## Pairwise overlap (n = 646 individuals in the union):
##                      n_jk alpha_jk
## methylation-clinical  606   0.9381
## methylation-other     617   0.9551
## clinical-other        633   0.9799
##   n_j: methylation = 618, clinical = 633, other = 645
##   alpha_jk = n_jk / n;  sd(alpha) = 0.0210 (cv = 0.022),  range 0.9381 - 0.9799
## 
## Top 4 variables per component (correlation with Y):
##   [methylation]
##     comp1: cg3 (+0.72), cg34 (+0.72), cg24 (+0.70), cg1 (+0.69)
##     comp2: cg9 (+0.72), cg22 (+0.69), cg21 (+0.64), cg23 (+0.63)
##   [clinical]
##     comp1: TC (+0.52), LDL (+0.46), HDL (+0.42), Gluc (-0.21)
##     comp2: HDL (-0.13), BMI (+0.08), waist (+0.07), TC (-0.07)
##   [other]
##     comp1: Gran (-0.80), CD4T (+0.79), Mono (-0.38), NK (+0.21)
##     comp2: NK (+0.82), CD4T (-0.42), Gran (-0.33), Mono (+0.02)

3.2 Individuals plot

plotIndividuals() projects the individuals onto the two shared components and returns a ggplot object you can further style. A grouping factor – a phenotype, a clinical category, a cluster – can be supplied to colour the points.

To pick a meaningful one, look at which clinical variables actually drive the shared components (their correlation with each component):

round(fit$corsY[["clinical"]], 3)

##        comp1  comp2
## TC     0.516 -0.066
## LDL    0.455  0.012
## HDL    0.416 -0.128
## BMI    0.060  0.080
## waist -0.113  0.074
## Gluc  -0.208  0.053

Cholesterol (LDL, total cholesterol) loads most on the first component, whereas BMI barely does – so LDL categories will show structure while BMI would not. We build LDL clinical categories (NCEP ATP III cut-points) and align them to the individuals of the analysis (those absent from the clinical table stay NA):

ldl <- X$clinical[, "LDL"]                        # real LDL values (mg/dL)
ldl <- ldl[rownames(fit$Y)]                       # align to the individual union
ldl_group <- cut(ldl, breaks = c(-Inf, 100, 160, Inf),
                 labels = c("optimal", "borderline", "high"))
table(ldl_group, useNA = "ifany")

## ldl_group
##    optimal borderline       high       <NA> 
##         81        409        143         13

plotIndividuals(fit, group = ldl_group)

The 95% confidence ellipses (one per group) make the shift along the first component clear – the optimal and high ellipses sit at opposite ends of it, with borderline between them – even though the individual points overlap heavily.

Contrast this with BMI, which (as the loadings above anticipated) barely aligns with the shared components: colouring by BMI category gives essentially concentric ellipses, i.e. no separation.

bmi <- X$clinical[, "BMI"][rownames(fit$Y)]
bmi_group <- cut(bmi, breaks = c(-Inf, 25, 30, Inf),
                 labels = c("normal", "overweight", "obese"))
plotIndividuals(fit, group = bmi_group, title = "Individuals by BMI category")

The ellipses can be switched off with ellipse = FALSE.

3.3 Variables plot

plotVariables() draws, for one table, the correlation of each variable with the two shared components as a correlation circle. Only the top variables (by correlation radius) are labelled, so it stays readable for omics-scale tables.

plotVariables(fit, table = "clinical", top = 6)

plotLoadings() is a complementary view: the signed correlation of the top variables with a single component as a bar chart (here it makes the cholesterol-vs-glucose split of component 1 explicit).

plotLoadings(fit, table = "clinical", comp = 1, top = 6)

3.4 Variance explained

plotAVE() shows how much of each table’s variance each shared component captures (outer-model AVE); plotScree() shows the component eigenvalues.

plotAVE(fit)

plotScree(fit)

3.5 Significant variables

getSignif() lists the variables significantly correlated with any component and topVars() returns the strongest ones.

head(getSignif(fit, pval.cut = 1e-3))

##   variable       table
## 1      cg1 methylation
## 2      cg2 methylation
## 3      cg3 methylation
## 4      cg4 methylation
## 5      cg5 methylation
## 6      cg6 methylation

topVars(fit, axis = 1, topN = 5)

## $methylation
## [1] "cg3"  "cg34" "cg24" "cg1" 
## 
## $clinical
## [1] "TC"  "LDL" "HDL" "BMI"
## 
## $other
## [1] "CD4T" "NK"   "Mono" "Gran"

Both read the p-values mgcca() stored beside the correlations. Those are the two-sided p-values of a correlation coefficient, and corTestP() is the small helper that produces them from a correlation and the number of individuals behind it – which is what you need when the correlation is one you computed yourself rather than one the fit carries:

r <- fit$corsY[["clinical"]][, 1]
p <- corTestP(r, nrow(X$clinical))      # n of the table the correlation came from
all.equal(p, fit$pval.cor[["clinical"]][, 1])

## [1] TRUE

signif(head(p, 4), 3)

##       TC      LDL      HDL      BMI 
## 2.15e-44 9.87e-34 6.95e-28 1.31e-01

3.6 The two plotting families

plotIndividuals(), plotVariables(), plotLoadings(), plotAVE(), plotScree(), plotScores() and plotBiplot() all return ggplot objects, so they can be restyled, composed or saved with the usual ggplot2 verbs. plotInds() and plotVars() are the package’s older base-graphics counterparts of the first two: they draw straight to the current device and take base graphical parameters, which is convenient inside a script that is already drawing with base graphics. They answer the same two questions, so there is no need to demonstrate both families in full.

plotInds(fit, group = ldl_group)

plotVars(fit, df = 2, nlab = 5)          # df = 2: the clinical table

4 Saving and reloading results

Every mgcca() run writes its results to the HDF5 file, and together with them a small provenance manifest stored as native HDF5 attributes: the method, lambda, the algebra route taken per table, the number of components, the table names, the eigenvalues and the package version/date. The manifest is attached to the FINAL_RESULTS group (run-level facts) and to each FINAL_RESULTS/corsY/<table> dataset (the per-table lambda and route_dual), so the results file is self-describing.

Because of this, a finished analysis can be reopened later from the file alone – no need to keep the R session, a saved .rds, or even to remember how it was run. mgcca_load() reads the manifest and the numeric blocks and rebuilds the full "mgcca" object:

# forget the in-memory object; everything needed is on disk
rm(fit)

fit <- mgcca_load(h5file)          # rebuilt straight from the HDF5 file
fit

## Generalized Canonical Correlation Analysis (mgcca)
##   tables      : 3  (methylation, clinical, other)
##   individuals : 646
##   components  : 2
##   eigenvalues : 0.7271, 0.6393
##   AVE (outer) : 0.1692, 0.13
##   spectral gap: 0.03843 (relative, external at L = 2)
##   scores      : available (per table)

The reconstructed provenance is available on the object’s descriptor:

desc <- attr(fit, "desc")
desc$method                        # "solve"

## [1] "solve"

desc$route                         # "auto"

## [1] "auto"

desc$route_dual                    # per-table: covariance vs Gram route

## methylation    clinical       other 
##       FALSE       FALSE       FALSE

desc$eig_values                    # shared-component eigenvalues

## [1] 0.7271469 0.6392769

All the plotting and post-analysis helpers work on the reloaded object exactly as before:

plotIndividuals(fit, group = ldl_group)

If you only need the numeric blocks (no provenance, no "mgcca" wrapper), the lower-level reader mgcca_results() returns them directly. It does not read the manifest, so it has to be told which tables to collect – either by handing it the descriptor mgcca(..., collect = FALSE) returns, or by naming them:

res <- mgcca_results(h5file, datasets = names(X), outputs = "Y")
dim(res$Y)

## [1] 646   2

identical(res$Y, fit$Y)      # same numbers, read back from the file

## [1] TRUE

outputs is subtractive: the datasets behind the components you do not ask for are never read, and every component you do ask for is the one a full collection would have produced.

5 How it scales to omics data

The same file-backed call (with a filename) works on omics-scale tables. Two things make this possible:

• Out-of-core, in HDF5. Import writes the raw tables to the file; scaling, cross-products, the eigendecomposition and the scores are all computed by BigDataStatMeth reading the data in blocks. On this path nothing forces a whole table into RAM. (The in-memory backend, filename = NULL, is the opposite trade-off: it holds the tables in RAM and is meant for tables that fit, not for omics-scale data.)
• Per-table algebra route (route = "auto"). For a table with more variables than individuals (p > n, e.g. methylation with hundreds of thousands of CpGs) mgcca uses the Gram route (an n x n matrix), never the p x p covariance matrix, which would be intractable. For p < n it uses the covariance route. route = "cov" or "dual" force one route for every table.

So the very same mgcca(..., method = "penalized", lambda = ...) that runs on a 100-variable subset also runs, unchanged, on the full methylation table (~485k features), completing in seconds without exhausting memory.

6 A live multi-omics example: miniACC

miniACC is a small MultiAssayExperiment that ships inside the MultiAssayExperiment package itself. It is a subset of the TCGA adrenocortical carcinoma (ACC) cohort prepared for examples and documentation: 92 patients, five assays cut down to a couple of hundred features each, and the cohort’s clinical table. It is not a data set to draw conclusions from – the features were selected to make a compact example, not to represent the assays – but it keeps the two properties this section needs, and it needs no download: real assays measured on overlapping but different sets of patients, and a real clinical annotation to read the components against.

Everything from here to the end of this section runs when the vignette is built. The full-scale version of the same analysis, on the complete ACC cohort, is in the appendix and is shown rather than run.

library(MultiAssayExperiment)
library(SummarizedExperiment)

data("miniACC", package = "MultiAssayExperiment")
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
## Functionality:
##  experiments() - obtain the ExperimentList instance
##  colData() - the primary/phenotype DataFrame
##  sampleMap() - the sample coordination DataFrame
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment
##  *Format() - convert into a long or wide DataFrame
##  assays() - convert ExperimentList to a SimpleList of matrices
##  exportClass() - save data to flat files

6.1 Assembling the tables

We keep two of the five assays – gene expression and the RPPA protein array. They differ in how many patients they cover, which is the case mgcca is for.

acc_mae <- miniACC[, , c("RNASeq2GeneNorm", "RPPAArray")]

mgccaImpute() fills the missing entries within each assay by k-nearest neighbours. That is a different problem from the missing individuals mgcca() handles natively: here it is scattered holes inside a table, there it is whole rows absent from a table. With impute.zero = TRUE the exact zeros are treated as missing first, which is the usual convention for these count-derived tables (a handful of near-empty rows fall back to the column mean, and impute::impute.knn() says so).

acc_mae <- mgccaImpute(acc_mae, method = "knn", k = 10, impute.zero = TRUE)

getTables() turns the assays into the list of tables mgcca() expects: assays are stored features-by-samples, so each one is transposed on the way out. Their row names are full TCGA barcodes, and the patient identifier – what the two tables have to agree on – is their first 12 characters.

tabs <- getTables(acc_mae)
class(tabs)

## [1] "ListMAE"

vapply(tabs, dim, integer(2))                      # individuals x variables

##      RNASeq2GeneNorm RPPAArray
## [1,]              79        46
## [2,]             198        33

for (a in seq_along(tabs))
  rownames(tabs[[a]]) <- substr(rownames(tabs[[a]]), 1, 12)

length(Reduce(union, lapply(tabs, rownames)))      # individuals in the analysis

## [1] 79

length(Reduce(intersect, lapply(tabs, rownames)))  # measured on both layers

## [1] 46

So the analysis runs on 79 patients, of whom 46 carry both layers and 33 carry expression only. Those 33 are not dropped and nothing is imputed for them: they contribute to the expression table and to nothing else.

Some assays arrive from curatedTCGAData as character matrices – the ACC methylation table is the usual example – and mgcca() needs numeric tables. miniACC carries no such assay, so we put the expression table into that shape to exercise the conversion step a downloaded table needs:

storage.mode(tabs[[1]]) <- "character"
tabs[[1]][1, 1]

## [1] "1487.0317"

tabs <- matrixChr2Num(tabs, index = 1)
tabs[[1]][1, 1]

## [1] 1487.032

6.2 Choosing the penalty

Both tables have more variables than individuals, so the fit is penalized and needs a lambda per table. Rather than picking one by hand, mgcca_select_lambda() chooses it by held-out cross-block agreement: it splits the individuals into K folds, fits on each training part and measures how well the tables agree on the held-out individuals.

sel <- mgcca_select_lambda(tabs, nfac = 2, K = 5, seed = 1)

## 9 of 9 candidates are within one standard error of the best; the penalty is not sharply identified. Report the region.

sel

## mgcca ridge selection by held-out cross-block agreement
##   criterion   : cv_crossblock_agreement (5-fold over individuals, rule = argmax)
##   components  : 2
##   gamma       : 1   (lambda_j = gamma * s_j)
##   lambda      : RNASeq2GeneNorm = 198, RPPAArray = 45
##   agreement   : 0.4995   (mean held-out cross-block RV)
##   near-optimal: 9 of 9 candidates within 1 SE
##   Instability is a diagnostic at this choice, never the way to make it.

Read the message as well as the number. On a sample this small the criterion is flat: every candidate lies within one standard error of the best, so the selected value is a reasonable point inside a wide region rather than a sharply identified optimum. That is a statement worth reporting alongside the fit, and it is the reason the function prints it instead of returning a bare number. seed = 1 fixes the fold split; without it the folds come from the session’s own random stream and the selection is not reproducible across sessions.

6.3 Integrating the two layers

h5acc <- tempfile(fileext = ".h5")

acc <- mgcca(tabs, filename = h5acc, nfac = 2, method = "penalized",
             lambda = sel$lambda, scores = TRUE, route = "auto",
             collect = TRUE)
acc

## Generalized Canonical Correlation Analysis (mgcca)
##   tables      : 2  (RNASeq2GeneNorm, RPPAArray)
##   individuals : 79
##   components  : 2
##   eigenvalues : 0.8468, 0.8117
##   AVE (outer) : 0.1141, 0.08617
##   spectral gap: 0.0658 (relative, external at L = 2)
##   scores      : available (per table)

summary(acc, top = 4)

## Generalized Canonical Correlation Analysis (mgcca)
##   tables      : 2  (RNASeq2GeneNorm, RPPAArray)
##   individuals : 79
##   components  : 2
##   eigenvalues : 0.8468, 0.8117
##   AVE (outer) : 0.1141, 0.08617
##   spectral gap: 0.0658 (relative, external at L = 2)
##   scores      : available (per table)
## 
## AVE per table (outer model):
##                  comp1  comp2
## RNASeq2GeneNorm 0.1081 0.0822
## RPPAArray       0.1498 0.1102
## 
## Pairwise overlap (n = 79 individuals in the union):
##                           n_jk alpha_jk
## RNASeq2GeneNorm-RPPAArray   46   0.5823
##   n_j: RNASeq2GeneNorm = 79, RPPAArray = 46
##   alpha_jk = n_jk / n;  sd(alpha) = 0.0000 (cv = 0.000),  range 0.5823 - 0.5823
## 
## Top 4 variables per component (correlation with Y):
##   [RNASeq2GeneNorm]
##     comp1: STMN1 (-0.78), RAD51 (-0.78), CDK1 (-0.78), TYMS (-0.76)
##     comp2: CAV1 (+0.82), NOTCH1 (+0.76), LCK (+0.75), AXL (+0.72)
##   [RPPAArray]
##     comp1: MSH6 (-0.70), PRDX1 (+0.69), TSC1 (-0.67), XBP1 (+0.66)
##     comp2: BRCA2 (+0.69), ERCC1 (-0.67), ACVRL1 (+0.66), FASN (-0.56)

The overlap report in the summary is the design in one line: 46 of the 79 individuals are shared by the two tables, and alpha_jk is that count over the union.

The results are also on disk, and mgcca_results() reads them back without the "mgcca" wrapper. It does not read the provenance manifest, so it has to be told which tables to collect:

res <- mgcca_results(h5acc, datasets = names(tabs), outputs = c("Y", "AVE"))
dim(res$Y)

## [1] 79  2

names(Filter(Negate(is.null), res))

## [1] "Y"   "AVE"

6.4 Reading the components

plotScores() places the individuals of one table on its own scores, coloured by whatever grouping you supply; plotBiplot() overlays the individuals and the most correlated variables of a table on the same axes.

ids    <- rownames(acc$Y)
cd     <- as.data.frame(colData(miniACC))
stage  <- factor(as.character(cd[ids, "pathologic_stage"]))
vital  <- factor(ifelse(cd[ids, "vital_status"] == 1, "deceased", "alive"))
gender <- factor(as.character(cd[ids, "gender"]))
table(stage, useNA = "ifany")

## stage
##   stage i  stage ii stage iii  stage iv      <NA> 
##         9        37        16        15         2

plotScores(acc, table = 1, group = vital)

plotBiplot(acc, table = 2, top = 6)

mgcca_associate() makes the same question quantitative: it regresses each shared component on a phenotype and reports the variance explained, its p-value and a multiplicity-adjusted version.

assoc <- mgcca_associate(acc, list(stage = stage, vital = vital,
                                   gender = gender))
assoc

##   phenotype component           R2            p  n        p.adj
## 1     stage         1 0.2295985021 2.523281e-04 77 7.569843e-04
## 2     stage         2 0.0424506190 3.635259e-01 77 5.452889e-01
## 3     vital         1 0.4061435170 2.703819e-10 79 1.622291e-09
## 4     vital         2 0.0107895687 3.622986e-01 79 5.452889e-01
## 5    gender         1 0.0007690258 8.083159e-01 79 8.083159e-01
## 6    gender         2 0.0022647380 6.770612e-01 79 8.083159e-01

The first component tracks vital status and tumour stage; gender, which is in the table precisely as a control, does not. That is the shape one hopes to see – the shared signal lining up with disease rather than with a demographic variable – and on 79 patients from a data set built for examples it is an illustration of the reading, not a finding.

6.5 Are the components themselves significant?

mgcca_permtest() breaks the correspondence between individuals across tables and refits, giving a null distribution for the eigenvalues. It refits once per permutation, so the number of permutations sets the cost – and also the resolution: the smallest attainable p-value is 1 / (nperm + 1).

set.seed(1)
pt <- mgcca_permtest(tabs, nperm = 9, nfac = 2, method = "penalized",
                     lambda = sel$lambda)
pt$p.value

## [1] 0.1 0.1

range(pt$null)

## [1] 0.7263080 0.8385205

Nine permutations keep the vignette quick and are too few to test anything: 0.1 is the smallest p-value this call can return, so what is shown here is the mechanism, not evidence. A real use takes nperm = 99 or more, which is why the argument defaults to 99.

6.6 Scoring individuals the fit never saw

predict() applies the fitted weights – with the centring and scaling of the training tables – to new rows of a table, so a second cohort can be placed on the same axes. Here the “new” rows are three individuals of the protein table, which simply reproduces their scores:

new_rows <- tabs$RPPAArray[1:3, , drop = FALSE]
predict(acc, newdata = list(RPPAArray = new_rows))

## $RPPAArray
##                   comp1      comp2
## TCGA-OR-A5J2  0.1359425 -0.1493703
## TCGA-OR-A5J3 -0.2464442 -0.7029714
## TCGA-OR-A5J6  1.0704808  0.8479712

The table names of newdata must match the fit’s, and each one must carry every variable used at fit time; tables the fit does not know are ignored.

The HDF5 file is a temporary one here. A real analysis keeps it: it is what mgcca_load() reopens and what the reliability functions read.

BigDataStatMeth::hdf5_close_all()
unlink(h5acc)

7 Auditing the availability design before fitting

Everything above starts from a fit. mgcca_audit() runs before one, and answers a different kind of question: what does the availability design itself support? It reads who is present in which block – presence and absence, never a data value – and never calls mgcca(), so nothing a fit computes can move because of it.

Take a four-block design over 85 participants. Sixty of them carry methylation, expression and proteins; the remaining twenty-five carry proteins and metabolites. Nobody carries methylation and metabolites together.

ids   <- sprintf("ind%02d", 1:85)
core  <- 1:60                          # methylation + expression + proteins
extra <- 61:85                         # proteins + metabolites

avail <- cbind(
  methylation = seq_along(ids) %in% core,
  expression  = seq_along(ids) %in% core,
  proteins    = rep(TRUE, length(ids)),
  metabolites = seq_along(ids) %in% extra)
rownames(avail) <- ids

cohort  <- rep(c("cohort A", "cohort B"), length.out = length(ids))
avail_w <- ifelse(seq_along(ids) %in% core, 1, 3)   # declared design weights

colSums(avail)

## methylation  expression    proteins metabolites 
##          60          60          85          25

L is the rank of the working model you intend to fit; the audit records it and compares it against the support counts descriptively. strata and weights are what you declare about the design, and are given in the audit’s own row order.

audit <- mgcca_audit(avail, L = 2, strata = cohort, weights = avail_w)
audit

## mgcca availability audit (pre-fit; the estimator is not called)
##   source      : availability table
##   individuals : 85 ; blocks : 4 ; declared rank L = 2
##   strata      : 2 declared (cohort A, cohort B)
##   weights     : declared (ESS reported in addition to the raw counts)
##   design lines: ess_min = 30 ; availability_min = 0.05 ; max_weight = 20
## 
## Support graph
##   components 1 (connected) ; edges 4 of 6 pairs
##   design-support fragility: 1 bridge, 1 articulation block
## 
## Pair codes
##   DIRECT_SUPPORT           4
##   MODEL_TRANSFER_ELIGIBLE  2
##   UNSUPPORTED              0
##   STRUCTURAL_ZERO          0
## 
## Pairs carrying an advisory
##       block1      block2                    code                  advisory
##  methylation metabolites MODEL_TRANSFER_ELIGIBLE LOW_SUPPORT;SINGLE_BRIDGE
##   expression metabolites MODEL_TRANSFER_ELIGIBLE LOW_SUPPORT;SINGLE_BRIDGE
##     proteins metabolites          DIRECT_SUPPORT LOW_SUPPORT;SINGLE_BRIDGE
## 
## Stratum-block cells carrying an advisory
##   stratum       block n_obs proportion ess max_weight
##  cohort A metabolites    13     0.3023  13          3
##  cohort B metabolites    12     0.2857  12          3
##                       advisory
##  SEVERE_REPRESENTATION_WARNING
##  SEVERE_REPRESENTATION_WARNING
## 
##   structural_checks: not_evaluated. Empirical rank and conditioning
##   are not checked here, so a MODEL_TRANSFER_ELIGIBLE pair is identified
##   conditional on the declared rank-L model and its required structural
##   conditions, which are not empirically checked by this audit.
##   The design lines above are declared lines with margin. Every value
##   that raised a code is reported beside it; nothing here is a verdict,
##   and no design makes this audit refuse to return one.

Four pairs are carried directly by participants observed in both blocks. Methylation–metabolites and expression–metabolites are not: they have no shared participant at all, but a route of supported pairs runs between them, so they are coded MODEL_TRANSFER_ELIGIBLE. Such a pair is identified conditional on the declared rank-L model and its required structural conditions, which are not empirically checked by this audit – $structural_checks says not_evaluated, and at this level it always will. The route itself, and the support of its weakest edge, are on the record:

audit$provenance$pairs[, c("block1", "block2", "code", "n_jk", "path",
                           "bottleneck")]

##        block1      block2                    code n_jk
## 1 methylation  expression          DIRECT_SUPPORT   60
## 2 methylation    proteins          DIRECT_SUPPORT   60
## 3 methylation metabolites MODEL_TRANSFER_ELIGIBLE    0
## 4  expression    proteins          DIRECT_SUPPORT   60
## 5  expression metabolites MODEL_TRANSFER_ELIGIBLE    0
## 6    proteins metabolites          DIRECT_SUPPORT   25
##                                     path bottleneck
## 1              methylation -> expression         60
## 2                methylation -> proteins         60
## 3 methylation -> proteins -> metabolites         25
## 4                 expression -> proteins         60
## 5  expression -> proteins -> metabolites         25
## 6                proteins -> metabolites         25

The whole connection to metabolites rests on the single proteins–metabolites edge and the 25 participants behind it. That is design-support fragility: \(\kappa_{jk}\) counts the edge-disjoint routes between two blocks, so \(\kappa = 1\) means one edge failing would leave the two sides unconnected. It is a statement about the design’s topology, and says nothing about how well anything is estimated.

audit$graph$bridges

##     block1      block2 n_jk s_e s_e_ess kappa is_bridge
## 4 proteins metabolites   25  25      25     1      TRUE

audit$graph$articulation

## methylation  expression    proteins metabolites 
##       FALSE       FALSE        TRUE       FALSE

summary() adds the availability and representation tables. The three design lines (ess_min, availability_min, max_weight) are declared lines carrying margin, set by the analyst and configurable; they are not thresholds at which something fails and they are not estimated from the data. The value that raised a code is always printed beside it, so how close to the line a cell sits stays visible: the two metabolites cells below sit at an effective sample size of 13.0 and 12.0 against a declared line of 30, which is not the same situation as one sitting at 29.8.

summary(audit)

## mgcca availability audit (pre-fit; the estimator is not called)
##   source      : availability table
##   individuals : 85 ; blocks : 4 ; declared rank L = 2
##   strata      : 2 declared (cohort A, cohort B)
##   weights     : declared (ESS reported in addition to the raw counts)
##   design lines: ess_min = 30 ; availability_min = 0.05 ; max_weight = 20
## 
## Support graph
##   components 1 (connected) ; edges 4 of 6 pairs
##   design-support fragility: 1 bridge, 1 articulation block
## 
## Pair codes
##   DIRECT_SUPPORT           4
##   MODEL_TRANSFER_ELIGIBLE  2
##   UNSUPPORTED              0
##   STRUCTURAL_ZERO          0
## 
## Pairs carrying an advisory
##       block1      block2                    code                  advisory
##  methylation metabolites MODEL_TRANSFER_ELIGIBLE LOW_SUPPORT;SINGLE_BRIDGE
##   expression metabolites MODEL_TRANSFER_ELIGIBLE LOW_SUPPORT;SINGLE_BRIDGE
##     proteins metabolites          DIRECT_SUPPORT LOW_SUPPORT;SINGLE_BRIDGE
## 
## Stratum-block cells carrying an advisory
##   stratum       block n_obs proportion ess max_weight
##  cohort A metabolites    13     0.3023  13          3
##  cohort B metabolites    12     0.2857  12          3
##                       advisory
##  SEVERE_REPRESENTATION_WARNING
##  SEVERE_REPRESENTATION_WARNING
## 
##   structural_checks: not_evaluated. Empirical rank and conditioning
##   are not checked here, so a MODEL_TRANSFER_ELIGIBLE pair is identified
##   conditional on the declared rank-L model and its required structural
##   conditions, which are not empirically checked by this audit.
##   The design lines above are declared lines with margin. Every value
##   that raised a code is reported beside it; nothing here is a verdict,
##   and no design makes this audit refuse to return one.
## 
## Availability by block
##        block  n n_j pi_hat_j
##  methylation 85  60   0.7059
##   expression 85  60   0.7059
##     proteins 85  85   1.0000
##  metabolites 85  25   0.2941
## 
## Availability by declared stratum
##   stratum       block n_stratum n_jc pi_hat_jc
##  cohort A methylation        43   30    0.6977
##  cohort A  expression        43   30    0.6977
##  cohort A    proteins        43   43    1.0000
##  cohort A metabolites        43   13    0.3023
##  cohort B methylation        42   30    0.7143
##  cohort B  expression        42   30    0.7143
##  cohort B    proteins        42   42    1.0000
##  cohort B metabolites        42   12    0.2857
## 
## Pairs
##       block1      block2                    code n_jk pi_jkc_min s_e s_e_ess
##  methylation  expression          DIRECT_SUPPORT   60     0.6977  60      60
##  methylation    proteins          DIRECT_SUPPORT   60     0.6977  60      60
##  methylation metabolites MODEL_TRANSFER_ELIGIBLE    0     0.0000   0       0
##   expression    proteins          DIRECT_SUPPORT   60     0.6977  60      60
##   expression metabolites MODEL_TRANSFER_ELIGIBLE    0     0.0000   0       0
##     proteins metabolites          DIRECT_SUPPORT   25     0.2857  25      25
##  kappa                                   path bottleneck
##      2              methylation -> expression         60
##      2                methylation -> proteins         60
##      1 methylation -> proteins -> metabolites         25
##      2                 expression -> proteins         60
##      1  expression -> proteins -> metabolites         25
##      1                proteins -> metabolites         25
## 
##   s_e is the raw joint count n_jk. s_e_ess is the additional ESS-based
##   support; it does not replace s_e.
##   bottleneck is the widest-route value: the largest support a single
##   route can guarantee on its weakest edge. kappa counts edge-disjoint
##   routes -- design-support fragility, not identification robustness.
##   n_jk >= L is descriptive support; it does not prove numerical rank L.
## 
## Representation
##   stratum       block n_stratum n_obs proportion   ess max_weight
##  cohort A methylation        43    30     0.6977 30.00          1
##  cohort A  expression        43    30     0.6977 30.00          1
##  cohort A    proteins        43    43     1.0000 32.39          3
##  cohort A metabolites        43    13     0.3023 13.00          3
##  cohort B methylation        42    30     0.7143 30.00          1
##  cohort B  expression        42    30     0.7143 30.00          1
##  cohort B    proteins        42    42     1.0000 31.57          3
##  cohort B metabolites        42    12     0.2857 12.00          3
##  weight_concentration                      advisory
##                0.0000                              
##                0.0000                              
##                0.2468                              
##                0.0000 SEVERE_REPRESENTATION_WARNING
##                0.0000                              
##                0.0000                              
##                0.2484                              
##                0.0000 SEVERE_REPRESENTATION_WARNING
## 
##   A finite entry in this table evidences positivity for that cell.
##   It does not prove it: an empty cell may be empty by chance or
##   impossible by design, and the table cannot tell the two apart.
## 
##   Pairs coded MODEL_TRANSFER_ELIGIBLE are identified conditional on the declared rank-2 model and its required structural conditions, which are not empirically checked by this audit. n_jk >= L and bottleneck >= L are descriptive support only: they do not prove numerical rank 2.

Two cells carry SEVERE_REPRESENTATION_WARNING: metabolites is thin in both cohorts. Note also what the table can and cannot say. A finite entry in it evidences positivity for that cell – somebody was observed there – but does not prove it, and an empty cell may be empty by chance or impossible by design; the table cannot tell the two apart.

A bare plot() draws the same information as a curated four-panel overview: the pair classification as the anchor, realised availability, representation, and a small support graph.

plot(audit)

How to read this figure. Start at the classification matrix on the left – it is the precise panel, and the other three explain it. Each cell is one pair of blocks, and its colour is the pair’s class and nothing else: blue means the pair is carried directly by participants observed in both blocks, amber that it is not, and travels instead through other blocks. Colour never encodes how much support a pair has, which is why a thinly supported blue cell looks exactly like a richly supported one – the amount is printed inside the cell. There is no traffic-light scale and no good/bad axis here: these are properties of a design, not verdicts on it. The small triangles and squares in the cell corners are advisories, the top-right panel shows how much of each block each declared stratum actually has, the bottom-right panel how much weight is behind it, and the little graph is the same design drawn as a network. Read as a whole, this composite tells you which relationships your fit will rest on directly, which will rest on a model, and where the design is thin.

Each panel is also available on its own, and comes back as a ggplot you can restyle or save. The matrix first, full size:

plot(audit, type = "matrix")

How to read this figure. Read the lower triangle cell by cell; the diagonal just carries each block’s own number of participants. A blue cell prints \(n_{jk}\), the number of participants observed in both blocks – that pair is carried by real joint observations. An amber cell prints B=, and that is the different and more demanding case: no participant carries both blocks, so the relationship is reachable only by travelling through others, and it is identified conditional on the declared rank-L model and its required structural conditions, which are not empirically checked by this audit. The B value is the weakest link of the route it travels through – the smallest joint count on any step of that path – so B=25 means the entire relationship leans on 25 participants somewhere along the way. A light grey cell means the pair has neither, and a darker grey cell with an X means no participant anywhere carries both blocks. The corner glyphs are advisories: a triangle where the support the pair leans on falls below the declared ess_min line, a square where a single edge carries the whole connection. What this decides for you: amber and grey cells are the relationships to be careful about when you read the fit, and a triangle tells you the number behind a cell is small even when its colour is the strongest one.

The graph draws the same design as a network, which makes its shape obvious at a glance.

plot(audit, type = "graph")

How to read this figure. Look at the node colours first: they are the connected components, so blocks sharing a colour can reach each other, through one step or several, and blocks of different colours cannot be related by this design at all. Then look at the edges. An edge is a directly supported pair; its width is the raw count of participants observed in both blocks and nothing else – never a strength, a quality or a confidence. A dashed purple edge is a single bridge: remove those participants and the two sides stop being connected. A solid blue edge has at least one alternative route around it. That distinction is design-support redundancy – how many independent ways the design has of connecting two blocks – and it is a statement about the shape of the design, not about how well anything is estimated. What this decides for you: a single bridge is the place where losing data would cost you most.

When a pair is reachable only through others, pair = draws the route it depends on.

plot(audit, type = "graph", pair = c("methylation", "metabolites"))

How to read this figure. The emphasised path is the route the audit recorded for this pair, and the faded elements are everything the route does not use. The label on one edge marks the bottleneck: the weakest step of the route, and therefore the number the whole relationship rests on. The subtitle carries the condition and is not decoration – a model-transfer-eligible pair is identified conditional on the declared rank-L model and its required structural conditions, which are not empirically checked by this audit, and if the audit recorded other equally wide routes the subtitle says so rather than letting this one look privileged. What this decides for you: whether you are willing to lean on that route, knowing which participants and which assumption it rests on.

Availability shows the same design from the participants’ side.

plot(audit, type = "availability")

How to read this figure. One row per declared stratum, one column per block, and the colour is the proportion of that stratum present in that block. The scale is pinned to the full range from 0 to 1, never stretched to the range that happens to occur, so a pale column is genuinely sparse rather than merely the palest thing in this particular design. The number in each cell is the count behind the proportion. What this decides for you: a column that is pale in one stratum and dark in another is a block your cohorts do not contribute to equally. Note also what the table cannot say. A filled cell evidences positivity for that combination – somebody was observed there – but does not prove it, and an empty cell may be empty by chance or impossible by design; this table cannot tell those two apart.

Representation asks the same question of the declared weights.

plot(audit, type = "representation")

How to read this figure. Two panels, the same rows in the same order in both, and deliberately not one shared axis: a count and a proportion do not live on the same scale. On the left, the effective sample size – what the declared weights leave of the raw count once their spread is taken into account – annotated with the size of the stratum it came from, because comparing effective counts across strata of different sizes is exactly the mistake this annotation exists to prevent. On the right, the share of each stratum present in each block. Bar colour is the declared stratum, the same colour in both panels, so a stratum can be followed across them. The dashed grey lines are the declared design lines: conventional lines set by the analyst, carrying margin, configurable, and not thresholds at which anything fails. A bar short of one is worth a look, not an alarm – and how far short it falls is visible precisely because everything here is reported as a value.

Finally, the question users of block-wise missing data ask first: which combinations of blocks do participants actually carry?

plot(audit, type = "patterns")

How to read this figure. Each row is one combination of blocks that somebody actually carries. On the left, filled dots mark which blocks that combination contains; on the right, the bar and its number count the participants who carry exactly that combination. Rows run most frequent first. What this decides for you: this is the design in its rawest form, and it usually explains the rest of the audit at once. Here there are two rows and no others – 60 participants carry methylation, expression and proteins, 25 carry proteins and metabolites – and that is the whole reason metabolites hangs off a single edge and reaches the other blocks only through proteins.

7.1 Auditing the object you are about to fit

mgcca_audit() also takes the data object itself – a list of blocks, a MultiAssayExperiment, or an HDF5 file – and derives the availability table through the package’s own readers, with the same union-of-IDs rule the estimator uses. Only the participant IDs are read; no block is materialised. Auditing the very object the fit will consume is what makes the audit a statement about the actual analysis.

On an HDF5-backed design, save_hdf5 = TRUE writes the audit into an availability_audit group of that same file, so the certificate travels with the data, and mgcca_audit_load() brings it back whole.

set.seed(1)
blocks <- list(
  methylation = matrix(rnorm(60 * 4), 60, 4,
                       dimnames = list(ids[core], paste0("cg", 1:4))),
  expression  = matrix(rnorm(60 * 3), 60, 3,
                       dimnames = list(ids[core], paste0("gene", 1:3))),
  proteins    = matrix(rnorm(85 * 3), 85, 3,
                       dimnames = list(ids, paste0("prot", 1:3))),
  metabolites = matrix(rnorm(25 * 2), 25, 2,
                       dimnames = list(ids[extra], paste0("met", 1:2))))

h5_audit <- tempfile(fileext = ".h5")
desc <- mgcca_import_hdf5(blocks, filename = h5_audit, overwriteFile = TRUE)

stored <- mgcca_audit(desc, L = 2, strata = cohort, weights = avail_w,
                      save_hdf5 = TRUE)

# the availability table derived from the blocks is the one built by hand
identical(unname(stored$table), unname(avail))

## [1] TRUE

# and the stored audit comes back whole
identical(mgcca_audit_load(h5_audit), stored)

## [1] TRUE

What is written is structured rather than opaque: the tables go in as datasets, block and participant names as their dimnames, the status codes as integers with the frozen vocabulary stored beside them as a legend, and the transfer routes as sequences of block indices. Any HDF5 browser can read it.

stored$saved$tables

##  [1] "table"                   "strata"                 
##  [3] "weights"                 "strata_levels"          
##  [5] "pair_codes"              "advisory_codes"         
##  [7] "adjacency"               "kappa"                  
##  [9] "pairs"                   "paths"                  
## [11] "representation"          "availability_by_block"  
## [13] "availability_by_stratum" "pairs_by_stratum"       
## [15] "graph"                   "design"

rownames(as.matrix(
  BigDataStatMeth::hdf5_matrix(h5_audit, "availability_audit/pair_codes")))

## [1] "DIRECT_SUPPORT"          "MODEL_TRANSFER_ELIGIBLE"
## [3] "UNSUPPORTED"             "STRUCTURAL_ZERO"

BigDataStatMeth::hdf5_close_all()
unlink(h5_audit)

8 Appendix: the full-scale TCGA walkthrough (not evaluated)

Caution: nothing in this appendix is evaluated when the vignette is built. The code is shown so it can be copied and run; the outputs and figures beneath it were recorded from a run on the full data and are reproduced verbatim, not recomputed here. Recorded output can drift from the code that produced it, which is why the live account of the same workflow is in section 6 above and this one is kept separate. Running it needs curatedTCGAData, TCGAutils, TxDb.Hsapiens.UCSC.hg19.knownGene and org.Hs.eg.db (all declared in Suggests), and downloads several hundred MB.

Adrenocortical carcinoma (ACC) is a rare, aggressive endocrine cancer. The TCGA ACC cohort provides, for 80 patients, two omics layers – RNA-seq gene expression (20,501 genes) and DNA methylation (485,577 CpGs) – together with rich clinical annotation (tumour stage, survival, treatment, …). Both tables have p >> n, so both take the Gram route and no variable-by-variable matrix is ever formed; the fit completes in about 30 seconds.

8.1 Download, impute and assemble the tables

library(curatedTCGAData); library(TCGAutils); library(MultiAssayExperiment)

mae <- curatedTCGAData("ACC", c("RNASeq2GeneNorm", "Methylation"),
                       dry.run = FALSE, version = "1.1.38")
mae <- mgccaImpute(mae, method = "knn", k = 10, rowmax = 0.5, colmax = 0.8,
                   remove.col = TRUE, impute.zero = TRUE)

# individuals x variables, aligned by the 12-char sample ID; methylation -> numeric
tables_list <- getTables(mae)
for (a in seq_along(tables_list))
  rownames(tables_list[[a]]) <- substr(rownames(tables_list[[a]]), 1, 12)
tables_list <- matrixChr2Num(tables_list, index = 1)

8.2 Integrate the two layers

sel <- mgcca_select_lambda(tables_list, nfac = 2, K = 5, seed = 1)

fit <- mgcca(tables_list, filename = tempfile(fileext = ".h5"), nfac = 2,
             method = "penalized", lambda = sel$lambda,
             scores = TRUE, route = "auto", collect = TRUE)

attr(fit, "desc")$route_dual   # both TRUE: methylation and RNA-seq both p >> n
summary(fit, top = 4)

Recorded summary() output of that run, which used lambda = c(0.75, 1):

Generalized Canonical Correlation Analysis (mgcca)
  tables      : 2  (ACC_RNASeq2GeneNorm-20160128, ACC_Methylation-20160128)
  individuals : 80
  components  : 2
  eigenvalues : 1, 1
  AVE (outer) : 0.04473, 0.04128
  scores      : available (per table)

AVE per table (outer model):
                              comp1  comp2
ACC_RNASeq2GeneNorm-20160128 0.0847 0.0692
ACC_Methylation-20160128     0.0430 0.0401

Top 4 variables per component (correlation with Y):
  [ACC_RNASeq2GeneNorm-20160128]
    comp1: ARHGDIB (-0.79), TBX3 (+0.74), TNFAIP8L2 (-0.73), FLI1 (-0.73)
    comp2: STK25 (-0.81), PRCC (-0.80), ARPC2 (-0.79), ZDHHC16 (-0.77)
  [ACC_Methylation-20160128]
    comp1: cg14278148 (-0.74), cg12653626 (-0.71), cg05796838 (-0.71), cg27641801 (-0.70)
    comp2: cg10815272 (-0.74), cg06495763 (+0.72), cg10469906 (-0.70), cg13301931 (-0.70)

8.3 Do the shared components track the disease?

Colouring the individuals by clinical group shows that the integrated axes are not arbitrary. Tumour stage grades along the components (early-stage ellipse up-left, stage IV down-right), and vital status separates survivors from deceased patients.

stage <- factor(mae$pathologic_stage)
vital <- factor(ifelse(mae$vital_status == "1", "deceased", "alive"))
plotIndividuals(fit, group = stage)   # figure 1
plotIndividuals(fit, group = vital)   # figure 2

mgcca_associate(fit, list(
  stage  = factor(mae$pathologic_stage),
  vital  = factor(mae$vital_status),
  gender = factor(mae$gender)))

Recorded output:

  phenotype component     R2       p   n   p.adj
1     stage         1 0.113   0.030  78   0.045
2     stage         2 0.154 5.9e-03  78   0.018
3     vital         1 0.064   0.023  80   0.045
4     vital         2 0.287 3.1e-07  80   1.9e-06
5    gender         1 0.031   0.115  80   0.138
6    gender         2 0.003   0.630  80   0.630

Both components are significantly associated with tumour stage, and the second component is strongly associated with survival, while gender is not – the shared signal reflects disease biology, not a technical or demographic artefact.

8.4 What drives the components

The methylation correlation circle shows the full cloud of 485k CpGs with the top-correlated ones labelled; plotLoadings() lists the strongest CpGs for the first component; and plotAVE() shows how much of each layer each component captures.

plotVariables(fit, table = "ACC_Methylation-20160128", top = 15)
plotLoadings(fit,  table = "ACC_Methylation-20160128", comp = 1, top = 15)
plotAVE(fit)

If you only need the results on disk (no in-memory object), drop collect = TRUE and read them back later with mgcca_results(h5file, datasets = names(tables_list)).

8.5 Subsetting by genomic ranges

genCoords() restricts a MultiAssayExperiment to the features falling in a genomic region. The region is a GRanges, so GenomicRanges has to be attached (it comes with TCGAutils); the methylation assay is identified by meth.index, and col.name names the rowData column holding its coordinates. As an example, features around the MEN1 locus, a tumour-suppressor gene at 11q13 related to this cancer:

library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)

# MEN1 +- 2kb
chr.coords <- GRanges("11:64550986-64578766")

# Annotate RNA-seq with ranges; methylation already carries CpG coordinates
mae.annotated <- TCGAutils::symbolsToRanges(mae)
mae.annotated[[3]] <- NULL                       # drop the unranged part

mae.subset <- genCoords(multiassayexperiment = mae.annotated,
                        gen.coords = chr.coords,
                        meth.index = 1, col.name = "Genomic_Coordinate")

mae.subset.imputed <- mgccaImpute(mae.subset, method = "knn", k = 10,
                                  rowmax = 0.5, colmax = 0.8,
                                  remove.col = TRUE, impute.zero = TRUE)

sub_tables <- getTables(mae.subset.imputed)
for (a in seq_along(sub_tables))
  rownames(sub_tables[[a]]) <- substr(rownames(sub_tables[[a]]), 1, 12)
sub_tables <- matrixChr2Num(sub_tables, index = 1)

sub_sel <- mgcca_select_lambda(sub_tables, nfac = 2, K = 5, seed = 1)
fit.region <- mgcca(sub_tables, filename = tempfile(fileext = ".h5"),
                    nfac = 2, method = "penalized", lambda = sub_sel$lambda,
                    scores = TRUE, collect = TRUE)
plotIndividuals(fit.region, group = factor(mae$vital_status))

9 Session information

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] MultiAssayExperiment_1.39.1 SummarizedExperiment_1.43.0
##  [3] Biobase_2.73.2              GenomicRanges_1.65.4       
##  [5] Seqinfo_1.3.2               IRanges_2.47.5             
##  [7] S4Vectors_0.51.10           BiocGenerics_0.59.13       
##  [9] generics_0.1.4              MatrixGenerics_1.25.0      
## [11] matrixStats_1.5.0           mgcca_0.99.7               
## [13] BiocStyle_2.41.0           
## 
## loaded via a namespace (and not attached):
##  [1] made4_1.87.0          gtable_0.3.6          impute_1.87.0        
##  [4] xfun_0.62             bslib_0.12.0          ggplot2_4.0.3        
##  [7] caTools_1.18.4        lattice_0.23-1        vctrs_0.7.3          
## [10] tools_4.6.1           bitops_1.1-0          tibble_3.3.1         
## [13] BiocBaseUtils_1.15.1  pkgconfig_2.0.3       KernSmooth_2.23-27   
## [16] Matrix_1.7-6          data.table_1.18.6.1   RColorBrewer_1.1-3   
## [19] S7_0.2.2              scatterplot3d_0.3-45  lifecycle_1.0.5      
## [22] compiler_4.6.1        farver_2.1.2          gplots_3.3.0         
## [25] tinytex_0.61          htmltools_0.5.9       sass_0.4.10          
## [28] RCurl_1.98-1.20       yaml_2.3.12           pillar_1.11.1        
## [31] jquerylib_0.1.4       MASS_7.3-66           DelayedArray_0.39.8  
## [34] cachem_1.1.0          magick_2.9.1          abind_1.4-8          
## [37] gtools_3.9.5          tidyselect_1.2.1      digest_0.6.39        
## [40] dplyr_1.2.1           purrr_1.2.2           bookdown_0.48        
## [43] labeling_0.4.3        ggthemes_7.0.0        ade4_1.7-24          
## [46] fastmap_1.2.0         grid_4.6.1            cli_3.6.6            
## [49] SparseArray_1.13.4    magrittr_2.0.5        S4Arrays_1.13.2      
## [52] dichromat_2.0-1       withr_3.0.3           scales_1.4.0         
## [55] rmarkdown_2.32        XVector_0.53.0        otel_0.3.0           
## [58] evaluate_1.0.5        knitr_1.52            rlang_1.3.0          
## [61] Rcpp_1.1.2            BigDataStatMeth_2.0.6 glue_1.8.1           
## [64] BiocManager_1.30.27   jsonlite_2.0.0        R6_2.6.1