interpretation.qmd

# Interpreting the integration results {#sec-interpretation}

```{r}
#| child: "_setup.qmd"
```

```{r loading-packages}
#| include: false

library(targets)
library(moiraine)
library(purrr)
library(dplyr)
library(ggplot2)
library(circlize)
library(patchwork)

## To avoid message when loading GGally
library(GGally)
```

```{r setup-visible}
#| eval: false

library(targets)
library(moiraine)

## For working with lists
library(purrr)

## For data-frames manipulation
library(dplyr)

## For colour palettes
library(ggplot2)
library(circlize)

## For manipulating patchworks of plots
library(patchwork)
```

In the previous chapters, we have applied different integration tools to our example multi-omics dataset. In each of these chapters, we have shown some method-specific visualisations that could be used to interpret the integration results. With the `moiraine` package, it is possible to use method-agnostic functions to inspect and interpret these results. These functionalities make use of the metadata in the `MultiDataSet` object to enrich plots and other summary outputs with information about the samples and the omics features, and allow users to obtain consistent visualisations across the integration methods used. 

<details>

<summary>`_targets.R` script</summary>

```{r previous-vignettes}
#| echo: false
#| results: asis

knitr::current_input() |> 
  get_previous_content() |> 
  cat()
```

</details>


## Dimension reduction methods {#sec-interpretation-dim-reduction}

Despite relying on very different statistical approaches, the different integration methods included in the pipeline all perform dimension reduction of the omics datasets through feature extraction. That is, they construct a small number of latent components/variables/dimensions (that we refer to as **latent dimensions** in the `moiraine` package) that capture as much information from the original datasets as possible. A dimension reduction approach typically returns, for each latent dimension constructed, two sets of values:

* **Features weight**: the contribution of the features from the different omics datasets to the latent dimension. All methods included in the pipeline construct latent dimensions as linear combinations of the original features, and therefore the features contribution is quantified by their weight in the linear combination for the corresponding latent dimension.

* **Samples score**: the projection of the samples onto the latent dimension.

In addition, the fraction or percentage of variance that each latent dimension explains in the different omics datasets is usually calculated.

Interpreting the results of a dimension reduction method involves:

* Understanding the source of the variation captured by each latent dimension: is a given latent dimension representing an important source of biological variation, such as effect of a treatment, or age of the samples? Or do they show a source of technical variation, for example highlighting a group of outlier samples with different omics profiles from the rest of the observations? Answering these questions allows us to identify which latent dimensions capture the biological phenomenon investigated, or whether there are some sources of noise that should be accounted for in follow-up experiments.

* Investigating which omics features are driving the latent dimensions: once we have identified some latent dimensions of interest, we can look at the features that contribute the most to these dimensions, in order to understand the molecular mechanisms or pathways involved. This is typically done after looking into the phenomenon captured by the latent dimensions, but can also help to identify it.

## Generating a standardised output {#sec-interpretation-standard-output}

### `get_output` function


In the `moiraine` package, the output of the different integration methods can be converted to a standardised output containing the pieces of information (features weight, samples score and percentage of variance explained) characteristic of dimension reduction methods, stored in a consistent format. This enables us to use functions for visualisation or analysis that can be applied to the results of any integration method, rather than having to implement one for each object type returned by the different integration packages.

The `get_output()` function transforms the output from any integration package included in `moiraine` into an `output_dimension_reduction` object, which is a list with three tibbles: `features_weight`, `samples_score` and `variance_explained`:

::: {.panel-tabset group="method"}
#### sPLS

::: {.targets-chunk}
```{targets get-output-spls}
tar_target(
  spls_output,
  get_output(spls_final_run)
)
```
:::

```{r show-output-spls}
tar_load(spls_output)
spls_output
```

#### sO2PLS

::: {.targets-chunk}
```{targets get-output-so2pls}
tar_target(
  so2pls_output,
  get_output(so2pls_final_run)
)
```
:::

```{r show-output-so2pls}
tar_load(so2pls_output)
so2pls_output
```

#### MOFA

::: {.targets-chunk}
```{targets get-output-mofa}
tar_target(
  mofa_output,
  get_output(mofa_trained)
)
```
:::

```{r show-output-mofa}
tar_load(mofa_output)
mofa_output
```

#### DIABLO

::: {.targets-chunk}
```{targets get-output-diablo}
tar_target(
  diablo_output,
  get_output(diablo_final_run)
)
```
:::

```{r show-output-diablo}
tar_load(diablo_output)
diablo_output
```
:::

The `features_weight` tibble contains one row per combination of feature and latent dimension. The ID of the features and the name of the dataset from which they originate are stored in the `feature_id` and `dataset` columns, respectively. The `latent_dimension` column gives the name of the latent dimension; this is a factor column. For each feature and latent dimension, the `weight` column shows the weight that was attributed to the feature for the corresponding latent dimension. In addition, the `importance` column contains the features **importance score**, which is computed as the absolute value of the features weight, divided by the maximum absolute weight across all features from the same omics dataset for the corresponding latent dimension. This importance score allows us to compare the contribution of the features across latent dimensions or integration methods, as the weight can be on different scales and thus cannot be directly compared. The importance scores range from 0 to 1. For any method performing feature selection (e.g. sPLS or DIABLO), features that were not selected for a given latent dimension are assigned a weight and importance score of 0.

