This script saves dot plots of the selected CpG, gene, and trait in muscle for Figure 4B.


Setup

Load packages

library(tidyverse)
library(ggrepel)
library(GenomicRanges)
library(ggpubr)
library(DNAmArray)
library(SummarizedExperiment)

CpG

Load top CpG results and methData

load('../GOTO_Data/GOTO_methData-filtered.Rdata')
load('../GOTO_Data/GOTO_results-top-muscle-adj.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_rows

Save selected CpGs

betas <- betas %>% dplyr::select(all_of(c("cg02649849",
                                          "cg14434922",
                                          "cg14504259")))
dim(betas)
## [1] 178   3
cpg_list <- colnames(betas)

Phenotype 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))
ID_list <- targets %>% dplyr::select(ID)

CpG figure

Bind data frames

limma_base <- cbind(targets,
               betas)
dim(limma_base)
## [1] 178  22

cg14504259

Make box plot

plot <- limma_base %>% 
  ggplot(aes(x = timepoint, y = cg14504259,
             fill = timepoint)) +
  geom_line(aes(group = IOP2_ID), color = '#DDDDDD', linewidth=1) +
  geom_dotplot(binaxis = "y", binwidth = 0.02,
               stackdir = "center",
               dotsize = 0.9, stroke = 2) +
  scale_fill_manual(values=c('#EE6677', '#66CCEE'))+
  ggtitle('') +
  ylab('') + xlab('') +
  scale_x_discrete(expand = c(0,0)) +
  theme(
    panel.background = element_rect(fill = 'white',
                                    color = 'black'),
    legend.position = 'none',
    plot.title =element_text(hjust=0.5),
    text=element_text(size=16))
print(plot)

Save

png('../GOTO_Data/Figures/Figure_4B-cg14504259.png')
print(plot)
dev.off()
## png 
##   2

cg14434922

Make box plot

plot <- limma_base %>% 
  ggplot(aes(x = timepoint, y = cg14434922,
             fill = timepoint)) +
  geom_line(aes(group = IOP2_ID), color = '#DDDDDD', linewidth=1) +
  geom_dotplot(binaxis = "y", binwidth = 0.015, 
               stackdir = "center",
               dotsize = 0.8, stroke = 2) +
  scale_fill_manual(values=c('#EE6677', '#66CCEE'))+
  ggtitle('') +
  ylab('') + xlab('') +
  scale_x_discrete(expand = c(0,0)) +
  theme(
    panel.background = element_rect(fill = 'white',
                                    color = 'black'),
    legend.position = 'none',
    plot.title =element_text(hjust=0.5),
    text=element_text(size=16))
print(plot)

Save

png('../GOTO_Data/Figures/Figure_4B-cg14434922.png')
print(plot)
dev.off()
## png 
##   2

cg02649849

Make box plot

plot <- limma_base %>% 
  ggplot(aes(x = timepoint, y = cg02649849,
             fill = timepoint)) +
  geom_line(aes(group = IOP2_ID), color = '#DDDDDD', linewidth=1) +
  geom_dotplot(binaxis = "y", binwidth = 0.015, 
               stackdir = "center",
               dotsize = 0.8, stroke = 2) +
  scale_fill_manual(values=c('#EE6677', '#66CCEE'))+
  ggtitle('') +
  ylab('') + xlab('') +
  scale_x_discrete(expand = c(0,0)) +
  theme(
    panel.background = element_rect(fill = 'white',
                                    color = 'black'),
    legend.position = 'none',
    plot.title =element_text(hjust=0.5),
    text=element_text(size=16))
print(plot)

Save

png('../GOTO_Data/Figures/Figure_4B-cg02649849.png')
print(plot)
dev.off()
## png 
##   2

Gene figure

ID list

ID_list <- targets %>% mutate(
  ID = paste0(IOP2_ID, '_', timepoint)
) %>% dplyr::select(ID)

Load RNAseq data

load('../GOTO_Data/DEGs/expData_adipose.RData')

goto_exp <- expData %>%  
  dplyr::select(IOP2_ID, timepoint, final_plate,
                       starts_with('ENS')) %>% 
  mutate(
  ID = str_c(IOP2_ID, '_', timepoint)
)

exp_df <- left_join(ID_list, goto_exp, by = 'ID')

dim(exp_df)
## [1]   178 41655

Gene symbols

ens2gene <- cinaR::grch37
name <- "ENSG00000064218"
df <- exp_df %>% filter(!is.na(timepoint))
df$IOP2_ID <- as.factor(df$IOP2_ID)
df$timepoint <- as.factor(df$timepoint)

Dot plot

