---
title: "Data Integration and Dimension Reduction with RGCCA"
author: "Vincent Guillemot"
format:
  html:
    code-fold: true
    code-tools: true
    code-copy: true
    code-block-border-left: true
    toc: true
    toc_depth: 2
    toc_float: true
    number_sections: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

library(dplyr)
library(tidyr)
library(ggplot2)
library(RGCCA)
library(fgsea)
library(ggbeeswarm)
library(patchwork)
# library(HGNChelper)
library(ggpubr)

theme_set(theme_bw(base_size = 13))

shared_doc <- "https://docs.google.com/document/d/1T1Sm1V6IQ8hX2TXz84b6AQ2HgSsYFFIJktZBUQTN1ww/edit?usp=sharing"

mini_corrplot <- function(mat) {
  if (is.null(dimnames(mat))) colnames(mat) <- sprintf("V%i", 1:ncol(mat))
  rownames(mat) <- colnames(mat)
  as_tibble(mat, rownames = "row") %>%
    mutate(row = factor(row, levels = rev(unique(row)))) %>%
    pivot_longer(-row, names_to = "col") %>%
    ggplot(aes(col, row, fill = factor(value))) +
    geom_tile(color = "white", show.legend = FALSE) +
    geom_text(aes(label = round(value, 2)), size = 6, color = "white") +
    scale_fill_brewer(palette = "Set2") +
    coord_equal() +
    labs(x = "", y = "")
}
```

# Before you start

This notebook complements two other workshop files:

  * `reminders.qmd` introduces the core statistical ideas needed before touching RGCCA, 
  * `more_theory.qmd` gives a compact theory-oriented overview of RGCCA and SGCCA, with parameters and practical links to the package.

# Introduction

## Goals

The goals of this part of the workshop are to:

- identify the kinds of biological questions that can be addressed with RGCCA and the kinds that require other methods,
- implement RGCCA on simple and realistic data integration problems,
- understand the main parameters of the method,
- learn how to choose reasonable parameter values,
- interpret the main graphical outputs,
- and use sparse extensions to highlight variables of interest.

## Why use an R Markdown notebook?

A notebook is portable, keeps code and interpretation in one place, and supports reproducible analyses.

I decided to use the `qmd` for the notebooks of the RGCCA workshop, but for all our uses and purposes it's almost like an "R Markdown" document.

## Types of examples used here

This document relies on two types of examples:

- small simulated datasets where the correlation structure is fully controlled,
- and TCGA breast cancer data, which are also used in other workshop sessions.

# The big picture

::: {.callout-tip appearance="simple" icon="true"}


RGCCA is a method for integrating several data blocks measured on the **same samples**.  
Each block can represent a different type of information, such as gene expression, methylation, proteomics, clinical variables, or any other set of features collected on the same individuals.

:::::: {.columns}

::: {.column}

The goal is to build a few **synthetic components** per block that:

- summarize the main variation within each block,
- capture the relationships between connected blocks,
- and help identify the variables that drive these shared patterns.

In practice, an RGCCA analysis requires four main ingredients:

- **Blocks**: the datasets to integrate;
- **A design matrix**: which tells the model which blocks are connected;
- **A mode of analysis**: exploratory or driven by a response block;
- **Regularization or sparsity parameters**: which control stability and variable selection.

:::

::: {.column}

A useful way to think about RGCCA is the following:

- **PCA** summarizes one block,
- **CCA** relates two blocks,
- **RGCCA** extends this logic to several connected blocks,
- **SGCCA** adds variable selection to focus on the most informative features.

At the end of the analysis, the main outputs are:

- block components for the samples,
- relationships between blocks,
- variable loadings or selected features,
- and graphical displays that help interpret the biological structure of the data.

:::

::::::

:::


# A quick tour of the package

The `RGCCA` package is available on CRAN.

```r
install.packages("RGCCA")
```

A development version with additional methods is available on GitHub. Unless you already feel comfortable with RGCCA, you don't need to go on GitHub.

```r
# remotes::install_github("rgcca-factory/RGCCA")
```

The main function is called `rgcca()`. Before using it on real data, it is worth reading the help page and the vignette.

```r
?RGCCA::rgcca
```

The examples that we will see while discussing will also be recorded in the [Shared workshop document](`r shared_doc`).

Below is the example from the manual, on Russett's data (from "Russett BM. Inequality and Instability: The Relation of Land Tenure to Politics. World Politics. 1964;16(3):442-454. doi:10.2307/2009581"). Run this example to check that everything works!

```{r russett}
data(Russett)
blocks <- list(
  agriculture = Russett[, seq(3)],
  industry = Russett[, 4:5],
  politic = Russett[, 6:11]
)