The `samples score` tibble contains for each sample (`sample_id`) and latent dimension (`latent_dimension`) the sample's coordinate for the corresponding latent dimension. 

The `variance_explained` tibble gives for each latent dimension (`latent_dimension`) the proportion of variance explained (`prop_var_expl`) for each dataset (`dataset`). The values in `prop_var_expl` are between 0 and 1.

In addition, the name of the integration method used to obtain these results is stored in the `method` attribute of the `output_dimension_reduction` object:

::: {.panel-tabset group="method"}
#### sPLS

```{r show-method-attr-spls}
attr(spls_output, "method")
```


#### sO2PLS

```{r show-method-attr-so2pls}
attr(so2pls_output, "method")
```

#### MOFA

```{r show-method-attr-mofa}
attr(mofa_output, "method")
```

#### DIABLO

```{r show-method-attr-diablo}
attr(diablo_output, "method")
```
:::


For convenience, the `get_latent_dimensions()` function can be used on an `output_dimension_reduction` object to see the names of the latent dimensions (the levels used for the `latent_dimension` column in each tibble):

::: {.panel-tabset group="method"}
#### sPLS

```{r show-latent-dimensions-spls}
get_latent_dimensions(spls_output)
```


#### sO2PLS

```{r show-latent-dimensions-so2pls}
get_latent_dimensions(so2pls_output)
```

#### MOFA

```{r show-latent-dimensions-mofa}
get_latent_dimensions(mofa_output)
```

#### DIABLO

```{r show-latent-dimensions-diablo}
get_latent_dimensions(diablo_output)
```
:::

::: {.callout-note}
## Other methods covered by `get_output`

Note that both PCA and sPLS-DA (the method used for supervised features preselection in @sec-prefiltering-supervised) are also dimension reduction methods. Therefore, the `get_output()` function also converts `pcaRes` objects (from `run_pca()` or `pcaMethods::pca()`) and `mixo_splsda` objects (from `run_splsda()` or `mixOmics::splsda()`) into `output_dimension_reduction` objects.
:::

### Averaging latent dimensions over datasets

While MOFA computes one score per sample for each latent dimension created, sPLS, sO2PLS and DIABLO all compute one score per dataset for each sample and latent dimension. For each latent dimension, the samples score obtained for the different datasets are then compared, to assess the agreement or covariation between datasets. Ideally, these scores should be highly correlated across datasets, since the methods aim at maximising the variation between datasets, but it is not always the case. However, when they are highly correlated, it becomes redundant to interpret each dataset-specific version of the latent dimensions.

Instead, the `mixOmics` authors proposed a solution for `DIABLO`, which is to construct a weighted average space in which to represent the samples: for each latent component, the samples score are averaged over the different datasets. The weight attributed to each dataset is determined by how well the corresponding dataset discriminates between the samples group of interest. This way, rather than looking at samples score for each dataset for any given latent component, we can look at an average of them.

The `get_output()` function uses this idea to construct, for the output of sPLS, sO2PLS and DIABLO a set of average samples score for each latent dimension, rather than returning a set of samples score per dataset. For DIABLO, the average is weighted as explained above, while for sPLS and sO2PLS each dataset is given equal weight in the average. This calculation can be disabled in the `get_output()` function to extract the dataset-specific samples score, by setting the `use_average_dimensions` parameter to `FALSE`. Note that this only affects the `samples_score` tibble in terms of dimensions, but the name of the latent dimensions will change to reflect the dataset to which they refer.

::: {.panel-tabset group="method"}
#### sPLS

::: {.targets-chunk}
```{targets get-output-spls-no-average}
tar_target(
  spls_output_no_average,
  get_output(spls_final_run, use_average_dimensions = FALSE)
)
```
:::

```{r show-output-spls-no-average}
tar_load(spls_output_no_average)

get_latent_dimensions(spls_output_no_average)

nrow(spls_output$samples_score)
nrow(spls_output_no_average$samples_score)
```

#### sO2PLS

::: {.targets-chunk}
```{targets get-output-so2pls-no-average}
tar_target(
  so2pls_output_no_average,
  get_output(so2pls_final_run, use_average_dimensions = FALSE)
)
```
:::

```{r show-output-so2pls-no-average}
tar_load(so2pls_output_no_average)

get_latent_dimensions(so2pls_output_no_average)

nrow(so2pls_output$samples_score)
nrow(so2pls_output_no_average$samples_score)
```

#### DIABLO

::: {.targets-chunk}
```{targets get-output-diablo-no-average}
tar_target(
  diablo_output_no_average,
  get_output(diablo_final_run, use_average_dimensions = FALSE)
)
```
:::

```{r show-output-diablo-no-average}
tar_load(diablo_output_no_average)
  
get_latent_dimensions(diablo_output_no_average)

nrow(diablo_output$samples_score)
nrow(diablo_output_no_average$samples_score)
```
:::

## Percentage of variance explained

The first step in interpreting the results of a dimension reduction method is to assess how much variance is explained by the different latent dimensions for each dataset. The function `plot_variance_explained()` takes as input a `output_dimension_reduction` object and constructs a dataset-specific screeplot:

::: {.panel-tabset group="method"}
#### sPLS

```{r plot-variance-explained-spls}
plot_variance_explained(spls_output)
```


#### sO2PLS

