This script calculates all three cAge clocks and one bAge clock using the dnaMethyAge package
Load required packages
library(tidyverse)
library(dnaMethyAge)
library(SummarizedExperiment)
Load data
load('../GOTO_Data/GOTO_targets-filtered.Rdata')
load('../GOTO_Data/Processing/GOTO_methData-unfiltered.Rdata')
Blood
methData <- methData_unfiltered[ , methData_unfiltered$tissue == 'fasted blood']
targets <- targets %>% filter(tissue == "fasted blood")
betas <- assay(methData)
methData
## class: SummarizedExperiment
## dim: 865859 196
## metadata(0):
## assays(1): beta
## rownames(865859): cg14817997 cg26928153 ... cg07587934 cg16855331
## rowData names(4): cpg chr start end
## colnames(196): 203527980080_R01C01 203527980080_R02C01 ...
## 203550300093_R07C01 203550300093_R08C01
## colData names(19): DNA_labnr IOP2_ID ... Basename smoke
dim(targets)
## [1] 196 102
Select clocks
availableClock()
## To understand what do these clocks represent for, please refer to 'https://github.com/yiluyucheng/dnaMethyAge'.
## [1] "HannumG2013" "HorvathS2013" "LevineM2018" "ZhangQ2019"
## [5] "ShirebyG2020" "YangZ2016" "ZhangY2017" "LuA2019"
## [9] "HorvathS2018" "DunedinPACE" "McEwenL2019" "CBL_specific"
## [13] "PCGrimAge" "PCHorvathS2013" "PCHannumG2013" "PCHorvathS2018"
## [17] "PCPhenoAge" "CBL_common" "Cortex_common" "epiTOC2"
## [21] "BernabeuE2023c"
av <- c("HorvathS2013", #cAge Horvath
"LevineM2018", #bAge PhenoAge
"ZhangQ2019", #cAge Zhang
"BernabeuE2023c") #cAge Bernabeu
info <- targets %>%
select(Sample=Basename, Age=age, Sex=sex) %>%
mutate(Sex = ifelse(Sex == 'male', 'Male', 'Female'))
dim(info)
## [1] 196 3
for(i in av){
age <- methyAge(betas,
clock=i,
age_info=info,
fit_method='Linear',
do_plot=TRUE)
age$Basename <- age$Sample
age <- full_join(age, targets,
by='Basename')
save(age,
file=paste0('../GOTO_Data/Clocks/Clock_Out/',
i,'_blood.Rdata'))
}
## Start data imputation ...
## Loading required package: impute
## Cluster size 20306 broken into 14924 5382
## Cluster size 14924 broken into 13163 1761
## Cluster size 13163 broken into 2379 10784
## Cluster size 2379 broken into 1477 902
## Done cluster 1477
## Done cluster 902
## Done cluster 2379
## Cluster size 10784 broken into 3908 6876
## Cluster size 3908 broken into 1594 2314
## Cluster size 1594 broken into 447 1147
## Done cluster 447
## Done cluster 1147
## Done cluster 1594
## Cluster size 2314 broken into 1154 1160
## Done cluster 1154
## Done cluster 1160
## Done cluster 2314
## Done cluster 3908
## Cluster size 6876 broken into 4011 2865
## Cluster size 4011 broken into 1755 2256
## Cluster size 1755 broken into 855 900
## Done cluster 855
## Done cluster 900
## Done cluster 1755
## Cluster size 2256 broken into 2228 28
## Cluster size 2228 broken into 2171 57
## Cluster size 2171 broken into 1127 1044
## Done cluster 1127
## Done cluster 1044
## Done cluster 2171
## Done cluster 57
## Done cluster 2228
## Done cluster 28
## Done cluster 2256
## Done cluster 4011
## Cluster size 2865 broken into 1914 951
## Cluster size 1914 broken into 1900 14
## Cluster size 1900 broken into 930 970
## Done cluster 930
## Done cluster 970
## Done cluster 1900
## Done cluster 14
## Done cluster 1914
## Done cluster 951
## Done cluster 2865
## Done cluster 6876
## Done cluster 10784
## Done cluster 13163
## Cluster size 1761 broken into 784 977
## Done cluster 784
## Done cluster 977
## Done cluster 1761
## Done cluster 14924
## Cluster size 5382 broken into 3772 1610
## Cluster size 3772 broken into 2479 1293
## Cluster size 2479 broken into 1378 1101
## Done cluster 1378
## Done cluster 1101
## Done cluster 2479
## Done cluster 1293
## Done cluster 3772
## Cluster size 1610 broken into 721 889
## Done cluster 721
## Done cluster 889
## Done cluster 1610
## Done cluster 5382
## Finished imputation.
## Normalization by adusted BMIQ ...
## Estimate running time for normalisation (BMIQ with fixed reference) by a single core: 24.2 minutes.
## Age acceleration is calculated as the residual resulting from a linear regression model which DNAm age is regressed on chronological age.
## Age acceleration is calculated as the residual resulting from a linear regression model which DNAm age is regressed on chronological age.
## Zscore Standardizing:
## Age acceleration is calculated as the residual resulting from a linear regression model which DNAm age is regressed on chronological age.
## Age acceleration is calculated as the residual resulting from a linear regression model which DNAm age is regressed on chronological age.
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] RPMM_1.25 cluster_2.1.4
## [3] impute_1.72.3 SummarizedExperiment_1.24.0
## [5] Biobase_2.58.0 GenomicRanges_1.46.1
## [7] GenomeInfoDb_1.34.9 IRanges_2.32.0
## [9] S4Vectors_0.36.2 BiocGenerics_0.44.0
## [11] MatrixGenerics_1.10.0 matrixStats_1.0.0
## [13] dnaMethyAge_0.1.0 forcats_0.5.2
## [15] stringr_1.5.0 dplyr_1.1.3
## [17] purrr_0.3.4 readr_2.1.2
## [19] tidyr_1.2.1 tibble_3.2.1
## [21] ggplot2_3.4.3 tidyverse_1.3.2
## [23] rmarkdown_2.16
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.6 sass_0.4.6 jsonlite_1.8.5
## [4] modelr_0.1.9 bslib_0.5.0 assertthat_0.2.1
## [7] highr_0.10 GenomeInfoDbData_1.2.9 googlesheets4_1.0.1
## [10] cellranger_1.1.0 yaml_2.3.7 lattice_0.21-8
## [13] pillar_1.9.0 backports_1.4.1 glue_1.6.2
## [16] digest_0.6.31 XVector_0.34.0 rvest_1.0.3
## [19] colorspace_2.1-0 htmltools_0.5.5 Matrix_1.5-4.1
## [22] pkgconfig_2.0.3 broom_1.0.1 haven_2.5.1
## [25] zlibbioc_1.44.0 scales_1.2.1 tzdb_0.4.0
## [28] timechange_0.2.0 googledrive_2.0.0 generics_0.1.3
## [31] ellipsis_0.3.2 cachem_1.0.8 withr_2.5.0
## [34] cli_3.6.1 magrittr_2.0.3 crayon_1.5.2
## [37] readxl_1.4.1 evaluate_0.21 fs_1.6.2
## [40] fansi_1.0.4 xml2_1.3.4 tools_4.2.2
## [43] hms_1.1.2 gargle_1.5.0 lifecycle_1.0.3
## [46] munsell_0.5.0 reprex_2.0.2 DelayedArray_0.24.0
## [49] compiler_4.2.2 jquerylib_0.1.4 rlang_1.1.1
## [52] grid_4.2.2 RCurl_1.98-1.12 rstudioapi_0.14
## [55] bitops_1.0-7 gtable_0.3.3 DBI_1.1.3
## [58] R6_2.5.1 lubridate_1.9.2 knitr_1.43
## [61] fastmap_1.1.1 utf8_1.2.3 stringi_1.7.12
## [64] parallel_4.2.2 vctrs_0.6.3 dbplyr_2.2.1
## [67] tidyselect_1.2.0 xfun_0.39
Clear
rm(list=ls())