politic <- as.factor(apply(Russett[, 9:11], 1, which.max))

# RGCCA with default values : Blocks are fully connected, factorial scheme
# tau = 1 for all blocks, one component per block.
fit_rgcca <- rgcca(blocks = blocks)

print(fit_rgcca)

plot(fit_rgcca, type = "weight", block = 1:3)

plot(fit_rgcca,
  type = "sample", block = 1:2,
  comp = rep(1, 2), resp = politic
)

```

In a few words, RGCCA is designed for **block-structured data**:

- a block is a coherent set of variables,
- all blocks are measured on the same samples,
- and blocks can be connected according to the scientific question.

Below is an example of design with 4 blocks of data (`A`, `B`, `P`, and `R`)[^1]. The goal of RGCCA is to summarize each block by components (that I noted `YA`, `Y_B`, `Y_P`, and `Y_R`). The links between the blocks represent the "design matrix", meaning the correlations between components that I want to impose. With this design, I'm telling RGCCA to maximize all the correlations between pairs of components except two (A -- R and B -- P). These links represent the "inner" model.

![](img/RGCCA_design.png)

[^1]: notice that all blocks have the same height (i.e. same number of observations), but different widths (i.e. different number of variables).

One of the strengths of RGCCA is that it places unsupervised and supervised multiblock approaches in the same general framework. In practice, this is useful when one analysis needs to combine exploratory integration, prediction-oriented modeling, and feature interpretation.



## Practical questions to keep in mind

Before running any model, ask the following:

- What is a biologically meaningful block in my study (e.g. containing groups of samples)?
- Are all blocks measured on exactly the same samples?
- Should all blocks be connected, or only some of them?
- Is my goal exploration, prediction, feature selection, or a mix of these?
- Do I need all variables, or do I want a sparse solution?

# Types of questions addressed in the workshop

In the [shared document](`r shared_doc`), I prepared a table of anonymized biological questions collected from a previous RGCCA workshop.

Your turn: based on these examples, formulate your own data integration question, even if you do not have your data with you. Throughout the workshop, we will return to these questions and connect them to concrete analytical choices.

A useful habit is to express each question in terms of:

- **samples**,
- **blocks**,
- **target outcome or response**, if any,
- **expected relationships between blocks**,
- and **expected outputs**, such as components, selected variables, or interpretable pathways.

# The parent method: Principal Component Analysis

Many questions answered with PCA are exploratory and unsupervised.

| Typical question | PCA output or visualization |
|------------------|-----------------------------|
| Are my data informative? | Scree plot |
| Is there a batch effect? | Sample map |
| Are groups of variables correlated? | Correlation circle |

PCA is the right starting point because RGCCA also builds low-dimensional components. The main difference is that PCA summarizes one block at a time, whereas RGCCA summarizes several blocks jointly while taking their relationships into account.

## A first toy example

Assume that we are performing (standardized) PCA on a dataset with 3 variables. The only information that I am giving you is the correlation matrix below. Based on this correlation matrix, what's the variance of the first component and what are the weights[^2]?

[^2]: recall that a component in PCA is a linear combination (which needs weights) of the variables

```{r corrplot easy}
mini_corrplot(matrix(0.8, 3, 3) + diag(0.2, 3))
```

Same exercise with this second matrix.

```{r corrplot less easy}
mini_corrplot(matrix(c(1, 0.8, 0, 0.8, 1, 0, 0, 0, 1), 3, 3))
```

> Bonus: in R, use `eigen()` to decompose these correlation matrices and check that your intuition was correct. PCA is based on the eigenvalue decomposition of the covariance or correlation matrix.

```{r solution eigen, collapse = TRUE}
eigen(matrix(0.8, 3, 3) + diag(0.2, 3))
eigen(matrix(c(1, 0.8, 0, 0.8, 1, 0, 0, 0, 1), 3, 3))
```

## Key intuition to carry into RGCCA

PCA teaches several ideas that remain essential for RGCCA:

- highly correlated variables tend to move together,
- a component is a weighted combination (summary) of variables,
- the first component captures the strongest dominant pattern,
- variables unrelated to the dominant pattern usually receive little weight,
- and visualization of samples and variables often reveals structure faster than a table of coefficients.

# Two or three tiny blocks

Let's simulate easy to handle data, on which we control everything, and ask a few questions to `RGCCA`.

## One block with several variables and one block with one variable

Suppose the data are organized into two blocks:

- a first block with variables `X1`, `X2`, `X3`, and `X4`,
- a second block with a single variable `V`.

Assume the variables follow the correlation structure below.

```{r xv}
R_xv <- matrix(
  data = c(
    1, 0.8, 0.8, 0, -0.85,
    0.8, 1, 0.8, 0, -0.85,
    0.8, 0.8, 1, 0, -0.85,
    0, 0, 0, 1, 0,
    -0.85, -0.85, -0.85, 0, 1),
  nrow = 5,
  ncol = 5,
  dimnames = list(
    c("X1", "X2", "X3", "X4", "V"),
    c("X1", "X2", "X3", "X4", "V")))