```{r plot-variance-explained-so2pls}
#| fig.height: 9

plot_variance_explained(so2pls_output, ncol = 1) +
  theme(axis.text.x = element_text(size = 9, angle = 30, hjust = 1))
```

#### MOFA

```{r plot-variance-explained-mofa}
#| fig.height: 9

plot_variance_explained(mofa_output, ncol = 1) +
  theme(axis.text.x = element_text(size = 9, angle = 30, hjust = 1))
```

#### DIABLO

```{r plot-variance-explained-diablo}
#| fig.height: 8

plot_variance_explained(diablo_output, ncol = 2)
```
:::

## Samples plot {#sec-interpretation-samples-scores}


To understand the phenomenon driving the variation represented by the different latent components, we can investigate the coordinates of the samples in the space spanned by these latent components, in combination with information that we have about the samples.

We will start by loading the `MultiDataSet` object and generating some custom colour palettes for different samples covariates of interest, that we will use in the rest of the section.


```{r load-mo-set-de}
tar_load(mo_set_complete)

## Choosing colour palettes for the samples covariates
palette_status <- scale_colour_brewer(palette = "Set1")
palette_feedlot <- scale_colour_brewer(palette = "Set2")
palette_geno_comp <- scale_colour_brewer(palette = "Dark2")
palette_rnaseq_batch <- scale_colour_brewer(palette = "Paired")

## For some plots will need both colour and fill
palette_status_fill <- scale_colour_brewer(
  palette = "Set1", 
  aesthetics = c("colour", "fill")
)
```

### Matrix of scatter plots {#sec-interpretation-samples-ggpairs}

The `plot_samples_score()` function shows the samples score in a matrix of scatter plots that represent all possible two-by-two combinations of the latent dimensions. The plots in the resulting matrix are redundant: the lower plots (below the diagonal) are just a rotated version of the upper plots (above the diagonal). The function takes as input an `output_dimension_reduction` object. In addition, it accepts a `MultiDataSet` object to extract information about the samples that can then be used to customise the plot. Since the plots are redundant, different properties of the samples can be shown in the upper and lower plots. For example, we will colour the samples according to disease status in the upper plots, while representing some other covariate in the lower plots. We will use the shape of the points to illustrate gender in both upper and lower plots. By default, all latent dimensions are included in the plots, but it is possible to select some of them by passing their names to the `latent_dimensions` argument. There are more options to further customise the plot, which are shown in the function's help.


::: {.panel-tabset group="method"}
#### sPLS

```{r spls-plot-samples-scores}
#| fig.height: 9

plot_samples_score(
  spls_output,
  mo_data = mo_set_complete,
  colour_upper = "status",
  scale_colour_upper = palette_status,
  shape_upper = "gender",
  colour_lower = "feedlot",
  scale_colour_lower = palette_feedlot
) +
  theme(legend.box = "vertical")
```

Component 1 clearly splits the control and BRD animals, while component 2 seems to separate two samples from the rest. There is no clear separation of the samples by feedlot or gender.

#### sO2PLS

Here we will focus on the joint components. Since there is only one, instead of a matrix of scatterplot, the function will display the samples score with a violin plot, and the properties selected for the upper plots will be used to customise the points:

```{r so2pls-plot-samples-scores}
#| fig.height: 6

plot_samples_score(
  so2pls_output,
  latent_dimensions = "joint component 1",
  mo_data = mo_set_complete,
  colour_upper = "status",
  scale_colour_upper = palette_status,
  shape_upper = "gender",
  colour_lower = "feedlot",
  scale_colour_lower = palette_feedlot
)
```

Joint component 1 separates the BRD and control animals, although some of them overlap.

#### MOFA

For clarity we will focus here on the first four factors:

```{r mofa-plot-samples-scores}
#| fig.height: 9

plot_samples_score(
  mofa_output,
  latent_dimensions = paste("Factor", 1:4),
  mo_data = mo_set_complete,
  colour_upper = "status",
  scale_colour_upper = palette_status,
  shape_upper = "gender",
  colour_lower = "geno_comp_cluster",
  scale_colour_lower = palette_geno_comp
) +
  theme(legend.box = "vertical")
```

Factor 1 clearly separates the control and BRD animals. Factor 2 separates samples from the genomics composition cluster K1 from the other two clusters. None of the factors are separating the animals by gender, indicating that this covariate does not have a strong effect in the omics measurements.

#### DIABLO

```{r diablo-plot-samples-scores}
#| fig.height: 9

plot_samples_score(
  diablo_output,
  mo_data = mo_set_complete,
  colour_upper = "status",
  scale_colour_upper = palette_status,
  shape_upper = "gender",
  colour_lower = "rnaseq_batch",
  scale_colour_lower = palette_rnaseq_batch
) +
  theme(legend.box = "vertical")
```

As DIABLO is a supervised method, it is no surprise that the first component clearly separates the two disease status groups. The second latent component separates samples according to the RNAseq batch. From this plot it is not clear what components 3 and 4 represent.
:::

### Scatter plot for pair of latent dimensions {#sec-interpretation-samples-pair}

We can also focus on a specific pair of latent dimensions and represent the samples in the space spanned by these two dimensions, with the `plot_samples_score_pair()` function. The two latent dimensions to plot are selected by passing their names to the function. This function also accepts a `MultiDataSet` object, whose samples metadata will be used to customise the plot.

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-plot-samples-score-pair}
plot_samples_score_pair(
  spls_output,
  paste("Component", 1:2),
  mo_data = mo_set_complete,
  colour_by = "status",
  shape_by = "gender"
) +
  palette_status
