---
title: "Exposome Data Integration with Omic Data"
author: "Carles Hernandez-Ferer and Juan R. Gonzalez"
date: "`r doc_date()`"
package: "`r pkg_ver('omicRexposome')`"
abstract: >
  This is an introductory guide to integration analysis between exposome and omics data with R package omicRexposome. The document illustrates two types of analysis: 1) Association analysis, that are performed between exposome and a single omic data-set; and 2) Integration analysis where multiple data-sets, including exposome data, are analysed at the same time.
vignette: >
  %\VignetteIndexEntry{Exposome Data Integration with Omic Data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output: 
  BiocStyle::html_document
---

```{r setup, include=FALSE}
BiocStyle::markdown()
knitr::opts_chunk$set(echo = TRUE, warnings=FALSE)
```

# Introduction

`omicRexposome` is an R package designed to work join with `rexposome`. The aim of `omicRexposome` is to perform analysis joining exposome and omic data with the goal to find the relationship between a single or set of exposures (external exposome) and the behavior of a gene, a group of CpGs, the level of a protein, etc. Also to provide a series of tools to analyse exposome and omic data using standard methods from Biocondcutor.

## Installation

`omicRexposome` is currently in development and not available from CRAN nor Bioconductor. Anyway, the package can be installed by using devtools R package and taking the source from **Bioinformatic Research Group in Epidemiology**'s GitHub repository.

This can be done by opening an R session and typing the following code:

```{r eval=FALSE}
devtools::install_github("isglobal-brge/omicRexposome")
```

User must take into account that this sentence do not install the packages' dependencies.

## Pipeline

Two different types of analyses can be done with `omicRexposome`:

| Analysis          | `omicRexposome` function |
|:------------------|:-------------------------|
| Association Study | `association`            |
| Integration Study | `crossomics`             |

Both association and integration studies are based in objects of class `MultiDataSet`. A `MultiDataSet` object is a contained for multiple layers of sample information. Once the exposome data and the omics data are encapsulated in a `MultiDataSet` the object can be used for both association and integration studies.

The method `association` requires a `MultiDataSet` object having to types of information: the exposome data from an `ExposomeSet` object and omic information from objects of class `ExpressionSet`, `MethylationSet`, `SummarizedExperiment` or others. `ExposomeSet` objects are created with functions `read_exposome` and `load_exposome` from `rexposome` R package (see next section *Loading Exposome Data*) and encapsulates exposome data. The method `crossomics` expects a `MultiDataSet` with any number of different data-sets (at last two). Compared with `association` method, `crossomics` do not requires an `ExposomeSet`.

## Exposome and Omic Data

In order to illustrate the capabilities of `omicRexposome` and the exposome-omic analysis pipeline, we will use the data from `BRGdata` package. This package includes different omic-sets including methylation, transcriptome and proteome data-sets and an exposome-data set.

# Analysis

`omicRexposome` and `MultiDataSet` R packages are loaded using the standard library command:

```{r load_omicRexposome, message=FALSE}
library(omicRexposome)
library(MultiDataSet)
```

## Association Studies

The association studies are performed using the method `association`. This method requires, at last four, augments:

 1. Argument `object` should be filled with a `MultiDataSet` object.
 2. Argument `formula` should be filled with an expression containing the covariates used to adjust the model.
 3. Argument `expset` should be filled with the name that the exposome-set receives in the `MultiDataSet` object.
 4. Argument `omicset` should be filled with the name that the omic-set receives in the `MultiDataSet` object.
 
The argument `formula` should follow the pattern: `~sex+age`. The method `association` will fill the formula placing the exposures in the `ExposomeSet`m between `~` and the covariates `sex+age`.

`association` implements the `limma` pipeline using `lmFit` and `eBayes` in the extraction methods from `MultiDataSet`. The method takes care of the missing data in exposures, outcomes and omics data and locating and is subsets both data-sets, exposome data and omic data, by common samples. The argument `method` allows to select the fitting method used in `lmFit`. By default it takes the value `"ls"`  for *least squares* but it can also takes `"robust"` for *robust regression*.

The following subsections illustrates the usage of `association` with different types of omics data: *methylome*, *transcriptome* and *proteome*.

### Exposome - Transcriptome Data Association

First we get the exposome data from `BRGdata` package that we will use in the whole section.

```{r expos_data, cache=FALSE}
data("brge_expo", package = "brgedata")
class(brge_expo)
```


The aim of this analysis is to perform an association test between the gene expression levels and the exposures. So the first point is to obtain the transcriptome data from the `brgedata` package.

```{r gexp_data, cache=FALSE}
data("brge_gexp", package = "brgedata")
```

The association studies between exposures and transcriptome are done in the same way that the ones with methylome. The method used is `association`, that takes as input an object of `MultiDataSet` class with both exposome and expression data.

```{r gexp_analysis, message=FALSE, warning=FALSE}
mds <- createMultiDataSet()
mds <- add_genexp(mds, brge_gexp)
mds <- add_exp(mds, brge_expo)

gexp <- association(mds, formula=~Sex+Age, 
    expset = "exposures", omicset = "expression")
```

We can have a look to the number of hits and the lambda score of each analysis with the methods `tableHits` and `tableLambda`, seen in the previous section.

```{r gexp_tables}
hit <- tableHits(gexp, th=0.001)
lab <- tableLambda(gexp)
merge(hit, lab, by="exposure")
```

Since most of all models have a lambda under one, we should consider use *Surrogate Variable Analysis*. This can be done using the same `association` method but by setting the argument `sva` to `"fast"` so the pipeline of `isva` and `SmartSVA` R packages is applied. If `sva` is set to `"slow"` the applied. pipeline is the one from `sva` R package.

```{r gexp_analysis_sva, message=FALSE, warning=FALSE}
gexp <- association(mds, formula=~Sex+Age, 
    expset = "exposures", omicset = "expression", sva = "fast")
```

We can re-check the results creating the same table than before:

```{r gexp_tables_sva}
hit <- tableHits(gexp, th=0.001)
lab <- tableLambda(gexp)
merge(hit, lab, by="exposure")
```

The objects of class `ResultSet` have a method called `plotAssociation` that allows to create QQ Plots (that are another useful way to see if there are some inflation/deflation in the P-Values).

```{r gexp_plot_qq}
gridExtra::grid.arrange(
    plotAssociation(gexp, rid="Ben_p", type="qq") + 
        ggplot2::ggtitle("Transcriptome - Pb Association"),
    plotAssociation(gexp, rid="BPA_p", type="qq") + 
        ggplot2::ggtitle("Transcriptome - THM Association"),
    ncol=2
)
```

Following this line, the same method `plotAssociation` can be used to create volcano plots.

```{r gexp_plot_volcano}
gridExtra::grid.arrange(
    plotAssociation(gexp, rid="Ben_p", type="volcano", tPV=-log10(1e-04)) + 
        ggplot2::ggtitle("Transcriptome - Pb Association"),
    plotAssociation(gexp, rid="BPA_p", type="volcano", tPV=-log10(1e-04)) + 
        ggplot2::ggtitle("Transcriptome - THM Association"),
    ncol=2
)
```

### Exposome - Proteome Data Association

The proteome data-set included in `brgedata` has 47 proteins for 90 samples.

```{r prot_data, cache=FALSE}
data("brge_prot", package="brgedata")
brge_prot
```

The association analysis between exposures and proteome is also done using `association`.

```{r prot_analysis, message=FALSE, warning=FALSE}
mds <- createMultiDataSet()
mds <- add_eset(mds, brge_prot, dataset.type  ="proteome")
mds <- add_exp(mds, brge_expo)

prot <- association(mds, formula=~Sex+Age,
    expset = "exposures", omicset = "proteome")
```

The `tableHits` indicates that no association was found between the 47 proteins and the exposures.

```{r prot_hits}
tableHits(prot, th=0.001)
```

This is also seen in the Manhattan plot for proteins that can be obtained from `plotAssociation`.

```{r prot_plot_volcano}
gridExtra::grid.arrange(
    plotAssociation(prot, rid="Ben_p", type="protein") + 
        ggplot2::ggtitle("Proteome - Cd Association") +
        ggplot2::geom_hline(yintercept = 1, color = "LightPink"),
    plotAssociation(prot, rid="NO2_p", type="protein") + 
        ggplot2::ggtitle("Proteome - Cotinine Association") +
        ggplot2::geom_hline(yintercept = 1, color = "LightPink"),
    ncol=2
)
```

**NOTE**: A real Manhattan plot can be draw with `plot` method for `ResultSet` objects by setting the argument `type` to `"manhattan"`.


```{r, rm_asc, echo=FALSE, message=FALSE}
rm(prot, gexp)
gc()
```

## Integration Analysis

`omicRexposome` allows to study the relation between exposures and omic-features from another perspective, different from the association analyses. The integration analysis can be done, in `omicRexposome` using *multi canonical correlation analysis* or using *multiple co-inertia analysis*. The first methods is implemented in R package `PMA` (CRAN) and the second in `omicade4` R package (Bioconductor). The two methods are encapsulated in the `crossomics` method. 

The differences between `association` and `crossomics` are that the first method test association between two complete data-sets, by removing the samples having missing values in any of the involved data-sets, and the second try to find latent relationships between two or more sets.

Hence, we need to explore the missing data in the exposome data-set. This can be done using the methods `plotMissings` and `tableMissings` from `rexposome` R package.

```{r missing_exp, message=FALSE}
library(rexposome)
plotMissings(brge_expo, set = "exposures")
```

From the plot we can see that more of the exposures have up to 25% of missing values. Hence the first step in the integration analysis is to avoid missing values. so, we perform a fast imputation on the exposures side:

```{r impute_exp, message=FALSE}
brge_expo <- imputation(brge_expo)
```

`crossomics` function expects to obtain the different data-sets in a single labelled-list, in the argument called `list`. The argument `method` from `crossomics` function can be set to `mcia` (for *multiple co-inertia analysis*) or to `mcca` (for *multi canonical correlation analysis*).

The following code shows how to perform the integration of the exposome and the proteome. The method `crossomics` request a `MultiDataSet` object as input, containing the data-set to be integrated.

```{r crossomics_mcia, message=FALSE, warning=FALSE}
mds <- createMultiDataSet()
mds <- add_genexp(mds, brge_gexp)
mds <- add_eset(mds, brge_prot, dataset.type = "proteome")
mds <- add_exp(mds, brge_expo)

cr_mcia <- crossomics(mds, method = "mcia", verbose = TRUE)
cr_mcia
```

As can be seen, `crossomics` returns an object of class `ResultSet`. In the integration process, the different data-sets are subset by common samples. This is done taking advantage of `MultiDataSet` capabilities. 

The same is done when method is set to `mcca`.

```{r crossomics_mcca, message=FALSE, warning=FALSE, results='hide'}
cr_mcca <- crossomics(mds, method = "mcca", permute=c(4, 2))
cr_mcca
```

We used an extra argument (`permute`) into the previous call to `crossomics` using *multi canonical correlation analysis*. This argument allows to set the internal argument corresponding to `permutations` and `iterations`, that are used to tune-up internal parameters.

When a `ResultSet` is generated using `crossomics` the methods `plotHits`, `plotLambda` and `plotAssociation` can **NOT** be used. But the `plotIntegration` will help us to understand what was done. This method allows to provide the colors to be used on the plots:

```{r integration_colors}
colors <- c("green", "blue", "red")
names(colors) <- names(mds)
```

The graphical representation of the results from a *multiple co-inertia analysis* is a composition of four different plots.

```{r plot_pe_mcia, messages=FALSE, warnings=FALSE}
plotIntegration(cr_mcia, colors=colors)
```

The first plot (first row, first column) is the samples space. It illustrates how the different data-sets are related in terms of intra-sample variability (each data-set has a different color). The second plot (first row, second column) shows the feature space. The features of each set are drawn on the same components so the relation between each data-set can be seen (the features are colored depending of the set were they belong).

The third plot (second row, first column) shows the inertia of each component. The two first plots are drawn on the first and second component. Finally, the fourth plot shows the behavior of the data-sets.


A radar plots is obtained when `plotIntegration` is used on a `ResultSet` created though *multi canonical correlation analysis*.

```{r plot_pe_mcca, messages=FALSE, warnings=FALSE, fig.height=8, fig.width=9}
plotIntegration(cr_mcca, colors=colors)
```

This plot shows the features of the three data-sets in the same 2D space.The relation between the features can be understood by proximity. This means that the features that clusters, or that are in the same quadrant are related and goes in a different direction than the features in the other quadrants.

```{r, rm_crs}
rm(cr_mcia, cr_mcca)
```

# Session info {.unnumbered}

```{r sessionInfo, echo=FALSE}
sessionInfo()
```