Script 31: DNAm-trait associations in adipose
Compiled: June 26, 2024
This script assesses if health improvements are associated with differential methylation in adipose tissue.
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-fat.Rdata')Save only adipose samples
methData <- methData[ , methData$tissue == 'fat']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] 178 230Phenotype data
Load data
load("../GOTO_Data/GOTO_targets-filtered.Rdata")Filter adipose
targets <- targets %>%
filter(tissue == 'fat')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 178 samples"Order uuid alphabetically
targets <- targets[order(targets$Basename),]
rownames(targets) <- targets$Basename
dim(targets)## [1] 178 19Order DNAm data
betas <- betas[rownames(betas) %in% targets$Basename, ]
betas <- betas[order(rownames(betas)), ]
dim(betas)## [1] 178 230Analysis
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] 230Expand 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 cg07219853 interleukin_6 -0.0211 0.00507 0.0000506 0.000506
## 2 cg04715165 interleukin_6 -0.0211 0.00518 0.0000734 0.000734
## 3 cg04543233 perc_body_fat -0.00454 0.00113 0.0000896 0.000896
## 4 cg06161697 perc_body_fat 0.00323 0.000805 0.0000922 0.000922
## 5 cg03160217 interleukin_6 -0.0120 0.00302 0.000102 0.00102
## 6 cg20907471 perc_body_fat -0.00616 0.00162 0.000202 0.00202Save
save(out,
file='../GOTO_Data/Traits/GOTO-fat_trait.Rdata')
write_csv(out, "../GOTO_Data/Tables/ST15.csv")Heatmap
Format data
out <- out %>%
mutate(trait = case_when(
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] "Interleukin-6" "Total Body Fat %" "Fasting insulin" "Leptin"
## [5] "WC" "Adiponectin" "BMI" "HDL size"
## [9] "HDL-C" "SBP"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.5, "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")
Save
png("../GOTO_Data/Figures/Figure_4A.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 stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] circlize_0.4.15 ComplexHeatmap_2.14.0
## [3] RColorBrewer_1.1-3 pheatmap_1.0.12
## [5] lmerTest_3.1-3 lme4_1.1-30
## [7] Matrix_1.5-4.1 SummarizedExperiment_1.24.0
## [9] MatrixGenerics_1.10.0 matrixStats_1.0.0
## [11] AnnotationHub_3.2.2 BiocFileCache_2.2.1
## [13] dbplyr_2.2.1 GenomicFeatures_1.46.5
## [15] AnnotationDbi_1.56.2 Biobase_2.58.0
## [17] GenomicRanges_1.46.1 GenomeInfoDb_1.34.9
## [19] IRanges_2.32.0 S4Vectors_0.36.2
## [21] BiocGenerics_0.44.0 forcats_0.5.2
## [23] stringr_1.5.0 dplyr_1.1.3
## [25] purrr_0.3.4 readr_2.1.2
## [27] tidyr_1.2.1 tibble_3.2.1
## [29] ggplot2_3.4.3 tidyverse_1.3.2
## [31] cinaR_0.2.3 edgeR_3.40.2
## [33] limma_3.54.2 rmarkdown_2.16
##
## loaded via a namespace (and not attached):
## [1] readxl_1.4.1 backports_1.4.1
## [3] plyr_1.8.8 splines_4.2.2
## [5] BiocParallel_1.32.6 digest_0.6.31
## [7] foreach_1.5.2 htmltools_0.5.5
## [9] magick_2.7.4 fansi_1.0.4
## [11] magrittr_2.0.3 memoise_2.0.1
## [13] cluster_2.1.4 doParallel_1.0.17
## [15] googlesheets4_1.0.1 tzdb_0.4.0
## [17] Biostrings_2.62.0 modelr_0.1.9
## [19] vroom_1.5.7 timechange_0.2.0
## [21] prettyunits_1.1.1 colorspace_2.1-0
## [23] blob_1.2.4 rvest_1.0.3
## [25] rappdirs_0.3.3 haven_2.5.1
## [27] xfun_0.39 crayon_1.5.2
## [29] RCurl_1.98-1.12 jsonlite_1.8.5
## [31] iterators_1.0.14 glue_1.6.2
## [33] gtable_0.3.3 gargle_1.5.0
## [35] zlibbioc_1.44.0 XVector_0.34.0
## [37] GetoptLong_1.0.5 DelayedArray_0.24.0
## [39] shape_1.4.6 scales_1.2.1
## [41] DBI_1.1.3 Rcpp_1.0.10
## [43] xtable_1.8-4 progress_1.2.2
## [45] clue_0.3-64 bit_4.0.5
## [47] httr_1.4.6 ellipsis_0.3.2
## [49] pkgconfig_2.0.3 XML_3.99-0.14
## [51] sass_0.4.6 locfit_1.5-9.8
## [53] utf8_1.2.3 reshape2_1.4.4
## [55] tidyselect_1.2.0 rlang_1.1.1
## [57] later_1.3.1 munsell_0.5.0
## [59] BiocVersion_3.16.0 cellranger_1.1.0
## [61] tools_4.2.2 cachem_1.0.8
## [63] cli_3.6.1 generics_0.1.3
## [65] RSQLite_2.2.17 broom_1.0.1
## [67] evaluate_0.21 fastmap_1.1.1
## [69] yaml_2.3.7 knitr_1.43
## [71] bit64_4.0.5 fs_1.6.2
## [73] KEGGREST_1.34.0 nlme_3.1-162
## [75] mime_0.12 xml2_1.3.4
## [77] biomaRt_2.50.3 compiler_4.2.2
## [79] rstudioapi_0.14 filelock_1.0.2
## [81] curl_5.0.1 png_0.1-8
## [83] interactiveDisplayBase_1.32.0 reprex_2.0.2
## [85] bslib_0.5.0 stringi_1.7.12
## [87] highr_0.10 lattice_0.21-8
## [89] nloptr_2.0.3 vctrs_0.6.3
## [91] pillar_1.9.0 lifecycle_1.0.3
## [93] BiocManager_1.30.21 GlobalOptions_0.1.2
## [95] jquerylib_0.1.4 bitops_1.0-7
## [97] httpuv_1.6.11 rtracklayer_1.54.0
## [99] R6_2.5.1 BiocIO_1.8.0
## [101] promises_1.2.0.1 codetools_0.2-19
## [103] boot_1.3-28.1 MASS_7.3-60
## [105] assertthat_0.2.1 rjson_0.2.21
## [107] withr_2.5.0 GenomicAlignments_1.30.0
## [109] Rsamtools_2.10.0 GenomeInfoDbData_1.2.9
## [111] parallel_4.2.2 hms_1.1.2
## [113] minqa_1.2.5 googledrive_2.0.0
## [115] Cairo_1.6-0 numDeriv_2016.8-1.1
## [117] shiny_1.7.2 lubridate_1.9.2
## [119] restfulr_0.0.15Clear
rm(list=ls())