Generalized Canonical Correlation with missing individuals - Examples
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.
Installation
mgcca is installed from Bioconductor:
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## [1] 646Running 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:
## 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)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):
## 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.053Cholesterol (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
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.
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.
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).
Variance explained
plotAVE() shows how much of each table’s variance each
shared component captures (outer-model AVE); plotScree()
shows the component eigenvalues.
Significant variables
getSignif() lists the variables significantly correlated
with any component and topVars() returns the strongest
ones.
## variable table
## 1 cg1 methylation
## 2 cg2 methylation
## 3 cg3 methylation
## 4 cg4 methylation
## 5 cg5 methylation
## 6 cg6 methylation## $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## TC LDL HDL BMI
## 2.15e-44 9.87e-34 6.95e-28 1.31e-01The 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.
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:
## [1] "solve"## [1] "auto"## methylation clinical other
## FALSE FALSE FALSE## [1] 0.7271469 0.6392769All the plotting and post-analysis helpers work on the reloaded object exactly as before:
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:
## [1] 646 2## [1] TRUEoutputs 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.
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
BigDataStatMethreading 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)mgccauses the Gram route (an n x n matrix), never the p x p covariance matrix, which would be intractable. Forp < nit 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.
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 filesAssembling 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.
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).
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.
## [1] "ListMAE"## RNASeq2GeneNorm RPPAArray
## [1,] 79 46
## [2,] 198 33for (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## [1] 46So 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:
## [1] "1487.0317"## [1] 1487.032Choosing 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.
## 9 of 9 candidates are within one standard error of the best; the penalty is not sharply identified. Report the region.## 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.
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)## 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:
## [1] 79 2## [1] "Y" "AVE"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
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.
## 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-01The 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.
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## [1] 0.7263080 0.8385205Nine 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.
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:
## $RPPAArray
## comp1 comp2
## TCGA-OR-A5J2 0.1359425 -0.1493703
## TCGA-OR-A5J3 -0.2464442 -0.7029714
## TCGA-OR-A5J6 1.0704808 0.8479712The 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.
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 25L 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.
## 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:
## 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 25The 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.
## block1 block2 n_jk s_e s_e_ess kappa is_bridge
## 4 proteins metabolites 25 25 25 1 TRUE## methylation expression proteins metabolites
## FALSE FALSE TRUE FALSEsummary() 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.
## 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.
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:
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.
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.
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.
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.
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?
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.
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## [1] TRUEWhat 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.
## [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"## [1] "DIRECT_SUPPORT" "MODEL_TRANSFER_ELIGIBLE"
## [3] "UNSUPPORTED" "STRUCTURAL_ZERO"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 @ref(live-miniacc) 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.
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)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)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.630Both 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.
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)).
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))Session information
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: Etc/UTC
## 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.12
## [9] generics_0.1.4 MatrixGenerics_1.25.0
## [11] matrixStats_1.5.0 mgcca_0.99.6
## [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.61 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 stringr_1.6.0
## [25] gplots_3.3.0 htmltools_0.5.9 sys_3.4.3
## [28] buildtools_1.0.0 sass_0.4.10 RCurl_1.98-1.20
## [31] yaml_2.3.12 pillar_1.11.1 jquerylib_0.1.4
## [34] MASS_7.3-66 DelayedArray_0.39.8 cachem_1.1.0
## [37] abind_1.4-8 gtools_3.9.5 tidyselect_1.2.1
## [40] digest_0.6.39 stringi_1.8.9 dplyr_1.2.1
## [43] purrr_1.2.2 labeling_0.4.3 maketools_1.3.2
## [46] ggthemes_6.0.0 ade4_1.7-24 fastmap_1.2.0
## [49] grid_4.6.1 cli_3.6.6 SparseArray_1.13.4
## [52] magrittr_2.0.5 S4Arrays_1.13.2 withr_3.0.3
## [55] scales_1.4.0 rmarkdown_2.32 XVector_0.53.0
## [58] otel_0.2.0 evaluate_1.0.5 knitr_1.52
## [61] rlang_1.3.0 Rcpp_1.1.2 BigDataStatMeth_2.0.6
## [64] glue_1.8.1 BiocManager_1.30.27 jsonlite_2.0.0
## [67] R6_2.6.1