mini_corrplot(R_xv)
```

My question is the following: what is the best combination of `X1`, `X2` and `X3` that is the most correlated with `V`?

RGCCA is well suited to answer this question because, like PCA, it creates components with lower dimension than the original blocks. But, unlike PCA, it tries to balance two goals:

1. each component should (optimally) summarize its own block,
2. and the components that we decide to connect should be as strongly related as possible.

Because RGCCA builds components, it also computes variable weights. In this example, we would expect:

- `X1`, `X2`, and `X3` to receive large weights of similar magnitude,
- `X4` to receive a weight close to zero,
- and `V` to define the second block component by itself.

Let us test our intuition with simulated (Gaussian multivariate) data.
To simulate this data, I use `MASS:mvrnorm`: this function is very similar to `rnorm`.

Below is a simple example showing the histogram of a random univariate gaussian sample ($\mathcal N(0, 1)$) of size 100.

```{r rnorm histogram}
tibble(x = rnorm(100)) %>%
  ggplot(aes(x)) + 
  geom_histogram(
    bins = 10, 
    fill = "#2b9b81",
    color = "white")
```

A multivariate sample needs a vector of theoretical means and a theoretical covariance matrix. Here we will generate a sample of size $n = 27$, of a $p = 5$-dimensional multivariate Gaussian with 0 mean and a covariance matrix that we previously defined. Remember your STATS courses: we could denote this multivariate variable like below.

$$
X = \mathcal{N}
  \left(
    \mu = \left[
      \begin{array}{c} 
        0\\0\\0\\0
      \end{array}
    \right],
    \Sigma = \left[
      \begin{array}{ccccc} 
        1 & 0.8 & 0.8 & 0 & -0.85\\
        0.8 & 1 & 0.8 & 0 & -0.85\\
        0.8 & 0.8 & 1 & 0 & -0.85\\
        0 & 0 & 0 & 1 & 0\\
        -0.85 & -0.85 & -0.85 & 0 & 1 \\
      \end{array}
      \right]
  \right)
$$


```{r two mini blocks}
set.seed(413)
dat_xv <- MASS::mvrnorm(n = 27, mu = rep(0, 5), Sigma = R_xv, empirical = TRUE)
blocks <- list(X = dat_xv[, 1:4], V = dat_xv[, 5, drop = FALSE])
```

Running RGCCA on this example requires at least one object: the list of "blocks". On your real dataset, you would start at defining this list with the blocks of data you want to analyze with RGCCA.

```{r rgcca two mini blocks}
fit.rgcca <- rgcca(blocks = blocks)
```

We will take a look at the following graphs:

- **AVE plot**: how much of the relevant structure is captured by the components?
- **Sample plot**: do the two block components show a clear relationship?
- **Loadings plot**: which variables drive each block component?

::: {.panel-tabset}

### AVE plot

The AVE shows that the first component of `V` explains its own block very well. I's not surprinsing because `V` contains only one variable.

It also shows that the first component of `X` explains $65\%$ of its own block, which is not unexpected, since its based on 3 out of its 4 variables, the component will explain at most $75\%$.

```{r}
plot(fit.rgcca, type = "ave")
```

### Samples

Here we see that the correlation between the two first components is negative. It's directly related to the correlation of $0.85$ between `X1` and `V`, `X2` and `V`, and `X3` and `V`.

```{r}
plot(fit.rgcca, type = "samples")
```

### Loadings

Finally, the loadings are as expected: 

  * high similar loadings for `X1`, `X2`, and `X3`,
  * 0 for `X4`,
  * and 1 for `V`.

```{r}
plot(fit.rgcca, type = "loadings")
```

:::

## Two simple blocks with two variables

Consider the new correlation matrix below for two blocks with two variables each: block `A` with `A1` and `A2`, and block `B` with `B1` and `B2`.

Before running the model, predict which variables will dominate the first component of each block. After that, simulate some data, run `rgcca` and test that your intuition was correct!

```{r ab}
R_ab <- matrix(
  data = c(
    1, 0.8, 0.9, 0,
    0.8, 1, 0.9, 0,
    0.9, 0.9, 1, 0,
    0, 0, 0, 1),
  nrow = 4,
  ncol = 4,
  dimnames = list(
    c("A1", "A2", "B1", "B2"),
    c("A1", "A2", "B1", "B2")))

