Script 46: Analysing epigenetic clocks
Compiled: June 29, 2024
This script analyses epigenetic age changes and outputs Supplementary Tables
Setup
Load packages
library(tidyverse)
library(lme4)
library(lmerTest)
library(SummarizedExperiment)Load targets
load('../GOTO_Data/GOTO_targets-filtered.Rdata')Create actual age (+0.25 post intervention)
targets$age <- ifelse(targets$timepoint==1,
targets$age + 0.25,
targets$age)Save muscle and blood targets separately
targets_blood <- targets %>% filter(tissue == 'fasted blood')
targets_muscle <- targets %>% filter(tissue == 'muscle')cAge correlations
Bernabeu CA
load('../GOTO_Data/Clocks/Clock_Out/BernabeuE2023c_blood.Rdata')
bern_CA <- age %>%
select(Basename=Sample, bern_CA = mAge)
targets_cage <- full_join(targets_blood,
bern_CA,
by='Basename')Correlation test
cor <- cor.test(targets_cage$bern_CA, targets_cage$age)
cor##
## Pearson's product-moment correlation
##
## data: targets_cage$bern_CA and targets_cage$age
## t = 33.349, df = 194, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8988419 0.9411840
## sample estimates:
## cor
## 0.9227509out <- data.frame(
clock = 'bern_CA',
r = cor$estimate,
low_ci = cor$conf.int[1],
upp_ci = cor$conf.int[2],
p = cor$p.value
)Horvath
load('../GOTO_Data/Clocks/Clock_Out/HorvathS2013_blood.Rdata')
horvath <- age %>% select(Basename=Sample, horvath = mAge)
targets_cage <- left_join(targets_cage, horvath, by='Basename')Correlation test
cor <- cor.test(targets_cage$horvath, targets_cage$age)
cor##
## Pearson's product-moment correlation
##
## data: targets_cage$horvath and targets_cage$age
## t = 14.858, df = 194, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6565406 0.7890357
## sample estimates:
## cor
## 0.729562out <- rbind(out, data.frame(
clock = 'Horvath',
r = cor$estimate,
low_ci = cor$conf.int[1],
upp_ci = cor$conf.int[2],
p = cor$p.value
))Zhang
load('../GOTO_Data/Clocks/Clock_Out/ZhangQ2019_blood.Rdata')
zhang <- age %>% select(Basename=Sample, zhang = mAge)
targets_cage <- left_join(targets_cage, zhang, by='Basename')Correlation test
cor <- cor.test(targets_cage$zhang, targets_cage$age)
cor##
## Pearson's product-moment correlation
##
## data: targets_cage$zhang and targets_cage$age
## t = 26.94, df = 194, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8545315 0.9145865
## sample estimates:
## cor
## 0.8882973out <- rbind(out, data.frame(
clock = 'Zhang',
r = cor$estimate,
low_ci = cor$conf.int[1],
upp_ci = cor$conf.int[2],
p = cor$p.value
))Multiple testing
out$padj <- p.adjust(out$p, method='fdr')out## clock r low_ci upp_ci p padj
## cor bern_CA 0.9227509 0.8988419 0.9411840 2.872001e-82 8.616003e-82
## cor1 Horvath 0.7295620 0.6565406 0.7890357 7.636092e-34 7.636092e-34
## cor2 Zhang 0.8882973 0.8545315 0.9145865 1.775479e-67 2.663219e-67write_csv(out, file='../GOTO_Data/Tables/ST22.csv')bAge predictions
Bernabeu BA and grimAge
bern_BA <- read_tsv('../GOTO_Data/Clocks/Clock_Out/bage_predictions.tsv')## Rows: 534 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): Sample
## dbl (6): bAge, bAge_Years, GrimAge, Age, Sex, Dead
## lgl (1): TTE
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.bern_BA <- bern_BA %>%
select(Basename=Sample, bern_BA = bAge_Years, grimAge = GrimAge)
targets_bage <- left_join(targets_blood, bern_BA, by='Basename')Linear mixed model - Bernabeu bAge
fit <- lmer(bern_BA ~ timepoint + age + sex +
plate + array_row + (1|IOP2_ID),
data = targets_bage)
summary(fit)$coefficients## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.66798092 4.89465260 88.01688 2.99673585 3.545228e-03
## timepoint1 -0.30923660 0.26389316 101.53494 -1.17182502 2.440097e-01
## age 0.78032959 0.07536386 87.17455 10.35416191 7.511724e-17
## sex1 -0.33980276 0.89469994 87.15471 -0.37979522 7.050213e-01
## plateM01004_03 2.10261160 1.60969648 87.93212 1.30621619 1.948856e-01
## plateM01004_04 -0.11115421 1.54471875 87.23191 -0.07195757 9.428004e-01
## plateM01004_05 -2.46015440 1.56704456 87.26122 -1.56993263 1.200500e-01
## plateM01004_06 -0.88437497 1.62804709 87.19333 -0.54321216 5.883698e-01
## plateM01004_07 1.64037724 2.07051965 87.26466 0.79225388 4.303611e-01
## plateM01004_08 -0.13444865 1.91731256 87.42903 -0.07012349 9.442555e-01
## plateM01044-01 2.83056718 2.08146693 87.12454 1.35989054 1.773733e-01
## array_row 0.04165137 0.15446383 162.32992 0.26965130 7.877710e-01out <- data.frame(clock = "bern_BA",
delta = summary(fit)$coefficients[2,1],
weeks = 52*summary(fit)$coefficients[2,1],
SE = summary(fit)$coefficients[2,2],
p = summary(fit)$coefficients[2,5])Linear mixed model - grimAge
fit <- lmer(grimAge ~ timepoint + age + sex +
plate + array_row + (1|IOP2_ID), data = targets_bage)
summary(fit)$coefficients## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 16.81234550 3.55967535 88.60635 4.72299967 8.660998e-06
## timepoint1 -0.66575254 0.15420609 104.26335 -4.31729074 3.607808e-05
## age 0.70790064 0.05486995 87.73465 12.90142663 5.662755e-22
## sex1 -3.51955167 0.65141970 87.71200 -5.40289413 5.560795e-07
## plateM01004_03 0.05035535 1.17075283 88.56990 0.04301109 9.657895e-01
## plateM01004_04 0.44768767 1.12457093 87.79697 0.39809643 6.915255e-01
## plateM01004_05 -1.26406084 1.14077844 87.82925 -1.10806866 2.708574e-01
## plateM01004_06 -1.38823274 1.18529760 87.75450 -1.17121028 2.446846e-01
## plateM01004_07 -0.40621187 1.50729155 87.83304 -0.26949788 7.881787e-01
## plateM01004_08 -0.64975366 1.39544561 88.01420 -0.46562450 6.426345e-01
## plateM01044-01 1.66222714 1.51555258 87.67877 1.09677959 2.757410e-01
## array_row 0.09620557 0.10278289 180.15453 0.93600764 3.505223e-01out <- rbind(out, data.frame(clock = "grimAge",
delta = summary(fit)$coefficients[2,1],
weeks = 52*summary(fit)$coefficients[2,1],
SE = summary(fit)$coefficients[2,2],
p = summary(fit)$coefficients[2,5]))phenoAge
load('../GOTO_Data/Clocks/Clock_Out/LevineM2018_blood.Rdata')
phenoAge <- age %>% select(Basename=Sample, phenoAge = mAge)
targets_bage <- left_join(targets_bage, phenoAge, by='Basename')Linear mixed model - phenoAge
fit <- lmer(phenoAge ~ timepoint + age + sex +
plate + array_row + (1|IOP2_ID), data = targets_bage)
summary(fit)$coefficients## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.43924184 6.28899308 88.36135 -0.06984295 9.444765e-01
## timepoint1 -0.77879587 0.37509487 101.16596 -2.07626374 4.040462e-02
## age 0.76030199 0.09678518 87.54321 7.85556218 9.393428e-12
## sex1 -2.13812814 1.14899524 87.52501 -1.86086772 6.612032e-02
## plateM01004_03 5.13738996 2.06812931 88.24541 2.48407579 1.487731e-02
## plateM01004_04 2.90494148 1.98385295 87.59664 1.46429274 1.466937e-01
## plateM01004_05 1.51162864 2.01255941 87.62382 0.75109765 4.546064e-01
## plateM01004_06 3.68351144 2.09082368 87.56085 1.76175135 8.160169e-02
## plateM01004_07 2.93318632 2.65917902 87.62702 1.10304207 2.730290e-01
## plateM01004_08 1.21450735 2.46264619 87.77942 0.49317168 6.231228e-01
## plateM01044-01 7.16559439 2.67302370 87.49702 2.68070739 8.779114e-03
## array_row -0.01495657 0.20508039 152.48197 -0.07293028 9.419572e-01out <- rbind(out, data.frame(clock = "phenoAge",
delta = summary(fit)$coefficients[2,1],
weeks = 52*summary(fit)$coefficients[2,1],
SE = summary(fit)$coefficients[2,2],
p = summary(fit)$coefficients[2,5]))MEAT
load('../GOTO_Data/Clocks/Clock_Out/MEAT.Rdata')
meat <- as.data.frame(colData(methData_age)) %>%
select(meat = DNAmage) %>%
rownames_to_column(var='Basename')
targets_muscle <- left_join(targets_muscle,
meat, by='Basename')Linear mixed model - MEAT
fit <- lmer(meat ~ timepoint + age + sex +
plate + array_row + (1|IOP2_ID),
data = targets_muscle)
summary(fit)$coefficients## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.3055361 4.88620439 68.64571 0.8811617 3.813052e-01
## timepoint1 -1.1125197 0.42299800 78.62866 -2.6300825 1.026640e-02
## age 0.7998256 0.07346348 68.60834 10.8873900 1.343227e-16
## sex1 1.5271132 0.78532777 68.25473 1.9445551 5.595055e-02
## plateM01004_03 1.0641267 1.26922821 76.45672 0.8384045 4.044171e-01
## plateM01004_04 1.9217653 1.16430267 78.56514 1.6505719 1.028184e-01
## plateM01004_05 2.3413503 1.13741965 68.96361 2.0584753 4.332477e-02
## plateM01004_06 0.4435046 1.25936483 68.51311 0.3521653 7.257953e-01
## plateM01004_07 -0.8746515 1.64277610 68.24823 -0.5324229 5.961624e-01
## plateM01004_08 2.3477638 1.68046093 68.40013 1.3970951 1.669010e-01
## plateM01044-01 -2.6189302 3.37150670 68.21464 -0.7767833 4.399720e-01
## array_row 0.1288987 0.15112074 93.05455 0.8529520 3.958761e-01out <- rbind(out, data.frame(clock = "MEAT",
delta = summary(fit)$coefficients[2,1],
weeks = 52*summary(fit)$coefficients[2,1],
SE = summary(fit)$coefficients[2,2],
p = summary(fit)$coefficients[2,5]))out$padj <- p.adjust(out$p, method='fdr')
out## clock delta weeks SE p padj
## 1 bern_BA -0.3092366 -16.08030 0.2638932 2.440097e-01 0.2440096777
## 2 grimAge -0.6657525 -34.61913 0.1542061 3.607808e-05 0.0001443123
## 3 phenoAge -0.7787959 -40.49739 0.3750949 4.040462e-02 0.0538728252
## 4 MEAT -1.1125197 -57.85102 0.4229980 1.026640e-02 0.0205328080write_csv(out, file='../GOTO_Data/Tables/ST23.csv')Health associations
grimAge
Linear mixed models
traits <- c('bmi',
'wc',
'perc_body_fat',
'fasting_insulin',
'leptin',
'adiponectin',
'interleukin_6',
'hdl_size_fasted',
'hdl_c',
'systolic_bp')for(i in traits){
fit <- lmer(grimAge ~ get(i) + age + sex +
plate + array_row + (1|IOP2_ID),
data = targets_bage)
summary(fit)$coefficients
out <- data.frame(clock = 'grim',
trait = i,
estimate = summary(fit)$coefficients[2,1],
se = summary(fit)$coefficients[2,2],
p = summary(fit)$coefficients[2,5]
)
if(i == 'bmi'){
res <- out
} else {
res <- rbind(res, out)
}
}MEAT
for(i in traits){
fit <- lmer(meat ~ get(i) + age + sex +
plate + array_row + (1|IOP2_ID),
data = targets_muscle)
summary(fit)$coefficients
out <- data.frame(clock = 'meat',
trait = i,
estimate = summary(fit)$coefficients[2,1],
se = summary(fit)$coefficients[2,2],
p = summary(fit)$coefficients[2,5]
)
res <- rbind(res, out)
}Multiple testing
res$padj <- p.adjust(res$p, method='fdr')
res <- res %>% arrange(padj)
res## clock trait estimate se p padj
## 1 grim bmi 0.298275858 0.08440412 0.0005210228 0.005210228
## 2 grim wc 0.073299153 0.01990162 0.0003252913 0.005210228
## 3 grim leptin 0.084663864 0.02684840 0.0018894368 0.012596245
## 4 grim perc_body_fat 0.101323268 0.04057253 0.0136891150 0.039111757
## 5 grim fasting_insulin 0.106354374 0.04113461 0.0106502388 0.039111757
## 6 grim interleukin_6 0.302687286 0.11794374 0.0112830503 0.039111757
## 7 grim hdl_size_fasted -3.220168552 1.28356099 0.0130301673 0.039111757
## 8 grim adiponectin -0.124870296 0.05717670 0.0302481546 0.075620386
## 9 grim hdl_c -0.904859426 0.56306187 0.1098078495 0.244017443
## 10 meat leptin 0.061894568 0.05352251 0.2498657719 0.499731544
## 11 meat hdl_c 1.104677171 1.07175815 0.3046700206 0.553945492
## 12 meat interleukin_6 0.197094319 0.21852485 0.3687356803 0.614559467
## 13 grim systolic_bp 0.006290725 0.01144853 0.5835009023 0.766642179
## 14 meat wc 0.018112016 0.04003315 0.6516458523 0.766642179
## 15 meat perc_body_fat 0.046493496 0.07977277 0.5609375028 0.766642179
## 16 meat fasting_insulin -0.039210239 0.07860347 0.6187073597 0.766642179
## 17 meat systolic_bp -0.009934952 0.02165237 0.6470550761 0.766642179
## 18 meat bmi 0.043891781 0.14302521 0.7595301782 0.843922420
## 19 meat adiponectin -0.014798179 0.08971859 0.8693466150 0.889945740
## 20 meat hdl_size_fasted 0.315509250 2.27562018 0.8899457400 0.889945740write_csv(out, file='../GOTO_Data/Tables/ST24.csv')Plots
targets_bage %>%
ggplot(aes(x=grimAge, y=bmi)) +
geom_smooth(method='lm', formula= y~x, se=FALSE, color='#66ccEE', size=1.5) +
geom_point(fill=ifelse(targets_bage$timepoint==0,'#ee6677', '#66CCEE'), color='black', size=2, shape=21, stroke=1) +
ggtitle('Associations between grimAge and BMI in GOTO') +
theme(
plot.title = element_text(hjust = 0.5)
)
Session Info
sessionInfo()## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Rocky Linux 8.10 (Green Obsidian)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.15.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] SummarizedExperiment_1.24.0 Biobase_2.58.0
## [3] GenomicRanges_1.46.1 GenomeInfoDb_1.34.9
## [5] IRanges_2.32.0 S4Vectors_0.36.2
## [7] BiocGenerics_0.44.0 MatrixGenerics_1.10.0
## [9] matrixStats_1.0.0 lmerTest_3.1-3
## [11] lme4_1.1-30 Matrix_1.5-4.1
## [13] forcats_0.5.2 stringr_1.5.0
## [15] dplyr_1.1.3 purrr_0.3.4
## [17] readr_2.1.2 tidyr_1.2.1
## [19] tibble_3.2.1 ggplot2_3.4.3
## [21] tidyverse_1.3.2 rmarkdown_2.16
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-162 bitops_1.0-7 fs_1.6.2
## [4] bit64_4.0.5 lubridate_1.9.2 httr_1.4.6
## [7] numDeriv_2016.8-1.1 tools_4.2.2 backports_1.4.1
## [10] bslib_0.5.0 utf8_1.2.3 R6_2.5.1
## [13] mgcv_1.8-42 DBI_1.1.3 colorspace_2.1-0
## [16] withr_2.5.0 tidyselect_1.2.0 bit_4.0.5
## [19] compiler_4.2.2 cli_3.6.1 rvest_1.0.3
## [22] xml2_1.3.4 DelayedArray_0.24.0 labeling_0.4.2
## [25] sass_0.4.6 scales_1.2.1 digest_0.6.31
## [28] minqa_1.2.5 XVector_0.34.0 pkgconfig_2.0.3
## [31] htmltools_0.5.5 highr_0.10 dbplyr_2.2.1
## [34] fastmap_1.1.1 rlang_1.1.1 readxl_1.4.1
## [37] rstudioapi_0.14 farver_2.1.1 jquerylib_0.1.4
## [40] generics_0.1.3 jsonlite_1.8.5 vroom_1.5.7
## [43] googlesheets4_1.0.1 RCurl_1.98-1.12 magrittr_2.0.3
## [46] GenomeInfoDbData_1.2.9 Rcpp_1.0.10 munsell_0.5.0
## [49] fansi_1.0.4 lifecycle_1.0.3 stringi_1.7.12
## [52] yaml_2.3.7 MASS_7.3-60 zlibbioc_1.44.0
## [55] grid_4.2.2 parallel_4.2.2 crayon_1.5.2
## [58] lattice_0.21-8 haven_2.5.1 splines_4.2.2
## [61] hms_1.1.2 knitr_1.43 pillar_1.9.0
## [64] boot_1.3-28.1 reprex_2.0.2 glue_1.6.2
## [67] evaluate_0.21 modelr_0.1.9 vctrs_0.6.3
## [70] nloptr_2.0.3 tzdb_0.4.0 cellranger_1.1.0
## [73] gtable_0.3.3 assertthat_0.2.1 cachem_1.0.8
## [76] xfun_0.39 broom_1.0.1 googledrive_2.0.0
## [79] gargle_1.5.0 timechange_0.2.0 ellipsis_0.3.2Clear
rm(list=ls())