```

The two disease status groups are well separated in this space.

#### sO2PLS

To showcase this function, we will look at the first two metabolomics-specific latent components:

```{r so2pls-plot-samples-score-pair}
plot_samples_score_pair(
  so2pls_output,
  paste("metabolome specific component", 1:2),
  mo_data = mo_set_complete,
  colour_by = "status",
  shape_by = "gender"
) +
  palette_status
```

These two specific components are not related to disease status or gender.

#### MOFA

```{r mofa-plot-samples-score-pair}
plot_samples_score_pair(
  mofa_output,
  paste("Factor", 1:2),
  mo_data = mo_set_complete,
  colour_by = "status",
  shape_by = "geno_comp_cluster"
) +
  palette_status
```

We see that these first two factors discriminate the disease status groups and the genomics composition clusters.

#### DIABLO

```{r diablo-plot-samples-score-pair}
plot_samples_score_pair(
  diablo_output,
  paste("Component", 1:2),
  mo_data = mo_set_complete,
  colour_by = "status",
  shape_by = "rnaseq_batch"
) +
  palette_status
```

We see that these first two components discriminate the disease status groups and the RNAseq batches.
:::

### Samples score vs covariate {#sec-interpretation-samples-covariate}

Lastly, we can plot the samples score for each latent dimension against a samples covariate of interest, with the `plot_samples_score_covariate()` function. This is useful to focus on one dimension at a time rather than pairs of them. By default, the function displays all latent dimensions, but we can focus on a subset of them through the `latent_dimensions` arguments. The function needs a `MultiDataSet` object to extract information about the samples. The plot generated will depend on whether the covariate of interest is categorical or continuous. We will show both options below. In addition to the main covariate of interest that will be displayed on the x-axis, the colour and shape of the points can be further customised according to other samples properties.

First, with a categorical covariate:

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-plot-samples-score-covariate-cat}
plot_samples_score_covariate(
  spls_output,
  mo_set_complete,
  "status",
  colour_by = "status"
) +
  palette_status_fill
```


#### sO2PLS

For clarity we will focus on the first joint component only.

```{r so2pls-plot-samples-score-covariate-cat}
plot_samples_score_covariate(
  so2pls_output,
  mo_set_complete,
  "status",
  colour_by = "status",
  latent_dimensions = "joint component 1"
) +
  palette_status_fill
```

#### MOFA

For clarity we will focus on the first two factors only.

```{r mofa-plot-samples-score-covariate-cat}
plot_samples_score_covariate(
  mofa_output,
  mo_set_complete,
  "status",
  colour_by = "status",
  shape_by = "geno_comp_cluster",
  latent_dimensions = paste("Factor", 1:2)
) +
  palette_status_fill
```

#### DIABLO

For clarity we will focus on the first two components only.

```{r diablo-plot-samples-score-covariate-cat}
plot_samples_score_covariate(
  diablo_output,
  mo_set_complete,
  "status",
  colour_by = "status",
  shape_by = "rnaseq_batch",
  latent_dimensions = paste("Component", 1:2)
) +
  palette_status_fill
```
:::


Now, with a continuous covariate:

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-plot-samples-score-covariate-cont}
plot_samples_score_covariate(
  spls_output,
  mo_set_complete,
  "day_on_feed",
  colour_by = "status"
) +
  palette_status_fill
```


#### sO2PLS

For clarity we will focus on the first joint component only.

```{r so2pls-plot-samples-score-covariate-cont}
plot_samples_score_covariate(
  so2pls_output,
  mo_set_complete,
  "day_on_feed",
  colour_by = "status",
  latent_dimensions = "joint component 1"
) +
  palette_status_fill
```

#### MOFA

For clarity we will focus on the first two factors only.

```{r mofa-plot-samples-score-covariate-cont}
plot_samples_score_covariate(
  mofa_output,
  mo_set_complete,
  "day_on_feed",
  colour_by = "status",
  shape_by = "geno_comp_cluster",
  latent_dimensions = paste("Factor", 1:2)
) +
  palette_status_fill
```

#### DIABLO

For clarity we will focus on the first two components only.

```{r diablo-plot-samples-score-covariate-cont}
plot_samples_score_covariate(
  diablo_output,
  mo_set_complete,
  "day_on_feed",
  colour_by = "status",
  shape_by = "rnaseq_batch",
  latent_dimensions = paste("Component", 1:2)
) +
  palette_status_fill
```
:::

## Feature plots

Next, we can have a look at the features contributing to each latent dimension.

### Features weight distribution {#sec-interpretation-features-distr}

We can start by looking at the distribution of features weight for each latent dimension across the datasets, with the `plot_features_weight_distr()` function. Again, the function takes an `output_dimension_reduction` object as first input; by default it will show all latent dimensions and datasets, but they can be specified through the `latent_dimensions` and `datasets` arguments. The function shows by default the distribution of the "signed" importance scores, i.e. the features weights normalised by the highest absolute weight value for the corresponding latent dimension and dataset, but we can also show the distribution of importance scores (which are absolute values) and non-transformed features weight (through the `features_metric` argument). Note that the function returns a [patchwork](https://patchwork.data-imaginist.com/) of ggplots.

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-plot-features-weight-distr}
plot_features_weight_distr(spls_output)
```

