This script calculates all three cAge clocks and one bAge clock using the dnaMethyAge package


Setup

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')

Subset

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

Calculate

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.


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] 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())