mini_corrplot(R_ab)
```

*Sploiler Alert*: Don't look below if you want to do it yourself!


```{r rgcca two mini blocks exo}
set.seed(413)

dat_ab <- MASS::mvrnorm(n = 27, mu = rep(0, 4), Sigma = R_ab, empirical = TRUE)
blocks_ab <- list(A = dat_ab[, 1:2], B = dat_ab[, 3:4])
fit.rgcca_ab <- RGCCA::rgcca(blocks = blocks_ab)

plot(fit.rgcca_ab, type = "ave")
plot(fit.rgcca_ab, type = "samples")
plot(fit.rgcca_ab, type = "loadings")
```

## Three tiny blocks, with one block of noise

Let's add a block to the mix (containing only noise).

```{r}
Z <- matrix(
  data = rnorm(27 * 3), 
  nrow = 27, ncol = 3, 
  dimnames = list(1:27, sprintf("Z%i", 1:3)))
```

Run RGCCA on three blocks, with the default design (all blocks are connected).

```{r}
blocks_abz <- list(A = dat_ab[, 1:2], B = dat_ab[, 3:4], Z = Z)
fit.rgcca_abz <- RGCCA::rgcca(blocks = blocks_abz)

plot(fit.rgcca_abz, type = "ave")
plot(fit.rgcca_abz, type = "samples", block = 1:2)
plot(fit.rgcca_abz, type = "samples", block = c(1, 3))
plot(fit.rgcca_abz, type = "samples", block = c(2, 3))
plot(fit.rgcca_abz, type = "loadings")

```

Let's use a fancy design, where blocks A and Z are disconnected. The difference is not striking, because we see what we want: the strong correlation between `A` and `B`.

```{r}
C <- matrix(
  data = c(
    0, 1, 0,
    1, 0, 1,
    0, 1, 0), 3, 3)

fit.rgcca_abz_fancy <- RGCCA::rgcca(blocks = blocks_abz, connection = C)

