Script 23: Dotplots of selected CpGs in muscle
Compiled: June 26, 2024
This script saves dot plots of the selected CpG, gene, and trait in muscle for Figure 2B.
Setup
Load packages
library(tidyverse)
library(ggrepel)
library(GenomicRanges)
library(ggpubr)
library(DNAmArray)
library(SummarizedExperiment)
library(limma)
library(edgeR)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 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("cg21005024")
dim(betas)## [1] 160 1Phenotype 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))
ID_list <- targets %>% dplyr::select(ID)Save selected CpG
targets$cg21005024 <- betas[,1]CpG Figure
Make box plot
plot <- targets %>%
ggplot(aes(x = timepoint, y = cg21005024,
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_2B-cg21005024.png')
print(plot)
dev.off()## png
## 2Gene Figure
Erik scripts
source('../GOTO_Data/RNAseq/goto.rnaseq.functions.R')Load RNAseq
pathIN_dat <- "../GOTO_Data/RNAseq/merge.gene.counts_biopet_13052016.RData"
pathIN_cov <- "../GOTO_Data/RNAseq/muscle_QC_covariates_filesv2.csv"
filt.samp <- "tissue_muscle|qc_sexswitch|qc_multdim2|qc_rep1|complete_pairs"
dat1 <- read.gotornaseq(pathIN_dat = pathIN_dat, pathIN_cov, filt.samp = filt.samp, type = 'voom-export', quiet = FALSE)## ||| PREPARING GOTO RNASEQ DATA
## || READING DATA
## | Loading RNASEQ .. OK!
## [555 samples x 56520 features]
## | Reading COVARIATES .. OK!
## [maintaining 174 samples x 117 features]
## | Merging data .. OK!
## [555 samples x 56637 features]
## || SUBSETTING SAMPLES
## | Subsetting SAMPLES on ['tissue_muscle']; PASS: 174 out of 555
## | Subsetting SAMPLES on ['qc_sexswitch']; PASS: 172 out of 174
## | Subsetting SAMPLES on ['qc_multdim2']; PASS: 168 out of 172
## | Subsetting SAMPLES on ['qc_rep1']; PASS: 168 out of 168
## | Subsetting SAMPLES on ['complete_pairs']; PASS: 162 out of 168
## || TRANSFORMING DATA
## | VOOM .. OK!
## | DONE!goto_exp <- dat1[["dat"]]
goto_exp <- goto_exp %>%
dplyr::select(sampID2, intervention, flowcell,
starts_with('ENS')) %>%
mutate(
ID = str_c(sampID2, '_', intervention)
)
exp_df <- left_join(ID_list, goto_exp, by = 'ID')
dim(exp_df)## [1] 160 56519ens2gene <- cinaR::grch37Print dotplot
name <- "ENSG00000106070"
df <- exp_df %>% filter(!is.na(intervention))
df$sampID2 <- as.factor(df$sampID2)
df$intervention <- as.factor(df$intervention)
plot <- df %>%
ggplot(aes(x = intervention, y = ENSG00000106070,
fill = intervention)) +
geom_line(aes(group = sampID2), color = '#DDDDDD', linewidth=1) +
geom_dotplot(binaxis = "y", binwidth = 0.12,
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_2B-GRB10.png')
print(plot)
dev.off()## png
## 2Trait Figure
Body fat
Print plot
plot <- targets %>%
ggplot(aes(x = timepoint, y = wc,
fill = timepoint)) +
geom_line(aes(group = IOP2_ID), color = '#DDDDDD', linewidth=1) +
geom_dotplot(binaxis = "y", binwidth = 3,
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)
WC
Save
png('../GOTO_Data/Figures/Figure_2B-wc.png')
print(plot)
dev.off()## png
## 2Print plot
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.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_2B-bodyfat.png')
print(plot)
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())