As sPLS performs feature selection, most of the weights are equal to 0.

#### sO2PLS

```{r so2pls-plot-features-weight-distr}
#| fig.height: 15

plot_features_weight_distr(so2pls_output) +
  ## showing two plots per column (from patchwork package)
  plot_layout(ncol = 2)
```

As sO2PLS performs feature selection for the joint components, most of the weights are equal to 0 for joint component 1.

#### MOFA

We will only look at factors 1 to 4 for clarity:

```{r mofa-plot-features-weight-distr}
#| fig.height: 15

plot_features_weight_distr(
  mofa_output,
  latent_dimensions = paste("Factor", 1:4)
) +
  ## showing one plot per column (from patchwork package)
  plot_layout(ncol = 1)
```

We can see for example for factor 1 that only a couple of metabolites have an importance score above 0.5, this small number of features are driving the variation captured by factor 1.

#### DIABLO

```{r diablo-plot-features-weight-distr}
#| fig.height: 15

plot_features_weight_distr(
  diablo_output
) +
  ## showing one plot per column (from patchwork package)
  plot_layout(ncol = 1)
```

As DIABLO performs feature selection, most of the weights are equal to zero.

:::

### Comparing features importance between latent components {#sec-interpretation-features-pair}

We can also compare the importance score given to the features between any two latent dimensions with the `plot_features_weight_pair()` function.

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-features-weight-pair}
plot_features_weight_pair(
  spls_output,
  paste("Component", 1:2),
  mo_data = mo_set_complete,
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  )
)
```

All genes and most metabolites selected for the two components are different.

#### sO2PLS

```{r so2pls-features-weight-pair}
plot_features_weight_pair(
  so2pls_output,
  paste("metabolome specific component", 1:2),
  mo_data = mo_set_complete,
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  )
)
```

The plot shows an interesting pattern: the weights of the compound for these two specific components are negatively correlated, except for a couple of outliers.

#### MOFA

```{r mofa-features-weight-pair}
plot_features_weight_pair(
  mofa_output,
  paste("Factor", 1:2),
  mo_data = mo_set_complete,
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  )
)
```

We can see that the genomic markers are mostly given weights of opposite signs for the first two MOFA factors (e.g. positive weight for factor 1 and negative weight for factor 2). For the transcriptomics and metabolomics datasets, the importance score are mostly uncorrelated between the first two factors.

#### DIABLO

```{r diablo-features-weight-pair}
plot_features_weight_pair(
  diablo_output,
  paste("Component", 1:2),
  mo_data = mo_set_complete,
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  )
)
```

We can see that different features were selected for the first two latent components from each omics dataset.

:::

These plots allow us to quickly check whether different latent dimensions capture different or related aspects of the data. The top 5 features according to their consensus importance score (see @sec-comparison-consensus) are highlighted in red.

By default, the function shows the signed importance score of the features, i.e. their importance score to which the sign of their weight was added. This can be changed through the `features_metric` argument of the function.


### Plotting top features weight  {#sec-interpretation-features-top-plot}

We can then investigate which features are assigned the highest weights for each latent dimension. These are the features driving the variation captured by the latent components. This is done with the `plot_top_features()` function. 

By default, the feature IDs will be used as labels in the plot, but by passing a `MultiDataSet` object to the function we can use a column from the features metadata table instead. To do so, we need to pass a named list to the `label_cols` argument: the names of the elements in the list should correspond to dataset names in the `MultiDataSet` object, and the value should be the name of the column in the features metadata table of the corresponding dataset to use as feature label. The number of features displayed is controlled through the `top_n` argument (the default value is 20). Note that if a lower number of features is selected for a certain dataset and latent dimension, only the selected features will be displayed. Again, the function returns a [patchwork](https://patchwork.data-imaginist.com/) of ggplots.


::: {.panel-tabset group="method"}
#### sPLS

```{r spls-plot-top-features}
#| fig.height: 8

plot_top_features(
  spls_output,
  mo_data = mo_set_complete,
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  )
)
```

We can see that citric acid and D-mannose have an importance score close to 1 for component 1. Citric acid has a negative weight, meaning that its abundance is negatively correlated with component 1, while D-mannose has a positive weight, so its abundance is positively correlated with component 1. For the second component, only two genes have an importance score above 0.5.


#### sO2PLS

We will focus on joint component 1:

```{r so2pls-plot-top-features}
plot_top_features(
  so2pls_output,
  mo_data = mo_set_complete,
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  ),
  latent_dimensions = "joint component 1"
)
```

While all top 20 genes have an importance score above 0.5 for joint component 1, only citric acid is given a high importance score for joint component 1 in the metabolomics dataset.

#### MOFA

We will focus on the first two factors here:

```{r mofa-plot-top-features}
#| fig.width: 10
#| fig.height: 8

plot_top_features(
  mofa_output,
  mo_data = mo_set_complete,
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  ),
  latent_dimensions = paste("Factor", 1:2)
)
```

For factor 1, citric acid is the only metabolite with an importance score above 0.75, while in the genomics dataset four markers have a high importance score, and most of the top 20 genes have an importance score aroung 0.7 or higher.

#### DIABLO

For clarity, we will focus on the first two latent components:

```{r diablo-plot-top-features}
#| fig.height: 8