plot(fit.rgcca_abz_fancy, type = "ave")
plot(fit.rgcca_abz_fancy, type = "samples", block = 1:2)
plot(fit.rgcca_abz_fancy, type = "samples", block = c(1, 3))
plot(fit.rgcca_abz_fancy, type = "samples", block = c(2, 3))
plot(fit.rgcca_abz_fancy, type = "loadings")

```

The message of this exercice on the design matrix is that if changing the design
has a very noticeable effect on the result, you should really dive into the data
to understand why.


# Understanding the main outputs

## Average Variance Explained and related summaries

As in PCA, RGCCA provides several ways to summarize what the extracted components capture. Three quantities are especially useful:

- **outer AVE**, which reflects how well the components summarize their own blocks,
- **inner AVE**, which reflects how strongly connected components are related to one another,
- **outer variances**, which are closer to the familiar variance summaries used in PCA.

These summaries are helpful because a good model is not only one that finds correlated components. It should also produce components that remain meaningful within each block.

## A simple rule for interpretation

- High **outer AVE** means strong within-block summarization.
- High **inner AVE** means strong between-block integration.
- A useful model usually balances both according to the scientific objective.

# Three less small blocks and more on the role of design

The design matrix is one of the most important ideas in RGCCA. It encodes which blocks are assumed to be connected.

This is not just a technical detail. It is how the scientific question enters the model.

## A larger simulated example

In this example, we have:

- 150 samples split into two known groups,
- block `A` with 7 variables and substantial within-block correlation,
- block `B` with 5 variables and substantial within-block correlation,
- a response block `G` describing group membership,
- variables that truly differ between groups,
- and variables that are correlated but not directly informative about the group structure.

The data are generated by the script `scripts/ABG.R` and stored in the following files:

- `data/simul/A.csv`
- `data/simul/B.csv`
- `data/simul/G.csv`

```{r load simul}
A <- read.csv("../data/simul/A.csv")
B <- read.csv("../data/simul/B.csv")
G <- read.csv("../data/simul/G.csv")
```

After loading the data, take a few minutes to inspect them:

- what are the dimensions of each block,
- do the empirical correlations match your expectations,
- are some variables obviously redundant,
- do the groups already separate in block-specific exploratory plots?

Two questions naturally arise here:

- What are the relationships between the blocks?
- Which variables matter for separating the groups?

Below, the third block is treated as a response block.

A word of caution: when using a response variable, the default design changes. It does not connect all the blocks together anymore. Instead, all the blocks are **only** connected to the response.

```{r rgcca three blocks}
blocks_abg <- list(A = A, B = B, G = G)

fit_rgcca_abg <- rgcca(
  blocks_abg,
  connection = matrix(1, 3, 3),
  ncomp = c(5, 5, 1),
  response = 3
)

plot(fit_rgcca_abg, "ave")

plot(fit_rgcca_abg, comp = 1)
plot(fit_rgcca_abg, comp = 2)
plot(fit_rgcca_abg, comp = 3)
plot(fit_rgcca_abg, type = "samples", response = G$G)
plot(fit_rgcca_abg, type = "samples", block = 1, response = G$G)
plot(fit_rgcca_abg, type = "samples", block = 2, response = G$G)
plot(fit_rgcca_abg, type = "cor_circle", block = 1) + coord_equal()
plot(fit_rgcca_abg, type = "cor_circle", block = 2) + coord_equal()
```

## Interpreting these plots

- The global sample map shows whether the integrated structure aligns with the response.
- Block-specific sample maps show which blocks contribute most to separation.
- Correlation circles help identify variables associated with each component.
- Component-specific plots help distinguish the first dominant pattern from later patterns.

## Important concept: supervised versus unsupervised RGCCA

When a block is declared as a response block, the analysis becomes partly supervised. This changes the interpretation:

- the model is no longer searching for correlated components (!!),
- it is also encouraged to align some components with the response,
- which makes the output more useful for discrimination or prediction.

# Comparing two design matrices

By default, RGCCA estimates components by maximizing relationships between connected blocks. But this may not match your scientific knowledge.

Suppose we do **not** want to force a direct link between blocks `A` and `B`. We can compare two designs:

- a fully connected design,
- and a more selective design in which a noise block is isolated from one of the biological blocks.

To make the contrast visible, we add a noise block `Z`. The graph below shows that the first component of block `Z` is not the same when we change the design. Explore in more details the difference between the two results.

```{r}
Z <- matrix(
  data = round(rnorm(10 * 150), 1),
  nrow = 150,
  ncol = 10,
  dimnames = list(1:150, sprintf("Z%02i", 1:10))
)

