MEAL 1.10.1

Experiments comprising different omic data analyses are becoming common. In this cases, the researchers are interested in analyzing the different omic types independently as well as finding association between the different layers.

In this vignette, we present how address to an analysis of methylation and expression data using *MEAL*. We will show how to analyze each data type independently and how to integrate the results of both analyses. We will also exemplify how to run an analysis when we are interested in our target region.

We will use the data from *brgedata* package. This package contains data from a spanish children cohort. We will work with the methylation and the expression data encapsulated in a `GenomicRatioSet`

and in an `ExpressionSet`

respectively. In our analysis, we will evaluate the effect of sex on methylation and expression. We will also deeply explore changes in a region of chromosome 11 around MMP3, a gene differently expressed in men and women (chr11:102600000-103300000).

Let’s start by loading the required packages and data:

```
library(MEAL)
library(brgedata)
library(MultiDataSet)
library(missMethyl)
library(minfi)
library(GenomicRanges)
library(ggplot2)
```

From the six objects of *brgedata*, we will focus on those containing methylation and expression data: brge_methy and brge_gexp.

brge_methy is a `GenomicRatioSet`

with methylation data corresponding to the Illumina Human Methylation 450K. It contains 392277 probes and 20 samples. The object contains 9 phenotypic variables (age, sex and cell types proportions):

```
data(brge_methy)
brge_methy
```

```
## class: GenomicRatioSet
## dim: 392277 20
## metadata(0):
## assays(1): Beta
## rownames(392277): cg13869341 cg24669183 ... cg26251715 cg25640065
## rowData names(14): Forward_Sequence SourceSeq ...
## Regulatory_Feature_Group DHS
## colnames(20): x0017 x0043 ... x0077 x0079
## colData names(9): age sex ... Mono Neu
## Annotation
## array: IlluminaHumanMethylation450k
## annotation: ilmn12.hg19
## Preprocessing
## Method: NA
## minfi version: NA
## Manifest version: NA
```

`colData(brge_methy)`

```
## DataFrame with 20 rows and 9 columns
## age sex NK Bcell
## <numeric> <factor> <numeric> <numeric>
## x0017 4 Female -5.81386004885394e-19 0.173507848676514
## x0043 4 Female 0.00164826459258012 0.182017163135389
## x0036 4 Male 0.0113186314286605 0.169017297114852
## x0041 4 Female 0.00850822150291592 0.0697157716842748
## x0032 4 Male 0 0.113977975025216
## ... ... ... ... ...
## x0018 4 Female 0.0170284419499822 0.0781697956157343
## x0057 4 Female 0 0.0797774108319114
## x0061 4 Female 0 0.164026603367578
## x0077 4 Female 0 0.11227306633587
## x0079 4 Female 0.0120147521727839 0.091365002951843
## CD4T CD8T Eos
## <numeric> <numeric> <numeric>
## x0017 0.204068686112657 0.100974056859719 0
## x0043 0.161712715578477 0.128772235710367 0
## x0036 0.155463683913806 0.12774174901185 0
## x0041 0.00732789124320035 0.0321860552286678 0
## x0032 0.222303988949171 0.0216089881163733 2.60503891663895e-18
## ... ... ... ...
## x0018 0.112735084403202 0.0679681631510725 0.00837769947260005
## x0057 0.111071701794334 0.00910488689520549 4.61887343424789e-18
## x0061 0.224202761557123 0.131252122776747 0.0085107393064941
## x0077 0.168056145550346 0.078405927598852 0.0613396950863682
## x0079 0.205830330328636 0.114753887947108 0
## Mono Neu
## <numeric> <numeric>
## x0017 0.0385653506040031 0.490936252606603
## x0043 0.049954190999866 0.491822214741198
## x0036 0.102710482176672 0.459631855532962
## x0041 0.0718279616580959 0.807749226994262
## x0032 0.0567246148936481 0.575614453750378
## ... ... ...
## x0018 0.0579534708330382 0.658600075099212
## x0057 0.101526022615966 0.686010320431362
## x0061 0.0382646576955094 0.429574723586544
## x0077 0.0583411264640133 0.515284487718098
## x0079 0.0750535276475847 0.513596840714489
```

brge_gexp is an `ExpressionSet`

with the expression data corresponding to an Affimetrix GeneChip Human Gene 2.0 ST Array. It contains 67528 features and 100 samples as well as two phenotypic variables (age and sex):

```
data(brge_gexp)
brge_gexp
```

```
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 67528 features, 100 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: x0001 x0002 ... x0139 (100 total)
## varLabels: age sex
## varMetadata: labelDescription
## featureData
## featureNames: TC01000001.hg.1 TC01000002.hg.1 ...
## TCUn_gl000247000001.hg.1 (67528 total)
## fvarLabels: transcript_cluster_id probeset_id ... notes (11 total)
## fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:
```

`lapply(pData(brge_gexp), table)`

```
## $age
##
## 4
## 99
##
## $sex
##
## Female Male
## 48 51
```

Finally, we will create a `GenomicRanges`

with our target region:

