Script 22: DNAm-trait associations in muscle
Compiled: June 26, 2024
This script assesses if changes in metabolic and muscle parameters change with DNAm at muscle CpGs.
Setup
Load packages
library(edgeR)
library(cinaR)
library(limma)
library(tidyverse)
library(GenomicFeatures)
library(AnnotationHub)
library(SummarizedExperiment)
library(lme4)
library(lmerTest)
library(pheatmap)
library(RColorBrewer)
library(ComplexHeatmap)
library(circlize)DNAm data
Load top CpG results and methData
load('../GOTO_Data/GOTO_methData-filtered.Rdata')
load('../GOTO_Data/GOTO_results-top-muscle-adj.Rdata')Save only muscle samples
methData <- methData[ , methData$tissue == 'muscle']Convert betas to data frame
beta_rows <- methData$Basename
betas <- as.data.frame(t(assay(methData)))
rownames(betas) <- beta_rowsSave only CpGs to be tested
betas <- betas %>% dplyr::select(any_of(top_cpgs$cpg))
dim(betas)## [1] 160 162Phenotype data
Load data
load("../GOTO_Data/GOTO_targets-filtered.Rdata")Filter muscle
targets <- targets %>%
filter(tissue == 'muscle')Save covariates
targets <- targets %>%
dplyr::select(IOP2_ID, timepoint, age, sex, smoke,
plate, array_row, Basename,
bmi, wc, perc_body_fat, fasting_insulin,
leptin, adiponectin, interleukin_6,
hdl_size_fasted, hdl_c, systolic_bp) %>%
mutate(ID = paste0(IOP2_ID, '_', timepoint))Ordering and filtering
Ensure DNAm data for all samples
targets <- targets %>% filter(Basename %in% rownames(betas))
print(paste0('We have complete covariate and DNAm data for ',
nrow(targets), ' samples'))## [1] "We have complete covariate and DNAm data for 160 samples"Order uuid alphabetically
targets <- targets[order(targets$Basename),]
rownames(targets) <- targets$Basename
dim(targets)## [1] 160 19Order DNAm data
betas <- betas[rownames(betas) %in% targets$Basename, ]
betas <- betas[order(rownames(betas)), ]
dim(betas)## [1] 160 162Analysis
Bind data frames
lm_df <- cbind(targets,
betas)Save traits
traits <- c('bmi',
'wc',
'perc_body_fat',
'fasting_insulin',
'leptin',
'adiponectin',
'interleukin_6',
'hdl_size_fasted',
'hdl_c',
'systolic_bp')Save CpGs
cpgs <- top_cpgs$cpg
length(cpgs)## [1] 162Expand pairwise comparisons
df <- expand.grid(traits, cpgs)
colnames(df) <- c('trait', 'cpg')Run models
for(i in 1:nrow(df)){
cpg <- as.character(df$cpg[i])
trait <- as.character(df$trait[i])
fit <- lmer(substitute(cpg ~ trait + age + sex + smoke +
plate + array_row + (1 | IOP2_ID),
list(cpg = as.name(cpg),
trait = as.name(trait))),
data=lm_df)
if(i == 1){
out <- data.frame(
cpg = cpg,
trait = trait,
es = coef(summary(fit))[2,1],
se = coef(summary(fit))[2,2],
p = coef(summary(fit))[2,5]
)
} else {
out <- rbind(out,
data.frame(
cpg = cpg,
trait = trait,
es = coef(summary(fit))[2,1],
se = coef(summary(fit))[2,2],
p = coef(summary(fit))[2,5]))
}
}Adjust p-values
out <- out %>% group_by(cpg) %>%
mutate(padj = p.adjust(p, method='fdr')) %>%
ungroup()Inspect top
out <- out %>% arrange(padj)
head(out)## # A tibble: 6 × 6
## cpg trait es se p padj
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cg27187848 wc 0.00227 0.000498 0.0000115 0.000115
## 2 cg16350675 leptin -0.00247 0.000575 0.0000412 0.000412
## 3 cg17411016 wc 0.00224 0.000567 0.000125 0.00125
## 4 cg25114611 leptin -0.00207 0.000534 0.000198 0.00198
## 5 cg13733654 wc -0.000701 0.000188 0.000290 0.00290
## 6 cg11835462 leptin -0.00226 0.000607 0.000330 0.00330Save
save(out,
file='../GOTO_Data/Traits/GOTO-muscle_traitMetabolic.Rdata')
write_csv(out, "../GOTO_Data/Tables/ST08.csv")Muscle based traits
Save previous
blood <- outAverage dominant grip
load("../GOTO_Data/Traits/avg_dom.Rdata")
grip <- df %>%
mutate(ID = paste0(
IOP2_ID, "_", ifelse(df$timepoint==T, 1, 0)
)) %>%
dplyr::select(ID, avg_dom = var)Merge
targets <- left_join(targets, grip, by="ID")PAX7
load("../GOTO_Data/Traits/pax7_per_fiber.Rdata")
pax7 <- df %>%
mutate(ID = paste0(
IOP2_ID, "_", ifelse(df$timepoint==T, 1, 0)
)) %>%
dplyr::select(ID, pax7 = var)Merge
targets <- left_join(targets, pax7, by="ID")Myonuclei per fiber
load("../GOTO_Data/Traits/myonuclei_per_fiber.Rdata")
myo <- df %>%
mutate(ID = paste0(
IOP2_ID, "_", ifelse(df$timepoint==T, 1, 0)
)) %>%
dplyr::select(ID, myo = var)Merge
targets <- left_join(targets, myo, by="ID")
targets <- targets %>% dplyr::select(ID, avg_dom, pax7, myo)
lm_df <- left_join(lm_df, targets, by="ID")Save traits
traits <- c('avg_dom',
'pax7',
'myo')Expand pairs
df <- expand.grid(traits, cpgs)
colnames(df) <- c('trait', 'cpg')Run models
for(i in 1:nrow(df)){
cpg <- as.character(df$cpg[i])
trait <- as.character(df$trait[i])
fit <- lmer(substitute(cpg ~ trait + age + sex + smoke +
plate + array_row + (1 | IOP2_ID),
list(cpg = as.name(cpg),
trait = as.name(trait))),
data=lm_df)
if(i == 1){
out <- data.frame(
cpg = cpg,
trait = trait,
es = coef(summary(fit))[2,1],
se = coef(summary(fit))[2,2],
p = coef(summary(fit))[2,5]
)
} else {
out <- rbind(out,
data.frame(
cpg = cpg,
trait = trait,
es = coef(summary(fit))[2,1],
se = coef(summary(fit))[2,2],
p = coef(summary(fit))[2,5]))
}
}Adjust p-values
out <- out %>% group_by(cpg) %>%
mutate(padj = p.adjust(p, method='fdr')) %>%
ungroup()Print top
out <- out %>% arrange(padj)
head(out)## # A tibble: 6 × 6
## cpg trait es se p padj
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cg09312464 pax7 -0.157 0.0239 0.00000000162 0.00000000485
## 2 cg15982707 pax7 -0.259 0.0397 0.00000000244 0.00000000733
## 3 cg24161080 pax7 -0.475 0.0743 0.00000000434 0.0000000130
## 4 cg17633300 pax7 -0.115 0.0182 0.00000000711 0.0000000213
## 5 cg11012616 pax7 -0.210 0.0333 0.00000000942 0.0000000283
## 6 cg17411016 pax7 -0.191 0.0304 0.0000000100 0.0000000301Save
save(out, file='../GOTO_Data/Traits/GOTO-muscle_traitMuscle.Rdata')
write_csv(out, "../GOTO_Data/Tables/ST09.csv")Heatmap
Format data
out <- rbind(blood, out)
out <- out %>%
mutate(trait = case_when(
trait == 'pax7' ~ "PAX7+ cells / fibre",
trait == "myo" ~ "Myonuclei / fibre",
trait == "avg_dom" ~ "Grip strength",
trait == "wc" ~ "WC",
trait == "perc_body_fat" ~ "Total Body Fat %",
trait == "leptin" ~ "Leptin",
trait == "hdl_size_fasted" ~ "HDL size",
trait == "bmi" ~ "BMI",
trait == "adiponectin" ~ "Adiponectin",
trait == "interleukin_6" ~ "Interleukin-6",
trait == "systolic_bp" ~ "SBP",
trait == "fasting_insulin" ~ "Fasting insulin",
trait == "hdl_c" ~ "HDL-C"))
unique(out$trait)## [1] "WC" "Leptin" "HDL size"
## [4] "Total Body Fat %" "Adiponectin" "Interleukin-6"
## [7] "BMI" "Fasting insulin" "SBP"
## [10] "HDL-C" "PAX7+ cells / fibre" "Myonuclei / fibre"
## [13] "Grip strength"Make df
heat_df <- reshape2::dcast(out, trait ~ cpg,
value.var = "padj")
heat_df <- heat_df %>% column_to_rownames(var='trait')
heat_df <- as.matrix(heat_df)Draw heatmap
hm <- draw(Heatmap(
heat_df,
column_title_gp = gpar(fontsize = 20, fontface = "bold"),
name = "padj",
col = colorRamp2(c(0, 0.001, 0.05, 0.1, 0.2, 1), c("#ab1226", "#e8213b", "#ef6678", "#b07f8b", "grey85", "grey95")),
show_row_dend = FALSE,
show_column_names = FALSE,
column_dend_height = unit(1, "cm"),
column_dend_side = "top",
width = ncol(heat_df)*unit(0.7, "mm"),
height = nrow(heat_df)*unit(3, "mm"),
column_km = 1, column_gap = unit(2, "mm"), border = TRUE,
row_names_gp = gpar(fontsize = 8),
heatmap_legend_param = list(
legend_direction = "horizontal",
legend_width = unit(3, "cm"),
border = 'black')),
heatmap_legend_side = "bottom")
Plot
print(hm)
Save
png("../GOTO_Data/Figures/Figure_2A.png")
print(hm)
dev.off()## png
## 2Session 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] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] circlize_0.4.15
## [2] ComplexHeatmap_2.14.0
## [3] RColorBrewer_1.1-3
## [4] pheatmap_1.0.12
## [5] clusterProfiler_4.2.2
## [6] AnnotationHub_3.2.2
## [7] BiocFileCache_2.2.1
## [8] dbplyr_2.2.1
## [9] cinaR_0.2.3
## [10] edgeR_3.40.2
## [11] ggpubr_0.4.0
## [12] GEOquery_2.62.2
## [13] MuSiC_0.2.0
## [14] nnls_1.4
## [15] gplots_3.1.3
## [16] plotly_4.10.1
## [17] SeuratObject_4.1.3
## [18] Seurat_4.3.0
## [19] gridExtra_2.3
## [20] lattice_0.21-8
## [21] bacon_1.22.0
## [22] ellipse_0.4.5
## [23] methylGSA_1.12.0
## [24] sva_3.42.0
## [25] genefilter_1.76.0
## [26] mgcv_1.8-42
## [27] nlme_3.1-162
## [28] limma_3.54.2
## [29] lmerTest_3.1-3
## [30] lme4_1.1-30
## [31] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [32] snpStats_1.44.0
## [33] survival_3.5-5
## [34] ggrepel_0.9.1
## [35] ggfortify_0.4.14
## [36] irlba_2.3.5.1
## [37] Matrix_1.5-4.1
## [38] omicsPrint_1.14.0
## [39] MASS_7.3-60
## [40] DNAmArray_2.0.0
## [41] pls_2.8-2
## [42] FDb.InfiniumMethylation.hg19_2.2.0
## [43] org.Hs.eg.db_3.14.0
## [44] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [45] GenomicFeatures_1.46.5
## [46] AnnotationDbi_1.56.2
## [47] IlluminaHumanMethylationEPICmanifest_0.3.0
## [48] minfi_1.40.0
## [49] bumphunter_1.36.0
## [50] locfit_1.5-9.8
## [51] iterators_1.0.14
## [52] foreach_1.5.2
## [53] Biostrings_2.62.0
## [54] XVector_0.34.0
## [55] SummarizedExperiment_1.24.0
## [56] Biobase_2.58.0
## [57] MatrixGenerics_1.10.0
## [58] matrixStats_1.0.0
## [59] GenomicRanges_1.46.1
## [60] GenomeInfoDb_1.34.9
## [61] IRanges_2.32.0
## [62] S4Vectors_0.36.2
## [63] BiocGenerics_0.44.0
## [64] BiocParallel_1.32.6
## [65] MethylAid_1.28.0
## [66] forcats_0.5.2
## [67] stringr_1.5.0
## [68] dplyr_1.1.3
## [69] purrr_0.3.4
## [70] readr_2.1.2
## [71] tidyr_1.2.1
## [72] tibble_3.2.1
## [73] ggplot2_3.4.3
## [74] tidyverse_1.3.2
## [75] rmarkdown_2.16
##
## loaded via a namespace (and not attached):
## [1] graphlayouts_0.8.1
## [2] pbapply_1.7-0
## [3] haven_2.5.1
## [4] vctrs_0.6.3
## [5] beanplot_1.3.1
## [6] blob_1.2.4
## [7] spatstat.data_3.0-1
## [8] later_1.3.1
## [9] nloptr_2.0.3
## [10] DBI_1.1.3
## [11] rappdirs_0.3.3
## [12] uwot_0.1.14
## [13] zlibbioc_1.44.0
## [14] MatrixModels_0.5-1
## [15] GlobalOptions_0.1.2
## [16] htmlwidgets_1.5.4
## [17] future_1.32.0
## [18] leiden_0.4.3
## [19] illuminaio_0.40.0
## [20] tidygraph_1.2.2
## [21] Rcpp_1.0.10
## [22] KernSmooth_2.23-21
## [23] promises_1.2.0.1
## [24] DelayedArray_0.24.0
## [25] magick_2.7.4
## [26] fs_1.6.2
## [27] fastmatch_1.1-3
## [28] digest_0.6.31
## [29] png_0.1-8
## [30] nor1mix_1.3-0
## [31] sctransform_0.3.5
## [32] scatterpie_0.1.8
## [33] cowplot_1.1.1
## [34] DOSE_3.20.1
## [35] ggraph_2.0.6
## [36] pkgconfig_2.0.3
## [37] GO.db_3.14.0
## [38] gridBase_0.4-7
## [39] spatstat.random_3.1-5
## [40] DelayedMatrixStats_1.16.0
## [41] minqa_1.2.5
## [42] reticulate_1.30
## [43] GetoptLong_1.0.5
## [44] xfun_0.39
## [45] bslib_0.5.0
## [46] zoo_1.8-12
## [47] tidyselect_1.2.0
## [48] reshape2_1.4.4
## [49] ica_1.0-3
## [50] viridisLite_0.4.2
## [51] rtracklayer_1.54.0
## [52] rlang_1.1.1
## [53] hexbin_1.28.3
## [54] jquerylib_0.1.4
## [55] glue_1.6.2
## [56] modelr_0.1.9
## [57] ggsignif_0.6.3
## [58] labeling_0.4.2
## [59] SparseM_1.81
## [60] httpuv_1.6.11
## [61] preprocessCore_1.60.2
## [62] reactome.db_1.77.0
## [63] DO.db_2.9
## [64] annotate_1.72.0
## [65] jsonlite_1.8.5
## [66] bit_4.0.5
## [67] mime_0.12
## [68] Rsamtools_2.10.0
## [69] stringi_1.7.12
## [70] spatstat.sparse_3.0-1
## [71] scattermore_0.8
## [72] spatstat.explore_3.1-0
## [73] yulab.utils_0.0.6
## [74] quadprog_1.5-8
## [75] bitops_1.0-7
## [76] cli_3.6.1
## [77] rhdf5filters_1.10.1
## [78] RSQLite_2.2.17
## [79] data.table_1.14.8
## [80] timechange_0.2.0
## [81] rstudioapi_0.14
## [82] GenomicAlignments_1.30.0
## [83] qvalue_2.26.0
## [84] listenv_0.9.0
## [85] miniUI_0.1.1.1
## [86] gridGraphics_0.5-1
## [87] readxl_1.4.1
## [88] lifecycle_1.0.3
## [89] htm2txt_2.2.2
## [90] munsell_0.5.0
## [91] cellranger_1.1.0
## [92] caTools_1.18.2
## [93] codetools_0.2-19
## [94] coda_0.19-4
## [95] MultiAssayExperiment_1.20.0
## [96] lmtest_0.9-40
## [97] missMethyl_1.28.0
## [98] xtable_1.8-4
## [99] ROCR_1.0-11
## [100] googlesheets4_1.0.1
## [101] BiocManager_1.30.21
## [102] abind_1.4-5
## [103] farver_2.1.1
## [104] parallelly_1.36.0
## [105] RANN_2.6.1
## [106] aplot_0.1.7
## [107] askpass_1.1
## [108] ggtree_3.2.1
## [109] BiocIO_1.8.0
## [110] RcppAnnoy_0.0.20
## [111] goftest_1.2-3
## [112] patchwork_1.1.2
## [113] cluster_2.1.4
## [114] future.apply_1.11.0
## [115] tidytree_0.4.0
## [116] ellipsis_0.3.2
## [117] prettyunits_1.1.1
## [118] lubridate_1.9.2
## [119] ggridges_0.5.4
## [120] googledrive_2.0.0
## [121] reprex_2.0.2
## [122] mclust_6.0.0
## [123] igraph_1.4.3
## [124] multtest_2.50.0
## [125] fgsea_1.20.0
## [126] gargle_1.5.0
## [127] spatstat.utils_3.0-3
## [128] htmltools_0.5.5
## [129] yaml_2.3.7
## [130] utf8_1.2.3
## [131] MCMCpack_1.6-3
## [132] interactiveDisplayBase_1.32.0
## [133] XML_3.99-0.14
## [134] withr_2.5.0
## [135] fitdistrplus_1.1-11
## [136] bit64_4.0.5
## [137] rngtools_1.5.2
## [138] doRNG_1.8.6
## [139] progressr_0.13.0
## [140] GOSemSim_2.20.0
## [141] memoise_2.0.1
## [142] evaluate_0.21
## [143] tzdb_0.4.0
## [144] curl_5.0.1
## [145] fansi_1.0.4
## [146] highr_0.10
## [147] tensor_1.5
## [148] cachem_1.0.8
## [149] deldir_1.0-9
## [150] rjson_0.2.21
## [151] rstatix_0.7.0
## [152] clue_0.3-64
## [153] tools_4.2.2
## [154] sass_0.4.6
## [155] magrittr_2.0.3
## [156] RCurl_1.98-1.12
## [157] car_3.1-0
## [158] ape_5.7-1
## [159] ggplotify_0.1.0
## [160] xml2_1.3.4
## [161] httr_1.4.6
## [162] assertthat_0.2.1
## [163] boot_1.3-28.1
## [164] globals_0.16.2
## [165] R6_2.5.1
## [166] Rhdf5lib_1.20.0
## [167] progress_1.2.2
## [168] KEGGREST_1.34.0
## [169] treeio_1.18.1
## [170] shape_1.4.6
## [171] gtools_3.9.4
## [172] statmod_1.5.0
## [173] BiocVersion_3.16.0
## [174] HDF5Array_1.22.1
## [175] rhdf5_2.42.1
## [176] splines_4.2.2
## [177] carData_3.0-5
## [178] ggfun_0.0.7
## [179] colorspace_2.1-0
## [180] generics_0.1.3
## [181] RobustRankAggreg_1.2.1
## [182] pillar_1.9.0
## [183] tweenr_2.0.2
## [184] sp_1.6-1
## [185] GenomeInfoDbData_1.2.9
## [186] plyr_1.8.8
## [187] gtable_0.3.3
## [188] rvest_1.0.3
## [189] restfulr_0.0.15
## [190] knitr_1.43
## [191] shadowtext_0.1.2
## [192] biomaRt_2.50.3
## [193] fastmap_1.1.1
## [194] Cairo_1.6-0
## [195] doParallel_1.0.17
## [196] quantreg_5.94
## [197] broom_1.0.1
## [198] openssl_2.0.6
## [199] scales_1.2.1
## [200] filelock_1.0.2
## [201] backports_1.4.1
## [202] RaggedExperiment_1.18.0
## [203] base64_2.0.1
## [204] vroom_1.5.7
## [205] enrichplot_1.14.2
## [206] mcmc_0.9-7
## [207] hms_1.1.2
## [208] ggforce_0.3.4
## [209] scrime_1.3.5
## [210] Rtsne_0.16
## [211] shiny_1.7.2
## [212] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.0
## [213] polyclip_1.10-4
## [214] numDeriv_2016.8-1.1
## [215] siggenes_1.68.0
## [216] lazyeval_0.2.2
## [217] crayon_1.5.2
## [218] downloader_0.4
## [219] sparseMatrixStats_1.10.0
## [220] viridis_0.6.2
## [221] reshape_0.8.9
## [222] compiler_4.2.2
## [223] spatstat.geom_3.2-1Clear
rm(list=ls())