blocks_abz <- list(A = A, B = B, Z = Z)

C_comp <- matrix(
  data = 1, nrow = 3, ncol = 3,
  dimnames = list(names(blocks_abz), names(blocks_abz))
)
diag(C_comp) <- 0

C_smart <- matrix(
  data = c(
    0, 1, 1,
    1, 0, 0,
    1, 0, 0
  ),
  nrow = 3,
  ncol = 3,
  dimnames = list(names(blocks_abz), names(blocks_abz))
)

fit_comp <- rgcca(blocks = blocks_abz, connection = C_comp, ncomp = 5)
fit_smart <- rgcca(blocks = blocks_abz, connection = C_smart, ncomp = 5)

data.frame(
  Z1_comp = fit_comp$Y$Z[, 1],
  Z1_smart = fit_smart$Y$Z[, 1],
  Y = G$G) %>%
  ggplot(aes(Z1_comp, Z1_smart, color = Y)) +
  geom_point()

```


RGCCA is not only about extracting components, it is also about *encoding* prior knowledge. A fully connected design can sometimes create relationships that are mathematically convenient but biologically unconvincing. A more selective design may produce components that are less global but more faithful to the study design.

# Real data example: TCGA-BRCA

Now that we have explored key RGCCA features on simulated data, we move to real breast cancer data from The Cancer Genome Atlas.

The blocks include:

- clinical data,
- gene expression,
- microRNA abundance,
- and protein measurements.

```{r}
cli <- read.csv("../data/TCGA/clinical.csv", row.names = 1)
mrna <- t(read.csv("../data/TCGA/mrna.csv", row.names = 1))
colnames(mrna) <- HGNChelper::checkGeneSymbols(
  colnames(mrna), unmapped.as.na = FALSE)$Suggested.Symbol
mirna <- t(read.csv("../data/TCGA/mirna.csv", row.names = 1))
prot <- t(read.csv("../data/TCGA/protein.csv", row.names = 1))
```



## Exploration and quick quality control

Before integrating the blocks, inspect each one separately.

Questions to ask:

- Are sample identifiers aligned across blocks?
- Are there outliers or distributional problems?
- Are some groups strongly unbalanced?
- Do we need filtering, transformation, or scaling?
- Are there missing values that need to be handled first?

> *Spoiler alert*: There are a few propositions below so don't look if you want to do it yourself. We will put the best ideas on the [shared document](`r shared_doc`).

First we check that rownames are the same in all 4 blocks.

```{r}
table(rownames(cli) == rownames(mrna))
table(rownames(cli) == rownames(mirna))
table(rownames(cli) == rownames(prot))
```

To evaluate the balanceness, we count the number of patients per
sub-type.

```{r}
cli %>% count(subtype)
```

Are there any missing values?

```{r}
sum(is.na(cli))
sum(is.na(mrna))
sum(is.na(mirna))
sum(is.na(prot))
```

Below is an example plot showing Gender, Age, and subtype from the clinical block.

```{r qc}
cli %>%
  ggplot(aes(subtype, patient.age_at_initial_pathologic_diagnosis)) +
  geom_violin(aes(color = subtype, fill = subtype), trim = FALSE, alpha = 0.5) +
  geom_beeswarm(aes(shape = patient.gender)) +
  # ggpubr::stat_pwc(hide.ns = TRUE, y = 100) +
  scale_shape_manual(values = c(4, 10)) +
  labs(
    x = "Cancer Subtype",
    y = "Age at Initial Diagnosis",
    fill = "",
    color = "",
    shape = "Gender"
  )