`targetRange <- GRanges("chr11:102600000-103300000")`

Let’s illustrate the use of this package by analyzing the effect of the sex in methylation and expression. First, a genome wide analysis will be performed and then the region funcionalities will be introduced. It should be noticed that methylation and expression analyses will include hypothesis testing and visualitzation.

This demonstration will show the main functions needed to perform a methylation analysis. A more exhaustive description of these and other auxiliary functions can be found at *MEAL* vignette.

The function `runPipeline`

can run different analyses to the same dataset: limma, DiffVar, bumphunter, blockFinder, DMRcate and RDA. This function accepts a `GenomicRatioSet`

, an `ExpressionSet`

or a `SummarizedExperiment`

, so it can analyze methylation and gene expression data. We can pass our variables of interest to function `runPipeline`

with the parameter `variable_names`

. `runPipeline`

also allows the inclusion of adjustment variables (such as cell count) with the parameter `covariable_names`

. When we are interested in a target region, we can apply RDA by passing a `GenomicRanges`

to the argument `range`

. Users can find a deeper explanation of this function in the **Methylation Analysis with MEAL** vignette.

As we are interested in evaluating the effect of sex, we will set `variable_names`

to “sex”. We will include cell counts as covaribles to adjust their effect and we will pass our target region to `range`

to compute RDA:

```
cellCounts <- colnames(colData(brge_methy))[3:9]
methRes <- runPipeline(set = brge_methy, variable_names = "sex", covariable_names = cellCounts, range = targetRange)
methRes
```

```
## Object of class 'ResultSet'
## . created with: runPipeline
## . sva: no
## . #results: 6 ( error: 1 )
## . featureData: 392277 probes x 19 variables
```

`names(methRes)`

```
## [1] "DiffMean" "DiffVar" "bumphunter" "blockFinder" "dmrcate"
## [6] "RDA"
```

The analysis generates a `ResultSet`

object containing the results of six different analyses: DiffMean, DiffVar, bumphunter, blockFinder, DMRcate and RDA. The six methods are based in the same linear model that we specified with the parameters `variable_names`

and `covariable_names`

. Consequently, the first step after running the pipeline will be evaluate how this model fits the data.

We will rely on the results of limma (DiffMean), as this method run the same linear model to all CpGs. A common way to evaluate the goodness of fit of these models is by making a QQplot on the p-values of the linear regression. This plot compares the distribution of p-values obtained in the analysis with a theoretical distribution. Our plot also shows the lambda, a measure of inflation. If lambda is higher than 1, the p-values are smaller than expected. If it is lower than 1, p-values are bigger than expected. In both causes, the most common cause is that we should include more variables in the model.

We can get a QQplot with the function `plot`

. When applied to a `ResultSet`

, we can make different plots depending on the method. We can select the name of the result with the parameter `rid`

and the kind of plot with the parameter `type`

. We will set `rid`

to “DiffMean” to use limma results and `type`

to “qq” to get a QQplot:

`plot(methRes, rid = "DiffMean", type = "qq")`

When there are few CpGs modicated, we expect points in the QQplot lying on the theoretical line. However, as all the CpGs from the sexual chromosomes will be differentially methylated between males and females, we obtained a S-shape. In this case, lambda is lower than 1 (0.851).

There are different ways to address this issue. An option is to include other covariates in our model. When we do not have more covariates, as it is our case, we can apply Surrogate Variable Analysis (SVA) to the data. SVA is a statistical technique that tries to determine hidden covariates or SV (Surrogate Variables) based on the measurements. This method is very useful to correct inflation when we do not have the variables that are missing in the model.

Our function `runPipeline`

includes the parameter `sva`

that runs SVA and includes these variables as covariates in the linear model. We do not need to include cell counts as covariates as their effect is be included in the SVs :

```
methRes <- runPipeline(set = brge_methy, variable_names = "sex", range = targetRange, sva = TRUE)
methRes
```

We have included two methods to analyze individual CpGs: DiffMean (limma) and DiffVar. DiffMean tests, for each CpG, if there is a difference between the means of the different groups. In our case, we tested, for each CpG, if the mean methylation in boys was different to the mean methylation in girls. On the other hand, DiffVar tests, for each CpG, if there is a difference in variance between the groups. In our case, we tested for each CpG if the variance was different in boys than in girls.

As both methods return results for each CpG, we can get an overview of the results using a Manhattan plot. A Manhattan plot represents the p-value of each CpG versus their location in the genome. This representation is useful to find genomic regions exhibiting differentiated behaviours.

We can get a Manhattan plot from a `ResultSet`

using the function `plot`

and setting `type`

to “manhattan”. We will make two plots, one for DiffMean results and the other for DiffVar. We will highlight the CpGs of our target region by passing a `GenomicRanges`

to highlight argument. `GenomicRanges`

passed to `plot`

should have chromosomes coded as number:

```
targetRangeNum <- GRanges("11:102600000-103300000")
plot(methRes, rid = "DiffMean", type = "manhattan", main = "Differences in Means", highlight = targetRangeNum)
```