plot <-df %>% 
  ggplot(aes(x = timepoint, y = ENSG00000064218,
             fill = timepoint)) +
  geom_line(aes(group = IOP2_ID), color = '#DDDDDD', linewidth=1) +
  geom_dotplot(binaxis = "y", binwidth = 0.5,
               stackdir = "center",
               dotsize = 0.8, stroke = 2) +
  scale_fill_manual(values=c('#EE6677', '#66CCEE'))+
  ggtitle('') +
  ylab('') + xlab('') +
  scale_x_discrete(expand = c(0,0)) +
  theme(
    panel.background = element_rect(fill = 'white',
                                    color = 'black'),
    legend.position = 'none',
    plot.title =element_text(hjust=0.5),
    text=element_text(size=16))
print(plot)

Save

png('../GOTO_Data/Figures/Figure_4B-DMRT3.png')
print(plot)
dev.off()
## png 
##   2

Trait figure

plot <- targets %>% 
  ggplot(aes(x = timepoint, y = perc_body_fat,
             fill = timepoint)) +
  geom_line(aes(group = IOP2_ID), color = '#DDDDDD', linewidth=1) +
  geom_dotplot(binaxis = "y", binwidth = 2,
               stackdir = "center",
               dotsize = 0.8, stroke = 2) +
  scale_fill_manual(values=c('#EE6677', '#66CCEE'))+
  ggtitle('') +
  ylab('') + xlab('') +
  scale_x_discrete(expand = c(0,0)) +
  theme(
    panel.background = element_rect(fill = 'white',
                                    color = 'black'),
    legend.position = 'none',
    plot.title =element_text(hjust=0.5),
    text=element_text(size=16))
print(plot)

Save