```

Below is the histogram of the values in each of the block.

```{r}
g_hist_mrna <- as_tibble(scale(mrna, scale = TRUE)) %>%
  pivot_longer(everything()) %>%
  ggplot(aes(value)) +
  geom_histogram(bins = 50, color = "white") +
  labs(x = "Gene Expression")

g_hist_mirna <- as_tibble(scale(mirna, scale = TRUE)) %>%
  pivot_longer(everything()) %>%
  ggplot(aes(value)) +
  geom_histogram(bins = 50, color = "white") +
  labs(x = "miRNA Expression")

g_hist_prot <- as_tibble(scale(prot, scale = TRUE)) %>%
  pivot_longer(everything()) %>%
  ggplot(aes(value)) +
  geom_histogram(bins = 50, color = "white") +
  labs(x = "Protein")

(g_hist_mrna + g_hist_mirna) /
  (g_hist_prot + plot_spacer())
```


## First multiblock model

The blocks are combined in a list, like earlier.

```{r}
blocks_tcga <- list(mRNA = mrna, miRNA = mirna, PROT = prot, Subtype = cli$subtype)
```

We ask first for a model predicting the subtype based on the information in the `rRNA`, `miRNA` and `PROT`. Beware that with this model, we do not ask for the shared structure between the omics blocks. Nonetheless, we get a very interesting structure of the component space: we get an interesting separation between the groups (not perfect though).

```{r}
fit_rgcca_tcga <- rgcca(blocks_tcga, response = 4, ncomp = c(5, 5, 5, 1))
plot(fit_rgcca_tcga, type = "ave")
plot(fit_rgcca_tcga, type = "samples", block = 1)
plot(fit_rgcca_tcga, type = "samples", block = 1, comp = 3:4)
plot(fit_rgcca_tcga, type = "samples", block = 1, comp = 4:5)
plot(fit_rgcca_tcga, type = "samples", block = 2)
plot(fit_rgcca_tcga, type = "samples", block = 2, comp = 3:4)
plot(fit_rgcca_tcga, type = "samples", block = 3)
plot(fit_rgcca_tcga, type = "samples", block = 3, comp = 3:4)
```

## Why preprocessing matters so much here

For omics data, preprocessing is often as important as the model itself. In practice, you may need:

- feature filtering to remove near-constant variables,
- log transformation or variance stabilization,
- sample matching across all blocks,
- scaling so that one block does not dominate only because of its measurement scale,
- and possibly batch correction before integration.

If these steps are ignored, the extracted components can mostly reflect technical structure rather than biology.

# From RGCCA to SGCCA

When the number of variables is (very) large, a dense solution can be difficult to interpret. This is where sparse RGCCA, usually referred to as SGCCA, becomes very useful.

The goal is not only to build integrated components, but also to retain a reduced subset of variables that contribute the most.

In SGCCA, there is a "sparsity" parameter that is directly linked to the number of zeros that we want to impose on the weights: a sparsity parameter close to 0 means that we want a lot of the weights to be null, a sparsity parameter close to 1 means that we do not want the weights to be forced towards zero.

## Tuning sparsity by cross-validation

We use cross-validation to estimate the optimal value for the sparsity parameter.
It consist in separating the data randomly into training, testing and validation sets
to test fairly different sparsity values.

![](img/diagram.png)

The result is shown below: plotting the result of `rgcca_cv` 

```{r sgcca cv, cache = TRUE}
set.seed(84251)

res.rgcca.cv <- rgcca_cv(
  blocks = blocks_tcga,
  ncomp = 3,
  method = "sgcca",
  prediction_model = "lda",
  metric = "Accuracy",
  response = 4,
  par_type = "sparsity",
  par_value = cbind(
    mRNA = seq(0.1, 1, by = 0.1),
    miRNA = seq(0.1, 1, by = 0.1),
    PROT = seq(0.1, 1, by = 0.1),
    Y = 1
  )
)

