---
title: "Chronological and gestational DNAm age estimation using different methylation-based clocks"
subtitle: "Dolors Pelegri and Juan R Gonzalez"
author: |
Institute for Global Health (ISGlobal), Barcelona, Spain
Bioinformatics Research Group in Epidemiolgy (BRGE)
http://brge.isglobal.org
date: "`r Sys.Date()`"
package: "`r pkg_ver('methylclock')`"
output:
BiocStyle::html_document:
number_sections: true
toc: yes
fig_caption: yes
vignette: >
%\VignetteIndexEntry{DNAm age using diffrent methylation clocks}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
bibliography: methylclock.bib
---
```{r setup_knitr, include=FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE,
cache=FALSE, comment = " ")
options(timeout = 300)
```
# Description of implemented clocks
This manual describes how to estimate chronological and gestational DNA
methylation (DNAm) age as well as biological age using different methylation
clocks. The package includes the following estimators:
## Chronological DNAm age (in years)
- **Horvath's clock**: It uses 353 CpGs described in @horvath2013dna. It was
trained using 27K and 450K arrays in samples from different tissues. Other
three different age-related biomarkers are also computed:
- **AgeAcDiff** (DNAmAge acceleration difference): Difference between
DNAmAge and chronological age.
- **IEAA** (Intrinsic Epigenetic Age Acceleration): Residuals obtained
after regressing DNAmAge and chronological age adjusted by cell counts.
- **EEAA** (Extrinsic Epigenetic Age Acceleration): Residuals obtained
after regressing DNAmAge and chronological age. This measure was also
known as DNAmAge acceleration residual in the first Horvath's paper.
Hannum's clock**: It uses 71 CpGs described in @hannum2013genome. It
was trained using 450K array in blood samples. Another are-related
biomarer is also computed:
- **AMAR** (Apparent Methylomic Aging Rate): Measure proposed in
@hannum2013genome computed as the ratio between DNAm age and the
chronological age.
- **BNN**: It uses Horvath's CpGs to train a Bayesian Neural Network
(BNN) to predict DNAm age as described in @alfonso2018.
- **Horvath's skin+blood clock (skinHorvath)**: Epigenetic clock for skin
and blood cells. It uses 391 CpGs described in @horvath2018epigenetic.
It was trained using 450K EPIC arrays in skin and blood sampels.
- **PedBE clock**: Epigenetic clock from buccal epithelial swabs. It's
intended purpose is buccal samples from individuals aged 0-20 years old.
It uses 84 CpGs described in @mcewen2019pedbe. The authors gathered 1,721
genome-wide DNAm profiles from 11 different cohorts with individuals aged
0 to 20 years old.
- **Wu's clock**: It uses 111 CpGs described in @wu2019dna. It is designed
to predict age in children. It was trained using 27K and 450K.
## Gestational DNAm age (in weeks)
- **Knight's clock**: It uses 148 CpGs described in @knight2016epigenetic.
It was trained using 27K and 450K arrays in cord blood samples.
- **Bohlin's clock**: It uses 96 CpGs described in @bohlin2016prediction.
It was trained using 450K array in cord blood samples.
- **Mayne's clock**: It uses 62 CpGs described in @mayne2017accelerated.
It was trained using 27K and 450K.
- **Lee's clocks**: Three different biological clocks described in
@lee2019placental are implemented. It was trained for 450K and EPIC
arrays in placenta samples.
- **RPC clock**: Robust placental clock (RPC). It uses 558 CpG sites.
- **CPC clock**: Control placental clock (CPC). It usses 546 CpG sites.
- **Refined RPC clock**: Useful for uncomplicated term pregnancies
(e.g. gestational age >36 weeks). It uses 396 CpG sites.
The biological DNAm clocks implemented in our package are:
- **Levine's clock** (also know as PhenoAge): It uses 513 CpGs described
in @levine2018epigenetic. It was trained using 27K, 450K and EPIC arrays
in blood samples.
- **Telomere Length's clock** (TL): It uses 140 CpGs described in @lu2019dna
It was trained using 450K and EPIC arrays in blood samples.
The main aim of this package is to facilitate the interconnection with R and
Bioconductor's infrastructure and, hence, avoiding submitting data to online
calculators. Additionally, `methylclock` also provides an unified way of
computing DNAm age to help downstream analyses.
# Getting started
The package depends on some R packages that can be previously installed into
your computer by:
```{r install_req_packages, eval=FALSE}
install.packages(c("tidyverse", "impute", "Rcpp"))
```
Then `methylclock` package is installed into your computer by executing:
```{r install_packages, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("methylclock")
```
The package is loaded into R as usual:
```{r load_package}
library(methylclockData)
library(methylclock)
```
These libraries are required to reproduce this document:
```{r load_others, eval=TRUE}
library(Biobase)
library(tibble)
library(impute)
library(ggplot2)
library(ggpmisc)
library(GEOquery)
```
# DNA Methylation clocks
The main function to estimate chronological and biological mDNA age is called
`DNAmAge` while the gestational DNAm age is estimated using `DNAmGA` function.
Both functions have similar input arguments. Next subsections detail some of
the important issues to be consider before computind DNAm clocks.
## Data format
The methylation data is given in the argument `x`. They can be either beta or
M values. The argument `toBetas` should be set to TRUE when M values are
provided. The `x` object can be:
- A **matrix** with CpGs in rows and individuals in columns having the name
of the CpGs in the rownames.
- A **data frame** or a **tibble** with CpGs in rows and individuals in columns
having the name of the CpGs in the first column (e.g. cg00000292, cg00002426,
cg00003994, ...) as required in the Horvath's DNA Methylation Age Calculator
website (https://dnamage.genetics.ucla.edu/home).
- A **GenomicRatioSet** object, the default method to encapsulate methylation
data in `minfi` Bioconductor package.
- An **ExpressionSet** object as obtained, for instance, when downloading
methylation data from GEO (https://www.ncbi.nlm.nih.gov/geo/).
## Data nomalization
In principle, data can be normalized by using any of the existing standard
methods such as QN, ASMN, PBC, SWAN, SQN, BMIQ (see a revision of those
methods in @wang2015systematic). `DNAmAge` function includes the BMIQ method
proposed by @teschendorff2012beta using Horvath's robust implementation that
basically consists of an optimal R code implementation and optimization
procedures. This normalization is recommended by Horvath since it improves
the predictions for his clock. This normalization procedure is very
time-consuming. In order to overcome these difficulties, we have parallelize
this process using `BiocParallel` library. This step is not mandatory, so that,
you can use your normalized data and set the argument `normalize` equal to
FALSE (default).
## Missing individual's data
All the implemented methods require complete cases. `DNAmAge` function has an
imputation method based on KNN implemented in the function `knn.impute` from
`impute` Bioconductor package. This is performed when missing data is present
in the CpGs used in any of the computed clocks. There is also another option
based on a fast imputation method that imputes missing values by the median of
required CpGs as recommended in @bohlin2016prediction. This is recommended when
analyzing 450K arrays since `knn.impute` for large datasets may be very time
consuming. Fast imputation can be performed by setting `fastImp=TRUE` which is
not the default value.
## Missing CpGs of DNAm clocks {#section-missingCpGs}
By default the package computes the different clocks when there are more than
80% of the required CpGs of each method. Nothing is required when having
missing CpGs since the main functions will return NA for those estimators
when this criteria is not meet. Let us use a test dataset (`TestDataset`)
which is available within the package to illustrate the type of information
we are obtaining:
```{r check}
# Get TestDataset data
TestDataset <- get_TestDataset()
cpgs.missing <- checkClocks(TestDataset)
```
```{r checkGA}
cpgs.missing.GA <- checkClocksGA(TestDataset)
```
The objects `cpgs.missing` and `cpgs.missing.GA` are lists having the missing
CpGs of each clock
```{r showMissing}
names(cpgs.missing)
```
We can see which are those CpGs for a given clock (for example Hannum) with
the function `commonClockCpgs` :
```{r showMissNames}
commonClockCpgs(cpgs.missing, "Hannum" )
commonClockCpgs(cpgs.missing.GA, "Bohlin" )
```
In Section \@ref(section-example) we describe how to change this 80% threshold.
## Cell counts
The EEAA method requires to estimate cell counts. We use the package `meffil`
(@min2018meffil) that provides some functions to estimate cell counts using
predefined datasets. This is performed by setting `cell.count=TRUE` (default
value). The reference panel is passed through the argument
`cell.count.reference`. So far, the following options are available:
- **"blood gse35069 complete"**: methylation profiles from
@reinius2012differential for purified blood cell types. It includes CD4T, CD8T,
Mono, Bcell, NK, Neu and Eos.
- **"blood gse35069"**: methylation profiles from @reinius2012differential for
purified blood cell types. It includes CD4T, CD8T, Mono, Bcell, NK and Gran.
- **"blood gse35069 chen"**: methylation profiles from @chen2017epigenome blood
cell types. It includes CD4T, CD8T, Mono, Bcell, NK, Neu and Eos.
- **"andrews and bakulski cord blood"**. Cord blood reference from
@bakulski2016dna. It includes Bcell, CD4T, CD8T, Gran, Mono, NK and nRBC.
- **"cord blood gse68456"** Cord blood methylation profiles from
@de2015nucleated. It includes CD4T, CD8T, Mono, Bcell, NK, Neu, Eos and RBC.
- **"gervin and lyle cord blood"** Cord blood reference generated by Kristina
Gervin and Robert Lyle, available at `miffil` package. It includes CD14, Bcell,
CD4T, CD8T, NK, Gran.
- **"saliva gse48472"**: Reference generated from the multi-tissue pannel from
@slieker2013identification. It includes Buccal, CD4T, CD8T, Mono, Bcell, NK,
Gran.
- **"guintivano dlpfc"**: Reference generated from @guintivano2013cell. It
includes dorsolateral prefrontal cortex, NeuN_neg and NeuN_pos.
- **"combined cord blood"**: References generated based in samples assayed by
Bakulski et al, Gervin et al., de Goede et al., and Lin et al. It includes
umbilical cord blood, Bcell, CD4T, CD8T, Gran, Mono, NK and nRBC
# Chronological and biological DNAm age estimation
Next we illustrate how to estimate the chronological DNAm age using several
datasets which aim to cover different data input formats.
**IMPORTANT NOTE**: On some systems we can find an error in the `DNAmAge()`
function when parameter `cell.count = TRUE`. This error is related to
`preprocessCore` package and can be fixed by disabling multi-threading
when installing the preprocessCore package using the command
```
BiocManager::install("preprocessCore",
configure.args = "--disable-threading",
force = TRUE)
```
## Data in Horvath's format (e.g. `csv` with CpGs in rows) {#section-example}
Let us start by reproducing the results proposed in @horvath2013dna. It uses
the format available in the file 'MethylationDataExample55.csv" from his
tutorial (available [here](https://dnamage.genetics.ucla.edu/home)). These data
are available at `methylclock` package. Although these data can be loaded into
R by using standard functions such as `read.csv` we hihgly recommend to use
functions from `tidiverse`, in particular `read_csv` from `readr` package.
The main reason is that currently researchers are analyzing Illumina 450K
or EPIC arrays that contains a huge number of CpGs that can take a long time
to be loaded when using basic importing R function. These functions import
`csv` data as tibble which is one of the possible formats of `DNAmAge` function
```{r load_horvath_example}
library(tidyverse)
MethylationData <- get_MethylationDataExample()
MethylationData
```
*IMPORTANT NOTE*: Be sure that the first column contains the CpG names.
Sometimes, your imported data look like this one (it can happen, for
instance, if the `csv` file was created in R without indicating
`row.names=FALSE`)
```
> mydata
# A tibble: 473,999 x 6
X1 Row.names BIB_15586_1X BIB_33043_1X EDP_5245_1X KAN_584_1X
<int> <chr> <dbl> <dbl> <dbl> <dbl>
1 1 cg000000~ 0.635 0.575 0.614 0.631
2 2 cg000001~ 0.954 0.948 0.933 0.950
3 3 cg000001~ 0.889 0.899 0.901 0.892
4 4 cg000001~ 0.115 0.124 0.107 0.123
5 5 cg000002~ 0.850 0.753 0.806 0.815
6 6 cg000002~ 0.676 0.771 0.729 0.665
7 7 cg000002~ 0.871 0.850 0.852 0.863
8 8 cg000003~ 0.238 0.174 0.316 0.206
```
If so, the first column must be removed before being used as the input
object in `DNAmAge` funcion. It can be done using `dplyr` function
```
> mydata2 <- select(mydata, -1)
# A tibble: 473,999 x 5
Row.names BIB_15586_1X BIB_33043_1X EDP_5245_1X KAN_584_1X
<chr> <dbl> <dbl> <dbl> <dbl>
1 cg000000~ 0.635 0.575 0.614 0.631
2 cg000001~ 0.954 0.948 0.933 0.950
3 cg000001~ 0.889 0.899 0.901 0.892
4 cg000001~ 0.115 0.124 0.107 0.123
5 cg000002~ 0.850 0.753 0.806 0.815
6 cg000002~ 0.676 0.771 0.729 0.665
7 cg000002~ 0.871 0.850 0.852 0.863
8 cg000003~ 0.238 0.174 0.316 0.206
```
In any case, if you use the object `mydata` that contains the CpGs in the
second column, you will see this error message:
```
> DNAmAge(mydata)
Error in DNAmAge(mydata) : First column should contain CpG names
```
Once data is in the proper format, DNAmAge can be estimated by simply:
```{r DNAmAge_horvath, warning=TRUE}
age.example55 <- DNAmAge(MethylationData)
age.example55
```
As mention in Section \@ref(section-missingCpGs) some clocks returns NA when
there are more than 80% of the required CpGs are missing as we can see when
typing
```{r show_cpg_miss}
missCpGs <- checkClocks(MethylationData)
```
Here we can observe that 72.1% of the required CpGs for SkinHorvath clock are
missing. We could estimate DNAm age using this clock just changing the argument
`min.perc` in `DNAmAge`. For example, we can indicate that the minimum amount
of required CpGs for computing a given clock should be 25%.
```{r DNAmAgemp_horvath, warning=TRUE}
age.example55.2 <- DNAmAge(MethylationData, min.perc = 0.25)
age.example55.2
```
In that case, we see as SkinHorvath clock is estimated (though it can be
observed that the estimation is not very accurate - this is why we considered
at least having 80% of the required CpGs).
By default all available clocks (Hovarth, Hannum, Levine, BNN, skinHorvath,...)
are estimated. One may select a set of clocks by using the argument `clocks`
as follows:
```{r DNAmAge_horvath_sel, warning=TRUE}
age.example55.sel <- DNAmAge(MethylationData, clocks=c("Horvath", "BNN"))
age.example55.sel
```
## Age acceleration
However, in epidemiological studies one is interested in assessing whether age
acceleration is associated with a given trait or condition. Three different
measures can be computed:
- **ageAcc**: Difference between DNAmAge and chronological age.
- **ageAcc2**: Residuals obtained after regressing chronological age and
DNAmAge (similar to IEAA).
- **ageAcc3**: Residuals obtained after regressing chronological age and
DNAmAge adjusted for cell counts (similar to EEAA).
All this estimates can be obtained for each clock when providing chronological
age through `age` argument. This information is normally provided in a
different file including different covariates (metadata or sample annotation
data). In this example data are available at 'SampleAnnotationExample55.csv'
file that is also available at `methylclock` package:
```{r covariates_horvath_example}
library(tidyverse)
path <- system.file("extdata", package = "methylclock")
covariates <- read_csv(file.path(path, "SampleAnnotationExample55.csv"))
covariates
```
In this case, chronological age is available at `Age` column:
```{r age_horvath_example}
age <- covariates$Age
head(age)
```
The different methylation clocks along with their age accelerated estimates
can be simply computed by:
```{r DNAmAge_horvath_cell, warning=TRUE}
age.example55 <- DNAmAge(MethylationData, age=age, cell.count=TRUE)
age.example55
```
By default, the argument `cell.count` is set equal to TRUE and, hence, can
be omitted. This implies that `ageAcc3` will be computed for all clocks.
In some occassions this can be very time consuming. In such cases one can
simply estimate DNAmAge, accAge and accAge2 by setting `cell.count=FALSE`.
NOTE: see section 3.5 to see the reference panels available to estimate
cell counts.
Then, we can investigate, for instance, whether the accelerated age is
associated with Autism. In that example we will use a non-parametric
test (NOTE: use t-test or linear regression for large sample sizes)
```{r compare_autistic}
autism <- covariates$diseaseStatus
kruskal.test(age.example55$ageAcc.Horvath ~ autism)
kruskal.test(age.example55$ageAcc2.Horvath ~ autism)
kruskal.test(age.example55$ageAcc3.Horvath ~ autism)
```
## Chronological age prediction using `ExpressionSet` data
One may be interested in assessing association between chronologial age and DNA
methylation age or evaluating how well chronological age is predicted by
DNAmAge. In order to illustrate this analysis we downloaded data from GEO
corresponding to a set of healthy individuals (GEO accession number GSE58045).
Data can be retrieved into R by using `GEOquery` package as an `ExpressionSet`
object that can be the input of our main function.
```{r get_gse58045, echo=FALSE}
# ff <- "c:/juan/CREAL/BayesianPrediction/Bayesian_clock/paper"
# load(file.path(ff, "data/GSE58045.Rdata"))
```
```{r get_geo_gse58045, eval=TRUE}
# To avoid connection buffer size
Sys.setenv("VROOM_CONNECTION_SIZE" = 131072 * 10)
# Download data
dd <- GEOquery::getGEO("GSE58045")
gse58045 <- dd[[1]]
# Restore connection buffer size
Sys.setenv("VROOM_CONNECTION_SIZE" = 131072)
```
```{r show_gse58045}
gse58045
```
The chronological age is obtained by using `pData` function from `Biobase`
package that is able to deal with `ExpressionSet` objects:
```{r age_gse58045}
pheno <- pData(gse58045)
age <- as.numeric(pheno$`age:ch1`)
```
And the different DNA methylation age estimates are obtained by using `DNAmAge`
function (NOTE: as there are missing values, the program automatically runs
`impute.knn` function to get complete cases):
```{r DNAmAge_gse58045, warning=TRUE}
age.gse58045 <- DNAmAge(gse58045, age=age)
age.gse58045
```
Figure \ref{fig:horvath_age} shows the correlation between DNAmAge obtained
from Horvath's method and the chronological age, while Figure \ref{fig:bnn_age}
depicts the correlation of a new method based on fitting a Bayesian Neural
Network to predict DNAmAge based on Horvath's CpGs.
```{r horvat_age}
plotDNAmAge(age.gse58045$Horvath, age)
```
```{r bnn_age}
plotDNAmAge(age.gse58045$BNN, age, tit="Bayesian Neural Network")
```
## Use of DNAmAge in association studies
Let us illustrate how to use DNAmAge information in association studies
(e.g case/control, smokers/non-smokers, responders/non-responders, ...).
GEO number GSE19711 contains transcriptomic and epigenomic data of a study
in lung cancer. Data can be retrieved into R by
```{r get_gse19711, echo=FALSE}
# ff <- "c:/juan/CREAL/BayesianPrediction/Bayesian_clock/paper"
# load(file.path(ff, "data/GSE19711.Rdata"))
```
```{r get_geo_gse19711, eval=TRUE}
# To avoid connection buffer size
Sys.setenv("VROOM_CONNECTION_SIZE" = 131072 * 10)
# Download data
dd <- GEOquery::getGEO("GSE19711")
gse19711 <- dd[[1]]
# Restore connection buffer size
Sys.setenv("VROOM_CONNECTION_SIZE" = 131072)
```
The object `gse19711`is an `ExpressionSet` that can contains CpGs and
phenotypic (e.g clinical) information
```{r show_gse19711}
gse19711
```
Let us imagine we are interested in comparing the accelerated age between
cases and controls. Age and case/control status information can be obtained by:
```{r get_case_control}
pheno <- pData(gse19711)
age <- as.numeric(pheno$`ageatrecruitment:ch1`)
disease <- pheno$`sample type:ch1`
table(disease)
disease[grep("Control", disease)] <- "Control"
disease[grep("Case", disease)] <- "Case"
disease <- factor(disease, levels=c("Control", "Case"))
table(disease)
```
The DNAmAge estimates of different methods is computed by
```{r DNAmAge_gse19711, warning=TRUE}
age.gse19711 <- DNAmAge(gse19711, age=age)
```
We can observe there are missing data. The funcion automatically impute those
using `impute.knn` function from `impute` package since complete cases are
required to compute the different methylation clocks. The estimates are:
```{r show_age.gse19711}
age.gse19711
```
The association between disease status and DNAmAge estimated using Horvath's
method can be computed by
```{r assoc_hpv}
mod.horvath1 <- glm(disease ~ ageAcc.Horvath ,
data=age.gse19711,
family="binomial")
summary(mod.horvath1)
mod.skinHorvath <- glm(disease ~ ageAcc2.Horvath ,
data=age.gse19711,
family="binomial")
summary(mod.skinHorvath)
mod.horvath3 <- glm(disease ~ ageAcc3.Horvath ,
data=age.gse19711,
family="binomial")
summary(mod.horvath3)
```
We do not observe statistical significant association between age acceleration
estimated using Horvath method and the risk of developing lung cancer. It is
worth to notice that Horvath's clock was created to predict chronological age
and the impact of age acceleration of this clock on disease may be limited.
On the other hand, Levine's clock aimed to distinguish risk between same-aged
individuals. Let us evaluate whether this age acceleration usin Levine's clock
is associated with lung cancer
```{r mod_levine}
mod.levine1 <- glm(disease ~ ageAcc.Levine , data=age.gse19711,
family="binomial")
summary(mod.levine1)
mod.levine2 <- glm(disease ~ ageAcc2.Levine , data=age.gse19711,
family="binomial")
summary(mod.levine2)
mod.levine3 <- glm(disease ~ ageAcc3.Levine , data=age.gse19711,
family="binomial")
summary(mod.levine3)
```
Here we observe as the risk of developing lung cancer increases
`r round((exp(coef(mod.levine1)[2]) - 1)*100,2)` percent per each unit in the
age accelerated variable (`ageAcc`). Similar conclusion is obtained when using
`ageAcc2` and `ageAcc3` variables.
In some occasions cell composition should be used to assess association. This
information is calculated in `DNAmAge` function and it can be incorporated in
the model by:
```{r assoc_cell}
cell <- attr(age.gse19711, "cell_proportion")
mod.cell <- glm(disease ~ ageAcc.Levine + cell, data=age.gse19711,
family="binomial")
summary(mod.cell)
```
Here we observe as the positive association disapears after adjusting
for cell counts.
## Use of DNAm age in children
```{r get_gse109446, echo=FALSE}
# ff <- "c:/juan/CREAL/BayesianPrediction/Bayesian_clock/paper"
# load(file.path(ff, "data/GSE109446.Rdata"))
```
```{r get_geo_109446, eval=TRUE}
dd <- GEOquery::getGEO("GSE109446")
gse109446 <- dd[[1]]
```
```{r age_gse109446, warning=TRUE}
controls <- pData(gse109446)$`diagnosis:ch1`=="control"
gse <- gse109446[,controls]
age <- as.numeric(pData(gse)$`age:ch1`)
age.gse <- DNAmAge(gse, age=age)
```
```{r plotClocks}
plotCorClocks(age.gse)
```
# Gestational DNAm age estimation
## Model predicion
Let us start by reproducing the example provided in @knight2016epigenetic as
a test data set (file 'TestDataset.csv'). It consists on 3 individuals whose
methylation data are available as supplementary data of their paper. The data
is also available at `methylclock` package as a data frame.
```{r load_3_inds}
TestDataset[1:5,]
```
The Gestational Age (in months) is simply computed by
```{r age_test, warning=TRUE}
ga.test <- DNAmGA(TestDataset)
ga.test
```
like in DNAmAge we can use the parameter `min.perc` to set the minimum missing
percentage.
The results are the same as those described in the additional file 7
of @knight2016epigenetic (link [here]
(https://static-content.springer.com/esm/art%3A10.1186%2Fs13059-016-1068-z/MediaObjects/13059_2016_1068_MOESM7_ESM.docx))
Let us continue by illustrating how to compute GA of real examples.
The PROGRESS cohort data is available in the additional file 8 of
@knight2016epigenetic. It is available at `methylclock` as a `tibble`:
```{r get_progress}
data(progress_data)
```
This file also contains different variables that are available in this
`tibble`.
```{r progressClin}
data(progress_vars)
```
The Clinical Variables including clinical assesment of gestational age
(EGA) are available at this `tibble`.
The Gestational Age (in months) is simply computed by
```{r age_progress, warning=TRUE}
ga.progress <- DNAmGA(progress_data)
ga.progress
```
We can compare these results with the clinical GA available in the variable EGA
```{r plot_progress}
plotDNAmAge(ga.progress$Knight, progress_vars$EGA,
tit="GA Knight's method",
clock="GA")
```
Figure 3b (only for PROGRESS dataset) in @knight2016epigenetic representing
the correlation between GA acceleration and birthweight can be reproduced by
```{r plotAcc}
library(ggplot2)
progress_vars$acc <- ga.progress$Knight - progress_vars$EGA
p <- ggplot(data=progress_vars, aes(x = acc, y = birthweight)) +
geom_point() +
geom_smooth(method = "lm", se=FALSE, color="black") +
xlab("GA acceleration") +
ylab("Birthweight (kgs.)")
p
```
Finally, we can also estimate the "accelerated gestational age" using two
of the the three different estimates previously described (`accAge`, `accAge2`)
by provinding information of gestational age through `age` argument. Notice
that in that case `accAge3` cannot be estimates since we do not have all the
CpGs required by the default reference panel to estimate cell counts for
gestational age which is "andrews and bakulski cord blood".
```{r acccelerated_ga, warning=TRUE}
accga.progress <- DNAmGA(progress_data,
age = progress_vars$EGA,
cell.count=FALSE)
accga.progress
```
One can also check which clocks can be estimated given the CpGs available
in the methylation data by
```{r check_clocks_2}
checkClocksGA(progress_data)
```
# Correlation among DNAm clocks
We can compute the correlation among biological clocks using the function
`plotCorClocks` that requires the package `ggplot2` and `ggpubr` to be
installed in your computer.
We can obtain, for instance, the correlation among the clocks estimated for
the healthy individuals study previosuly analyze (GEO accession number
GSE58045) by simply executing:
```{r plotCorClockHealth}
plotCorClocks(age.gse58045)
```
# References
```{sessioninfo}
sessionInfo()
```