plot_top_features(
  diablo_output,
  mo_data = mo_set_complete,
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  ),
  latent_dimensions = paste("Component", 1:2)
)
```

For component 1, we can see three markers from the genomics dataset and three genes with a high importance score. From the metabolomics dataset, less than 20 features were retained.

:::

### Extracting top features {#sec-interpretation-features-top-table}

It is also possible to extract information about the top contributing features as a table with the function `get_top_features()`. The function offers two options:

-   Return for each dataset the top N features contributing to the latent dimensions of interest. N is specified through the `n_features` argument;

-   Return for each dataset features whose importance score is at least equal to a certain threshold. The threshold is specified through the `min_importance` argument.

The first option is preferred when many features contribute to a given factor. The second option has been implemented for cases where some latent dimensions are driven by a small number of features (so for example in the top 10 contributing features, only 3 features would actually have a high importance score). By default, the function returns the top contributing features from all datasets for all latent dimensions, but we can focus on some datasets and/or latent dimensions of interest by passing their names to the `datasets` and `latent_dimensions` arguments of the function. Additionally, we can pass to the function a `MultiDataSet` object (through the `mo_data` argument); the function will extract information about the features from the features metadata.

We illustrate the two possible options below. For the second option, we add features information to the table:

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-get-top-features}
get_top_features(spls_output, n_features = 3) |> 
  head()

get_top_features(
  spls_output,
  min_importance = 0.8,
  mo_data = mo_set_complete
) |> 
  head()
```

#### sO2PLS

```{r so2pls-get-top-features}
get_top_features(so2pls_output, n_features = 3) |> 
  head()

get_top_features(
  so2pls_output,
  min_importance = 0.8,
  mo_data = mo_set_complete
) |> 
  head()
```

#### MOFA

```{r mofa-get-top-features}
get_top_features(mofa_output, n_features = 3) |> 
  head()

get_top_features(
  mofa_output,
  min_importance = 0.8,
  mo_data = mo_set_complete
) |> 
  head()
```

#### DIABLO

```{r diablo-get-top-features}
get_top_features(diablo_output, n_features = 3) |> 
  head()

get_top_features(
  diablo_output,
  min_importance = 0.8,
  mo_data = mo_set_complete
) |> 
  head()
```
:::



### Extracting the selected features {#sec-interpretation-features-selected}

For methods that perform feature selection (i.e. sPLS, sO2PLS and DIABLO), we can extract information about the selected features with the `get_selected_features()` function. As with `get_top_features()`, if a `MultiDataSet` object is passed to the function, it will extract information about the features from the features metadata tables and return them. If applied to the results of a method that does not perform feature selection, all features will be returned.


::: {.panel-tabset group="method"}
#### sPLS

```{r spls-get-selected-features}
selected_features <- get_selected_features(spls_output) 

dim(selected_features)
head(selected_features)
```

#### sO2PLS

```{r so2pls-get-selected-features}
selected_features <- get_selected_features(so2pls_output) 

dim(selected_features)
head(selected_features)
```

#### MOFA

```{r mofa-get-selected-features}
selected_features <- get_selected_features(mofa_output) 

dim(selected_features)
head(selected_features)
```

#### DIABLO

```{r diablo-get-selected-features}
selected_features <- get_selected_features(diablo_output) 

dim(selected_features)
head(selected_features)
```
:::

### Features measurements {#sec-interpretation-features-meas}

Finally, we can investigate further some features of interest using the functions shown in @sec-inspecting-plot-data. For the example, we will extract the top 3 contributing features from each dataset for the first latent dimension generated with the methods:

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-pull-top-features}
spls_top_features <- get_top_features(
  spls_output,
  n_features = 3,
  latent_dimensions = "Component 1"
) |> 
  pull(feature_id)

spls_top_features
```


#### sO2PLS

```{r so2pls-pull-top-features}
so2pls_top_features <- get_top_features(
  so2pls_output,
  n_features = 3,
  latent_dimensions = "joint component 1"
) |> 
  pull(feature_id)

so2pls_top_features
```

#### MOFA

```{r mofa-pull-top-features}
mofa_top_features <- get_top_features(
  mofa_output,
  n_features = 3,
  latent_dimensions = "Factor 1"
) |> 
  pull(feature_id)

mofa_top_features
```

#### DIABLO

```{r diablo-pull-top-features}
diablo_top_features <- get_top_features(
  diablo_output,
  n_features = 3,
  latent_dimensions = "Component 1"
) |> 
  pull(feature_id)

diablo_top_features
```
:::

We can use the `plot_data_heatmap()` function to visualise the measurements for the features across the samples as a heatmap. By passing a `MultiDataSet` object to the function, we can add information about the samples and the features around the heatmap, and change the label used for the features.

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-plot-top-features-heatmap}
#| fig.width: 10
#| fig.height: 4

plot_data_heatmap(
  mo_set_complete,
  spls_top_features,
  center = TRUE,
  scale = TRUE,
  show_column_names = FALSE,
  only_common_samples = TRUE,
  samples_info = c("status", "day_on_feed"),
  features_info = c("de_status"),
  colours_list = list(
    "status" = c("Control" = "gold", "BRD" = "lightblue"),
    "day_on_feed" = colorRamp2(c(5, 70), c("white", "pink3")),
    "de_status" = c("downregulated" = "deepskyblue", 
                    "upregulated" = "chartreuse3")
  ),
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  ),
  truncate = 20
)
```