plot(res.rgcca.cv)
```

## Fit a sparse model using the selected values

```{r}
fit_sgcca_tcga_best <- rgcca(
  blocks = blocks_tcga,
  response = 4,
  ncomp = c(3, 3, 3, 1),
  method = "sgcca",
  sparsity = c(0.4, 0.4, 0.4, 1)
)

plot(fit_sgcca_tcga_best, type = "samples", block = 1)
plot(fit_sgcca_tcga_best, type = "samples", block = 2)
plot(fit_sgcca_tcga_best, type = "samples", block = 3)

plot(fit_sgcca_tcga_best, type = "samples", block = c(1, 2))
plot(fit_sgcca_tcga_best, type = "samples", block = c(1, 3))

plot(fit_sgcca_tcga_best, type = "ave")
```

## How to interpret sparse results

A sparse solution is especially useful when you want:

- shorter lists of candidate biomarkers,
- clearer biological interpretation,
- more focused downstream enrichment analysis,
- or a predictive model that is easier to communicate.

But sparsity also introduces choices and trade-offs:

- too much sparsity can miss relevant variables,
- too little sparsity can reduce interpretability,
- and the best tuning often depends on the downstream goal.

# Downstream interpretation with pathway analysis

A common next step is to use the weights from one block to rank variables and perform enrichment analysis.

```{r fgsea}
annot_table <- read.csv("../data/TCGA/annotation_table_mrna.csv", row.names = 1)
kegg_pathways <- readRDS("../data/TCGA/kegg_pathways.rds")

a1 <- fit_sgcca_tcga_best$a$mRNA[, 1]
all(annot_table$symbols_new == gsub("^mRNA_", "", names(a1)))
names(a1) <- annot_table$entrez


kegg_hsa_name <- kegg_pathways %>%
  dplyr::group_by(pathway_id, pathway_name) %>%
  summarize(n = length(entrez_gene))

# sum(names(a1) %in% kegg_pathways$entrez_gene)

fgsea_res <- fgsea::fgseaMultilevel(
  pathways = split(kegg_pathways$entrez_gene, kegg_pathways$pathway_id),
  stats = a1^2,
  minSize = 15,
  maxSize = 500,
  eps = 0
)

fgsea_res <- fgsea_res |>
  left_join(kegg_hsa_name, by = c("pathway" = "pathway_id")) %>%
  relocate(pathway_name, .after = 1)

fgsea_res %>% 
  filter(padj <= 0.05) %>%
  ggplot(aes(x = padj, y = reorder(pathway_name, padj))) +
  scale_y_discrete(labels = scales::label_wrap(20)) +
  geom_point(aes(size = n)) + 
  labs(x = "Adjusted p-value", y = "")
```

RGCCA and SGCCA help identify multiblock structure, but biological interpretation usually requires one more layer. Pathway analysis, gene-set analysis, and annotation-driven summaries help connect components back to mechanisms.

# Common pitfalls and good habits

## Typical pitfalls

- Treating all blocks as equally informative without checking quality.
- Connecting all blocks by default without biological justification.
- Using too many components and over-interpreting late components.
- Ignoring scaling and preprocessing differences across blocks.
- Interpreting weights without checking stability.
- Confusing a good visual separation with robust prediction.

## Good habits

- Start with block-wise exploratory plots.
- Make the design matrix explicit and justify it.
- Compare a dense RGCCA solution with a sparse SGCCA solution.
- Inspect both sample space and variable space.
- Use cross-validation when tuning sparsity for prediction-related goals.
- Connect selected variables to pathway or annotation analyses.

# Suggested workflow for beginners

1. Define the samples and blocks.
2. Check sample alignment and data quality.
3. Explore each block with simple graphics.
4. Decide whether the question is exploratory, supervised, or predictive.
5. Write down a first design matrix based on biology.
6. Run a simple RGCCA model.
7. Inspect sample plots, loadings, and AVE summaries.
8. Adjust the design or number of components if needed.
9. Move to SGCCA if interpretability or variable selection becomes important.
10. Use downstream interpretation tools on selected features.


