Script 36: Volcano plot of blood results
Compiled: June 27, 2024
This script creates a volcano plot of the DNAm results in blood.
Setup
Load packages
library(tidyverse)
library(ggrepel)
library(GenomicRanges)
library(ggpubr)
library(DNAmArray)Load data
load("../GOTO_Data/GOTO_results-full-blood.Rdata")
dim(limma_base)## [1] 755777 7Save top CpGs
load("../GOTO_Data/GOTO_results-top-blood.Rdata")
cpg_list <- top_cpgs$cpg
length(cpg_list)## [1] 441Volcano plot
Save limits
sig_df <- top_cpgs %>%
filter(
padj_fdr <= 0.05
)
if(nrow(sig_df) >= 1){
sig_limit <- max(sig_df$p)
} else {
sig_limit <- 10E-07
}
min <- as.numeric(min(
abs(limma_base$beta),
na.rm=TRUE))
max <- as.numeric(max(
abs(limma_base$beta),
na.rm=TRUE))
p_max <- as.numeric(-log10(min(limma_base$p, na.rm=TRUE))) + 2
min## [1] 7.781955e-10max## [1] 0.03539525p_max## [1] 10.68585Print volcano plot
plot <- limma_base %>%
ggplot(aes(
x = beta,
y = -log10(p)
)) +
geom_hline(
yintercept = -log10(sig_limit),
linetype = "dashed",
linewidth = 1,
color = "#BBBBBB"
) +
xlim(-0.13, 0.13) +
ylim(0,13) +
geom_point(
color = ifelse(limma_base$p > sig_limit,
"#BBBBBB", ifelse(
abs(limma_base$beta) >= 0.05, "#4477AA", '#4477AA')),
size = ifelse(limma_base$p > sig_limit,
2, ifelse(
abs(limma_base$beta) >= 0.05, 2, 2)),
) +
ggtitle("") +
ylab(bquote(-log[10]~"p")) +
xlab("Effect Size") +
theme(
axis.text = element_text(
size=14,
color="#1B2021"),
axis.title = element_text(
size=16,
hjust=0.5,
color="#1B2021"),
panel.background = element_rect(
fill="white"),
panel.border = element_rect(
color="#1B2021",
fill=NA),
panel.grid.major = element_line(
color="grey95"),
panel.grid.minor = element_line(
color="grey95"),
plot.background = element_rect(
fill="white"))
print(plot)
Save
png('../GOTO_Data/Figures/Figure_5B.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] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggpubr_0.4.0
## [2] ggrepel_0.9.1
## [3] lmerTest_3.1-3
## [4] lme4_1.1-30
## [5] Matrix_1.5-4.1
## [6] IlluminaHumanMethylation450kmanifest_0.4.0
## [7] IlluminaHumanMethylationEPICmanifest_0.3.0
## [8] FlowSorted.BloodExtended.EPIC_1.1.1
## [9] FlowSorted.Blood.EPIC_1.12.1
## [10] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [11] nlme_3.1-162
## [12] quadprog_1.5-8
## [13] genefilter_1.76.0
## [14] ExperimentHub_2.2.1
## [15] AnnotationHub_3.2.2
## [16] BiocFileCache_2.2.1
## [17] dbplyr_2.2.1
## [18] MuSiC_0.2.0
## [19] nnls_1.4
## [20] lubridate_1.9.2
## [21] DNAmArray_2.0.0
## [22] pls_2.8-2
## [23] FDb.InfiniumMethylation.hg19_2.2.0
## [24] org.Hs.eg.db_3.14.0
## [25] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [26] GenomicFeatures_1.46.5
## [27] AnnotationDbi_1.56.2
## [28] minfi_1.40.0
## [29] bumphunter_1.36.0
## [30] locfit_1.5-9.8
## [31] iterators_1.0.14
## [32] foreach_1.5.2
## [33] Biostrings_2.62.0
## [34] XVector_0.34.0
## [35] SummarizedExperiment_1.24.0
## [36] Biobase_2.58.0
## [37] MatrixGenerics_1.10.0
## [38] matrixStats_1.0.0
## [39] GenomicRanges_1.46.1
## [40] GenomeInfoDb_1.34.9
## [41] IRanges_2.32.0
## [42] S4Vectors_0.36.2
## [43] BiocGenerics_0.44.0
## [44] forcats_0.5.2
## [45] stringr_1.5.0
## [46] dplyr_1.1.3
## [47] purrr_0.3.4
## [48] readr_2.1.2
## [49] tidyr_1.2.1
## [50] tibble_3.2.1
## [51] ggplot2_3.4.3
## [52] tidyverse_1.3.2
## [53] cli_3.6.1
## [54] htmltools_0.5.5
## [55] rlang_1.1.1
## [56] rmarkdown_2.16
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.3 tidyselect_1.2.0
## [3] RSQLite_2.2.17 grid_4.2.2
## [5] BiocParallel_1.32.6 cinaR_0.2.3
## [7] munsell_0.5.0 codetools_0.2-19
## [9] preprocessCore_1.60.2 withr_2.5.0
## [11] colorspace_2.1-0 filelock_1.0.2
## [13] highr_0.10 knitr_1.43
## [15] rstudioapi_0.14 ggsignif_0.6.3
## [17] labeling_0.4.2 GenomeInfoDbData_1.2.9
## [19] MCMCpack_1.6-3 farver_2.1.1
## [21] bit64_4.0.5 rhdf5_2.42.1
## [23] coda_0.19-4 vctrs_0.6.3
## [25] generics_0.1.3 xfun_0.39
## [27] timechange_0.2.0 R6_2.5.1
## [29] illuminaio_0.40.0 bitops_1.0-7
## [31] rhdf5filters_1.10.1 cachem_1.0.8
## [33] reshape_0.8.9 DelayedArray_0.24.0
## [35] assertthat_0.2.1 vroom_1.5.7
## [37] promises_1.2.0.1 BiocIO_1.8.0
## [39] scales_1.2.1 googlesheets4_1.0.1
## [41] gtable_0.3.3 mcmc_0.9-7
## [43] MatrixModels_0.5-1 splines_4.2.2
## [45] rstatix_0.7.0 rtracklayer_1.54.0
## [47] gargle_1.5.0 GEOquery_2.62.2
## [49] htm2txt_2.2.2 broom_1.0.1
## [51] abind_1.4-5 BiocManager_1.30.21
## [53] yaml_2.3.7 reshape2_1.4.4
## [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 RColorBrewer_1.1-3
## [63] siggenes_1.68.0 Rcpp_1.0.10
## [65] plyr_1.8.8 sparseMatrixStats_1.10.0
## [67] progress_1.2.2 zlibbioc_1.44.0
## [69] RCurl_1.98-1.12 prettyunits_1.1.1
## [71] openssl_2.0.6 haven_2.5.1
## [73] fs_1.6.2 magrittr_2.0.3
## [75] data.table_1.14.8 SparseM_1.81
## [77] reprex_2.0.2 googledrive_2.0.0
## [79] mime_0.12 hms_1.1.2
## [81] evaluate_0.21 xtable_1.8-4
## [83] XML_3.99-0.14 mclust_6.0.0
## [85] readxl_1.4.1 compiler_4.2.2
## [87] biomaRt_2.50.3 crayon_1.5.2
## [89] minqa_1.2.5 later_1.3.1
## [91] tzdb_0.4.0 DBI_1.1.3
## [93] MASS_7.3-60 rappdirs_0.3.3
## [95] boot_1.3-28.1 car_3.1-0
## [97] pkgconfig_2.0.3 GenomicAlignments_1.30.0
## [99] numDeriv_2016.8-1.1 xml2_1.3.4
## [101] annotate_1.72.0 bslib_0.5.0
## [103] rngtools_1.5.2 multtest_2.50.0
## [105] beanplot_1.3.1 rvest_1.0.3
## [107] doRNG_1.8.6 scrime_1.3.5
## [109] digest_0.6.31 base64_2.0.1
## [111] cellranger_1.1.0 edgeR_3.40.2
## [113] DelayedMatrixStats_1.16.0 restfulr_0.0.15
## [115] curl_5.0.1 shiny_1.7.2
## [117] Rsamtools_2.10.0 quantreg_5.94
## [119] nloptr_2.0.3 rjson_0.2.21
## [121] lifecycle_1.0.3 jsonlite_1.8.5
## [123] Rhdf5lib_1.20.0 carData_3.0-5
## [125] askpass_1.1 limma_3.54.2
## [127] fansi_1.0.4 pillar_1.9.0
## [129] lattice_0.21-8 KEGGREST_1.34.0
## [131] fastmap_1.1.1 httr_1.4.6
## [133] survival_3.5-5 interactiveDisplayBase_1.32.0
## [135] glue_1.6.2 png_0.1-8
## [137] BiocVersion_3.16.0 bit_4.0.5
## [139] stringi_1.7.12 sass_0.4.6
## [141] HDF5Array_1.22.1 blob_1.2.4
## [143] memoise_2.0.1Clear
rm(list=ls())