png('../GOTO_Data/Figures/Figure_4B-bodyfat.png')
print(plot)
dev.off()
## png 
##   2

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] parallel  grid      stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] DNAmArray_2.0.0                        
##  [2] pls_2.8-2                              
##  [3] FDb.InfiniumMethylation.hg19_2.2.0     
##  [4] org.Hs.eg.db_3.14.0                    
##  [5] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
##  [6] minfi_1.40.0                           
##  [7] bumphunter_1.36.0                      
##  [8] locfit_1.5-9.8                         
##  [9] iterators_1.0.14                       
## [10] foreach_1.5.2                          
## [11] Biostrings_2.62.0                      
## [12] XVector_0.34.0                         
## [13] ggpubr_0.4.0                           
## [14] ggrepel_0.9.1                          
## [15] circlize_0.4.15                        
## [16] ComplexHeatmap_2.14.0                  
## [17] RColorBrewer_1.1-3                     
## [18] pheatmap_1.0.12                        
## [19] lmerTest_3.1-3                         
## [20] lme4_1.1-30                            
## [21] Matrix_1.5-4.1                         
## [22] SummarizedExperiment_1.24.0            
## [23] MatrixGenerics_1.10.0                  
## [24] matrixStats_1.0.0                      
## [25] AnnotationHub_3.2.2                    
## [26] BiocFileCache_2.2.1                    
## [27] dbplyr_2.2.1                           
## [28] GenomicFeatures_1.46.5                 
## [29] AnnotationDbi_1.56.2                   
## [30] Biobase_2.58.0                         
## [31] GenomicRanges_1.46.1                   
## [32] GenomeInfoDb_1.34.9                    
## [33] IRanges_2.32.0                         
## [34] S4Vectors_0.36.2                       
## [35] BiocGenerics_0.44.0                    
## [36] forcats_0.5.2                          
## [37] stringr_1.5.0                          
## [38] dplyr_1.1.3                            
## [39] purrr_0.3.4                            
## [40] readr_2.1.2                            
## [41] tidyr_1.2.1                            
## [42] tibble_3.2.1                           
## [43] ggplot2_3.4.3                          
## [44] tidyverse_1.3.2                        
## [45] cinaR_0.2.3                            
## [46] edgeR_3.40.2                           
## [47] limma_3.54.2                           
## [48] rmarkdown_2.16                         
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.3                    tidyselect_1.2.0             
##   [3] RSQLite_2.2.17                BiocParallel_1.32.6          
##   [5] munsell_0.5.0                 preprocessCore_1.60.2        
##   [7] codetools_0.2-19              withr_2.5.0                  
##   [9] colorspace_2.1-0              filelock_1.0.2               
##  [11] highr_0.10                    knitr_1.43                   
##  [13] rstudioapi_0.14               ggsignif_0.6.3               
##  [15] labeling_0.4.2                GenomeInfoDbData_1.2.9       
##  [17] farver_2.1.1                  bit64_4.0.5                  
##  [19] rhdf5_2.42.1                  vctrs_0.6.3                  
##  [21] generics_0.1.3                xfun_0.39                    
##  [23] timechange_0.2.0              R6_2.5.1                     
##  [25] doParallel_1.0.17             illuminaio_0.40.0            
##  [27] clue_0.3-64                   reshape_0.8.9                
##  [29] rhdf5filters_1.10.1           bitops_1.0-7                 
##  [31] cachem_1.0.8                  DelayedArray_0.24.0          
##  [33] assertthat_0.2.1              promises_1.2.0.1             
##  [35] BiocIO_1.8.0                  scales_1.2.1                 
##  [37] vroom_1.5.7                   googlesheets4_1.0.1          
##  [39] gtable_0.3.3                  Cairo_1.6-0                  
##  [41] rlang_1.1.1                   genefilter_1.76.0            
##  [43] GlobalOptions_0.1.2           splines_4.2.2                
##  [45] rtracklayer_1.54.0            rstatix_0.7.0                
##  [47] GEOquery_2.62.2               gargle_1.5.0                 
##  [49] htm2txt_2.2.2                 broom_1.0.1                  
##  [51] BiocManager_1.30.21           yaml_2.3.7                   
##  [53] reshape2_1.4.4                abind_1.4-5                  
##  [55] modelr_0.1.9                  backports_1.4.1              
##  [57] httpuv_1.6.11                 tools_4.2.2                  
##  [59] nor1mix_1.3-0                 ellipsis_0.3.2               
##  [61] jquerylib_0.1.4               siggenes_1.68.0              
##  [63] Rcpp_1.0.10                   plyr_1.8.8                   
##  [65] sparseMatrixStats_1.10.0      progress_1.2.2               
##  [67] zlibbioc_1.44.0               RCurl_1.98-1.12              
##  [69] prettyunits_1.1.1             openssl_2.0.6                
##  [71] GetoptLong_1.0.5              haven_2.5.1                  
##  [73] cluster_2.1.4                 fs_1.6.2                     
##  [75] magrittr_2.0.3                data.table_1.14.8            
##  [77] magick_2.7.4                  reprex_2.0.2                 
##  [79] googledrive_2.0.0             hms_1.1.2                    
##  [81] mime_0.12                     evaluate_0.21                
##  [83] xtable_1.8-4                  XML_3.99-0.14                
##  [85] mclust_6.0.0                  readxl_1.4.1                 
##  [87] shape_1.4.6                   compiler_4.2.2               
##  [89] biomaRt_2.50.3                crayon_1.5.2                 
##  [91] minqa_1.2.5                   htmltools_0.5.5              
##  [93] later_1.3.1                   tzdb_0.4.0                   
##  [95] lubridate_1.9.2               DBI_1.1.3                    
##  [97] MASS_7.3-60                   rappdirs_0.3.3               
##  [99] boot_1.3-28.1                 car_3.1-0                    
## [101] cli_3.6.1                     quadprog_1.5-8               
## [103] pkgconfig_2.0.3               GenomicAlignments_1.30.0     
## [105] numDeriv_2016.8-1.1           xml2_1.3.4                   
## [107] annotate_1.72.0               bslib_0.5.0                  
## [109] rngtools_1.5.2                multtest_2.50.0              
## [111] beanplot_1.3.1                rvest_1.0.3                  
## [113] scrime_1.3.5                  doRNG_1.8.6                  
## [115] digest_0.6.31                 base64_2.0.1                 
## [117] cellranger_1.1.0              DelayedMatrixStats_1.16.0    
## [119] restfulr_0.0.15               curl_5.0.1                   
## [121] shiny_1.7.2                   Rsamtools_2.10.0             
## [123] rjson_0.2.21                  nloptr_2.0.3                 
## [125] lifecycle_1.0.3               nlme_3.1-162                 
## [127] jsonlite_1.8.5                Rhdf5lib_1.20.0              
## [129] carData_3.0-5                 askpass_1.1                  
## [131] fansi_1.0.4                   pillar_1.9.0                 
## [133] lattice_0.21-8                survival_3.5-5               
## [135] KEGGREST_1.34.0               fastmap_1.1.1                
## [137] httr_1.4.6                    interactiveDisplayBase_1.32.0
## [139] glue_1.6.2                    png_0.1-8                    
## [141] BiocVersion_3.16.0            bit_4.0.5                    
## [143] stringi_1.7.12                sass_0.4.6                   
## [145] HDF5Array_1.22.1              blob_1.2.4                   
## [147] memoise_2.0.1

Clear

rm(list=ls())