#### sO2PLS

```{r so2pls-plot-top-features-heatmap}
#| fig.width: 10
#| fig.height: 4

plot_data_heatmap(
  mo_set_complete,
  so2pls_top_features,
  center = TRUE,
  scale = TRUE,
  show_column_names = FALSE,
  only_common_samples = TRUE,
  samples_info = c("status", "day_on_feed"),
  features_info = c("de_status"),
  colours_list = list(
    "status" = c("Control" = "gold", "BRD" = "lightblue"),
    "day_on_feed" = colorRamp2(c(5, 70), c("white", "pink3")),
    "de_status" = c("downregulated" = "deepskyblue", 
                    "upregulated" = "chartreuse3")
  ),
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  ),
  truncate = 20
)
```

#### MOFA

```{r mofa-plot-top-features-heatmap}
#| fig.width: 10
#| fig.height: 6

plot_data_heatmap(
  mo_set_complete,
  mofa_top_features,
  center = TRUE,
  scale = TRUE,
  show_column_names = FALSE,
  only_common_samples = TRUE,
  samples_info = c("status", "day_on_feed"),
  features_info = c("de_status"),
  colours_list = list(
    "status" = c("Control" = "gold", "BRD" = "lightblue"),
    "day_on_feed" = colorRamp2(c(5, 70), c("white", "pink3")),
    "de_status" = c("downregulated" = "deepskyblue", 
                    "upregulated" = "chartreuse3")
  ),
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  ),
  truncate = 20
)
```

#### DIABLO

```{r diablo-plot-top-features-heatmap}
#| fig.width: 10
#| fig.height: 6

plot_data_heatmap(
  mo_set_complete,
  diablo_top_features,
  center = TRUE,
  scale = TRUE,
  show_column_names = FALSE,
  only_common_samples = TRUE,
  samples_info = c("status", "day_on_feed"),
  features_info = c("de_status"),
  colours_list = list(
    "status" = c("Control" = "gold", "BRD" = "lightblue"),
    "day_on_feed" = colorRamp2(c(5, 70), c("white", "pink3")),
    "de_status" = c("downregulated" = "deepskyblue", 
                    "upregulated" = "chartreuse3")
  ),
  label_cols = list(
    "rnaseq" = "Name",
    "metabolome" = "name"
  ),
  truncate = 20
)
```
:::

Alternatively, we can display the features’ measurements against some samples covariate, with the `plot_data_covariate()` function:

::: {.panel-tabset group="method"}
#### sPLS

```{r spls-plot-top-features-covariate}
plot_data_covariate(
  mo_set_complete,
  "status",
  spls_top_features,
  only_common_samples = TRUE,
  label_cols = list(
    "rnaseq" = "Name", 
    "metabolome" = "name"
  )
)
```

#### sO2PLS

```{r so2pls-plot-top-features-covariate}
plot_data_covariate(
  mo_set_complete,
  "status",
  so2pls_top_features,
  only_common_samples = TRUE,
  label_cols = list(
    "rnaseq" = "Name", 
    "metabolome" = "name"
  )
)
```

#### MOFA

```{r mofa-plot-top-features-covariate}
plot_data_covariate(
  mo_set_complete,
  "status",
  mofa_top_features,
  only_common_samples = TRUE,
  label_cols = list(
    "rnaseq" = "Name", 
    "metabolome" = "name"
  )
)
```

#### DIABLO

```{r diablo-plot-top-features-covariate}
plot_data_covariate(
  mo_set_complete,
  "status",
  diablo_top_features,
  only_common_samples = TRUE,
  label_cols = list(
    "rnaseq" = "Name", 
    "metabolome" = "name"
  )
)
```
:::

## Recap -- targets list

For convenience, here is the list of targets that we created in this section:

<details>

<summary>Targets list for results interpretation</summary>

In the targets script, before the targets list, we will add the following code to create the colour palettes used in some of the plots:

::: {.targets-chunk}
```{targets targets-preamble}
## In the preamble of the _targets.R script

## Choosing colour palettes for the samples covariates
palette_status <- scale_colour_brewer(palette = "Set1")
palette_feedlot <- scale_colour_brewer(palette = "Set2")
palette_geno_comp <- scale_colour_brewer(palette = "Dark2")
palette_rnaseq_batch <- scale_colour_brewer(palette = "Paired")
```
:::

::: {.panel-tabset group="method"}
#### sPLS

::: {.targets-chunk}
```{targets spls-recap-targets-list}
list(
  ## Generating standardised output
  tar_target(
    spls_output,
    get_output(spls_final_run)
  ),
  
  ## Percentage of variance explained
  tar_target(
    spls_plot_variance_explained,
    plot_variance_explained(spls_output)
  ),
  
  ## Samples score matrix plot
  tar_target(
    spls_samples_scores_plot,
    plot_samples_score(
      spls_output,
      mo_data = mo_set_complete,
      colour_upper = "status",
      scale_colour_upper = palette_status,
      shape_upper = "gender",
      colour_lower = "feedlot",
      scale_colour_lower = palette_feedlot
    ) +
      theme(legend.box = "vertical")
  ),
  
  ## Distribution of features weight
  tar_target(
    spls_features_weight_distribution,
    plot_features_weight_distr(spls_output)
  ),
  
  ## Plot of top contributing features
  tar_target(
    spls_top_features_plot,
    plot_top_features(
      spls_output,
      mo_data = mo_set_complete,
      label_cols = list(
        "rnaseq" = "Name",
        "metabolome" = "name"
      )
    )
  ),
  
  ## Table of selected features
  tar_target(
    spls_selected_features,
    get_selected_features(spls_output)
  )
)
```
:::

#### sO2PLS

::: {.targets-chunk}
```{targets so2pls-recap-targets-list}
list(
  ## Generating standardised output
  tar_target(
    so2pls_output,
    get_output(so2pls_final_run)
  ),
  
  ## Percentage of variance explained
  tar_target(
    so2pls_plot_variance_explained,
    plot_variance_explained(so2pls_output, ncol = 1) +
      theme(axis.text.x = element_text(size = 9, angle = 30, hjust = 1))
  ),
  
  ## Samples score matrix plot
  tar_target(
    so2pls_samples_scores_plot,
    plot_samples_score(
      so2pls_output,
      latent_dimensions = "joint component 1",
      mo_data = mo_set_complete,
      colour_upper = "status",
      scale_colour_upper = palette_status,
      shape_upper = "gender"
    )
  ),
  
  ## Distribution of features weight
  tar_target(
    so2pls_features_weight_distribution,
    plot_features_weight_distr(so2pls_output) +
      plot_layout(ncol = 2)
  ),
  
  ## Plot of top contributing features
  tar_target(
    so2pls_top_features_plot,
    plot_top_features(
      so2pls_output,
      mo_data = mo_set_complete,
      label_cols = list(
        "rnaseq" = "Name",
        "metabolome" = "name"
      )
    )
  ),
  
  ## Table of selected features
  tar_target(
    so2pls_selected_features,
    get_selected_features(
      so2pls_output,
      latent_dimensions = "joint component 1"
    )
  )
)
```
:::

#### MOFA

::: {.targets-chunk}
```{targets mofa-recap-targets-list}
list(
  ## Generating standardised output
  tar_target(
    mofa_output,
    get_output(mofa_trained)
  ),
  
  ## Percentage of variance explained
  tar_target(
    mofa_plot_variance_explained,
    plot_variance_explained(mofa_output, ncol = 1) +
      theme(axis.text.x = element_text(size = 9, angle = 30, hjust = 1))
  ),
  
  ## Samples score matrix plot
  tar_target(
    mofa_samples_scores_plot,
    plot_samples_score(
      mofa_output,
      latent_dimensions = paste("Factor", 1:4),
      mo_data = mo_set_complete,
      colour_upper = "status",
      scale_colour_upper = palette_status,
      shape_upper = "gender",
      colour_lower = "geno_comp_cluster",
      scale_colour_lower = palette_geno_comp
    ) +
      theme(legend.box = "vertical")
  ),
  
  ## Distribution of features weight
  tar_target(
    mofa_features_weight_distribution,
    plot_features_weight_distr(
      mofa_output,
      latent_dimensions = paste("Factor", 1:4)
    ) +
      plot_layout(ncol = 1)
  ),
  
  ## Plot of top contributing features
  tar_target(
    mofa_top_features_plot,
    plot_top_features(
      mofa_output,
      mo_data = mo_set_complete,
      label_cols = list(
        "rnaseq" = "Name",
        "metabolome" = "name"
      ),
      latent_dimensions = paste("Factor", 1:2)
    )
  ),
  
  ## Table of top contributing features
  tar_target(
    mofa_top_features,
    get_top_features(
      mofa_output,
      min_importance = 0.8,
      mo_data = mo_set_complete
    )
  )
)
```
:::

#### DIABLO

::: {.targets-chunk}
```{targets diablo-recap-targets-list}
list(
  ## Generating standardised output
  tar_target(
    diablo_output,
    get_output(diablo_final_run)
  ),
  
  ## Percentage of variance explained
  tar_target(
    diablo_plot_variance_explained,
    plot_variance_explained(diablo_output, ncol = 2)
  ),
  
  ## Samples score matrix plot
  tar_target(
    diablo_samples_scores_plot,
    plot_samples_score(
      diablo_output,
      mo_data = mo_set_complete,
      colour_upper = "status",
      scale_colour_upper = palette_status,
      shape_upper = "gender",
      colour_lower = "rnaseq_batch",
      scale_colour_lower = palette_rnaseq_batch
    ) +
      theme(legend.box = "vertical")
  ),
  
  ## Distribution of features weight
  tar_target(
    diablo_features_weight_distribution,
    plot_features_weight_distr(
      diablo_output
    ) +
      plot_layout(ncol = 1)
  ),
  
  ## Plot of top contributing features
  tar_target(
    diablo_top_features_plot,
    plot_top_features(
      diablo_output,
      mo_data = mo_set_complete,
      label_cols = list(
        "rnaseq" = "Name",
        "metabolome" = "name"
      ),
      latent_dimensions = paste("Component", 1:2)
    )
  ),
  
  ## Table of top contributing features
  tar_target(
    diablo_selected_features,
    get_selected_features(diablo_output) 
  )
)
```
:::

:::

</details>