Correlates of ADOS and methylation - compartment analysis
Mark Ziemann
2023-11-08
Source: https://github.com/markziemann/asd_meth
Introduction
Here I will be running a comparison of differential methylation data from guthrie and fresh blood samples.
Here are the files that I’m using:
limma_guthrie_ADOS.csv
limma_blood_ADOS.csv
suppressPackageStartupMessages({
library("parallel")
library("mitch")
library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
source("https://raw.githubusercontent.com/markziemann/ART_methylation/master/meth_functions.R")
library("data.table")
library("kableExtra")
library("eulerr")
library("RIdeogram")
library("GenomicRanges")
library("tictoc")
})
source("meth_functions.R")Load Limma and RUV data for manhattan plot
First I will generate plots for limmma analysis.
limma_guthrie_ADOS.csv
limma_blood_ADOS.csv
# blood at assessment
top <- read.csv("limma_blood_ADOS.csv")
nrow(top)## [1] 802647top <- subset(top,P.Value<1e-2)
nrow(top)## [1] 3573-log10(min(top$P.Value))## [1] 4.641884top$chr <- as.integer(gsub("chr","",top$chr))
top$snp <- paste(top$Name,top$UCSC_RefGene_Name)
top$snp <- sapply(strsplit(top$snp,";"),"[[",1)
up <- subset(top,logFC>0)
dn <- subset(top,logFC<0)
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 ,main="Blood at assessment limma hypermethylated")
manhattan(x=dn,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 , main="Blood at assessment limma hypomethylated")
pdf("manhat_limma_bl.pdf")
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 ,main="Blood at assessment limma hypermethylated")
manhattan(x=dn,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 , main="Blood at assessment limma hypomethylated")
dev.off()## png
## 2# neonatal guthrie card
top <- read.csv("limma_guthrie_ADOS.csv")
nrow(top)## [1] 790658top <- subset(top,P.Value<1e-2)
nrow(top)## [1] 3052-log10(min(top$P.Value))## [1] 5.396928top$chr <- as.integer(gsub("chr","",top$chr))
top$snp <- paste(top$Name,top$UCSC_RefGene_Name)
top$snp <- sapply(strsplit(top$snp,";"),"[[",1)
up <- subset(top,logFC>0)
dn <- subset(top,logFC<0)
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 ,main="Neonatal Guthrie card limma hypermethylated",
annotateTop=FALSE)
manhattan(x=dn,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 , main="Neonatal Guthrie card limma hypomethylated",
annotateTop=FALSE)
pdf("manhat_limma_gu.pdf")
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 ,main="Neonatal Guthrie card limma hypermethylated",
annotateTop=FALSE)
manhattan(x=dn,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 , main="Neonatal Guthrie card limma hypomethylated",
annotateTop=FALSE)
dev.off()## png
## 2Compartment enrichment
Looking for enrichments in different genomic compartments.
Guthrie cards.
par(mfrow=c(2,1))
# guthrie
dma1 <- read.csv("limma_guthrie_ADOS.csv")
if (nrow(subset(dma1,P.Value <1e-4)) > 100 ) {
comp <- compartment_enrichment(dma1)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp
make_forest_plots_up(comp)
make_forest_plots_dn(comp)
}Blood at assessment.
# blood
dma2 <- read.csv("limma_blood_ADOS.csv")
if (nrow(subset(dma2,P.Value <1e-4)) > 100 ) {
comp <- compartment_enrichment(dma2)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp
make_forest_plots_up(comp)
make_forest_plots_dn(comp)
}Compartments top 1000
Guthrie cards.
par(mfrow=c(2,1))
# guthrie
comp <- compartment_enrichment2(dma1)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI upperCI
## Intergenic 43381 23 0.8389085 0.48706693 0.5258383 1.276341
## 3'UTR 22433 18 1.2948872 0.27561695 0.7597998 2.072428
## 5'UTR 102937 51 0.7643920 0.07889085 0.5592212 1.025627
## Body 333689 227 1.1814543 0.08494842 0.9751244 1.431904
## TSS1500 115845 81 1.1443742 0.27418805 0.8876763 1.460454
## TSS200 75104 38 0.7904865 0.18865183 0.5511121 1.104070
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Intergenic 43381 15 0.5761750 0.0308946406 0.3193848 0.9621765
## 3'UTR 22433 16 1.2300273 0.3975719033 0.6960289 2.0232371
## 5'UTR 102937 50 0.8078831 0.1629589719 0.5884785 1.0884038
## Body 333689 234 1.4685899 0.0001316616 1.2021212 1.7968522
## ExonBnd 6713 4 1.0201702 0.8003323414 0.2765425 2.6334465
## TSS1500 115845 63 0.9182427 0.5942944466 0.6905059 1.2041218
## TSS200 75104 27 0.5880403 0.0050933600 0.3824177 0.8693314make_forest_plots_up(comp)
make_forest_plots_dn(comp)
Blood at assessment.
# blood
comp <- compartment_enrichment2(dma2)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI upperCI
## Intergenic 43528 26 1.4859557 0.06636519 0.953038174 2.226311
## 3'UTR 22758 5 0.5224343 0.18233626 0.168373824 1.233495
## 5'UTR 104442 51 1.2171947 0.21605071 0.881706511 1.652061
## Body 339559 125 0.8070945 0.07043460 0.635256389 1.023158
## ExonBnd 6862 1 0.3496924 0.53965923 0.008825083 1.964196
## TSS1500 117509 53 1.1091360 0.48032288 0.807675308 1.498631
## TSS200 75457 33 1.0634901 0.70520650 0.716808409 1.531364
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Intergenic 43528 25 0.7417994 0.151754452 0.4755129 1.1077977
## 3'UTR 22758 23 1.3413967 0.176996771 0.8425370 2.0355856
## 5'UTR 104442 80 1.0063621 0.951579367 0.7834649 1.2785220
## Body 339559 289 1.2517246 0.009742042 1.0534633 1.4881382
## ExonBnd 6862 10 1.9314000 0.045232660 0.9203789 3.5825126
## TSS1500 117509 67 0.7126767 0.009075108 0.5433220 0.9220956
## TSS200 75457 47 0.8001268 0.162291260 0.5800384 1.0806029make_forest_plots_up(comp)
make_forest_plots_dn(comp)
Looking for enrichments in proximity to CGI contexts
Guthrie cards.
par(mfrow=c(2,1))
if (nrow(subset(dma1,P.Value <1e-4)) > 100 ) {
# guthrie
comp <- cgi_enrichment(dma1)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp
make_forest_plots_up(comp)
make_forest_plots_dn(comp)
}Blood at assessment.
if (nrow(subset(dma2,P.Value <1e-4)) > 100 ) {
# blood
comp <- cgi_enrichment(dma2)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp
make_forest_plots_up(comp)
make_forest_plots_dn(comp)
}CGI enrichments in top 1000 probes
Guthrie cards.
par(mfrow=c(2,1))
# guthrie
comp <- cgi_enrichment2(dma1)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI upperCI
## Island 150182 80 0.7914758 0.05507706 0.6154100 1.006919
## OpenSea 442427 291 1.0411421 0.65617216 0.8708430 1.246143
## Shelf 55491 30 0.8262023 0.34108662 0.5512779 1.195275
## Shore 142558 110 1.2472827 0.04379063 1.0006516 1.543986
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 150182 83 0.8717702 0.27330781 0.6799065 1.106451
## OpenSea 442427 255 0.8576492 0.09192572 0.7153376 1.028657
## Shelf 55491 48 1.4423707 0.02080876 1.0475126 1.946126
## Shore 142558 103 1.2132335 0.08770710 0.9665051 1.511589make_forest_plots_up(comp)
make_forest_plots_dn(comp)
Blood at assessment.
# blood
comp <- cgi_enrichment2(dma2)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI upperCI
## Island 150186 57 0.9523840 0.773732545 0.7020288 1.272907
## OpenSea 452360 167 0.8620659 0.192891313 0.6872938 1.082060
## Shelf 56174 16 0.7062707 0.187293830 0.3985352 1.166031
## Shore 143927 77 1.4686262 0.004221713 1.1209174 1.905982
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 150186 60 0.4181603 6.780836e-13 0.3153031 0.5456317
## OpenSea 452360 478 1.8064067 2.870589e-13 1.5304745 2.1381921
## Shelf 56174 41 0.8485316 3.299684e-01 0.6030135 1.1643462
## Shore 143927 104 0.8219508 6.497237e-02 0.6605231 1.0143919make_forest_plots_up(comp)
make_forest_plots_dn(comp)
Session information
For reproducibility
sessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## 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
##
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats4 parallel stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] RIdeogram_0.2.2
## [2] ggplot2_3.4.3
## [3] dplyr_1.1.3
## [4] tictoc_1.2
## [5] kableExtra_1.3.4
## [6] data.table_1.14.8
## [7] ENmix_1.36.03
## [8] doParallel_1.0.17
## [9] qqman_0.1.9
## [10] RCircos_1.2.2
## [11] beeswarm_0.4.0
## [12] forestplot_3.1.3
## [13] abind_1.4-5
## [14] checkmate_2.2.0
## [15] reshape2_1.4.4
## [16] gplots_3.1.3
## [17] eulerr_7.0.0
## [18] GEOquery_2.68.0
## [19] RColorBrewer_1.1-3
## [20] IlluminaHumanMethylation450kmanifest_0.4.0
## [21] topconfects_1.16.0
## [22] DMRcatedata_2.18.0
## [23] ExperimentHub_2.8.1
## [24] AnnotationHub_3.8.0
## [25] BiocFileCache_2.8.0
## [26] dbplyr_2.3.3
## [27] DMRcate_2.14.1
## [28] limma_3.56.2
## [29] missMethyl_1.34.0
## [30] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [31] R.utils_2.12.2
## [32] R.oo_1.25.0
## [33] R.methodsS3_1.8.2
## [34] plyr_1.8.8
## [35] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [36] minfi_1.46.0
## [37] bumphunter_1.42.0
## [38] locfit_1.5-9.8
## [39] iterators_1.0.14
## [40] foreach_1.5.2
## [41] Biostrings_2.68.1
## [42] XVector_0.40.0
## [43] SummarizedExperiment_1.30.2
## [44] Biobase_2.60.0
## [45] MatrixGenerics_1.12.3
## [46] matrixStats_1.0.0
## [47] GenomicRanges_1.52.0
## [48] GenomeInfoDb_1.36.3
## [49] IRanges_2.34.1
## [50] S4Vectors_0.38.1
## [51] BiocGenerics_0.46.0
## [52] mitch_1.12.0
##
## loaded via a namespace (and not attached):
## [1] DSS_2.48.0 ProtGenerics_1.32.0
## [3] bitops_1.0-7 httr_1.4.7
## [5] webshot_0.5.5 dynamicTreeCut_1.63-1
## [7] tools_4.3.2 doRNG_1.8.6
## [9] backports_1.4.1 utf8_1.2.3
## [11] R6_2.5.1 HDF5Array_1.28.1
## [13] lazyeval_0.2.2 Gviz_1.44.1
## [15] rhdf5filters_1.12.1 permute_0.9-7
## [17] withr_2.5.0 prettyunits_1.1.1
## [19] GGally_2.1.2 gridExtra_2.3
## [21] grImport2_0.2-0 base64_2.0.1
## [23] preprocessCore_1.62.1 cli_3.6.1
## [25] sass_0.4.7 readr_2.1.4
## [27] genefilter_1.82.1 askpass_1.2.0
## [29] Rsamtools_2.16.0 systemfonts_1.0.4
## [31] foreign_0.8-85 siggenes_1.74.0
## [33] illuminaio_0.42.0 svglite_2.1.1
## [35] dichromat_2.0-0.1 scrime_1.3.5
## [37] BSgenome_1.68.0 readxl_1.4.3
## [39] impute_1.74.1 rstudioapi_0.15.0
## [41] RSQLite_2.3.1 generics_0.1.3
## [43] BiocIO_1.10.0 gtools_3.9.4
## [45] Matrix_1.6-1.1 interp_1.1-4
## [47] fansi_1.0.4 lifecycle_1.0.3
## [49] edgeR_3.42.4 yaml_2.3.7
## [51] rhdf5_2.44.0 blob_1.2.4
## [53] promises_1.2.1 crayon_1.5.2
## [55] lattice_0.22-5 echarts4r_0.4.5
## [57] GenomicFeatures_1.52.2 annotate_1.78.0
## [59] KEGGREST_1.40.0 pillar_1.9.0
## [61] knitr_1.44 beanplot_1.3.1
## [63] rjson_0.2.21 codetools_0.2-19
## [65] glue_1.6.2 vctrs_0.6.3
## [67] png_0.1-8 cellranger_1.1.0
## [69] gtable_0.3.4 assertthat_0.2.1
## [71] cachem_1.0.8 xfun_0.40
## [73] S4Arrays_1.0.6 mime_0.12
## [75] survival_3.5-7 statmod_1.5.0
## [77] interactiveDisplayBase_1.38.0 ellipsis_0.3.2
## [79] nlme_3.1-163 bit64_4.0.5
## [81] bsseq_1.36.0 progress_1.2.2
## [83] filelock_1.0.2 bslib_0.5.1
## [85] nor1mix_1.3-0 KernSmooth_2.23-22
## [87] rpart_4.1.21 colorspace_2.1-0
## [89] DBI_1.1.3 Hmisc_5.1-1
## [91] nnet_7.3-19 tidyselect_1.2.0
## [93] bit_4.0.5 compiler_4.3.2
## [95] curl_5.0.2 rvest_1.0.3
## [97] htmlTable_2.4.1 xml2_1.3.5
## [99] RPMM_1.25 DelayedArray_0.26.7
## [101] rtracklayer_1.60.1 scales_1.2.1
## [103] caTools_1.18.2 quadprog_1.5-8
## [105] rappdirs_0.3.3 stringr_1.5.0
## [107] digest_0.6.33 rmarkdown_2.24
## [109] htmltools_0.5.6 pkgconfig_2.0.3
## [111] jpeg_0.1-10 base64enc_0.1-3
## [113] sparseMatrixStats_1.12.2 highr_0.10
## [115] fastmap_1.1.1 ensembldb_2.24.0
## [117] rlang_1.1.1 htmlwidgets_1.6.2
## [119] shiny_1.7.5 DelayedMatrixStats_1.22.6
## [121] jquerylib_0.1.4 jsonlite_1.8.7
## [123] BiocParallel_1.34.2 mclust_6.0.0
## [125] VariantAnnotation_1.46.0 RCurl_1.98-1.12
## [127] magrittr_2.0.3 Formula_1.2-5
## [129] GenomeInfoDbData_1.2.10 Rhdf5lib_1.22.1
## [131] munsell_0.5.0 Rcpp_1.0.11
## [133] stringi_1.7.12 zlibbioc_1.46.0
## [135] MASS_7.3-60 org.Hs.eg.db_3.17.0
## [137] deldir_1.0-9 splines_4.3.2
## [139] multtest_2.56.0 hms_1.1.3
## [141] rngtools_1.5.2 geneplotter_1.78.0
## [143] biomaRt_2.56.1 BiocVersion_3.17.1
## [145] XML_3.99-0.14 evaluate_0.21
## [147] calibrate_1.7.7 latticeExtra_0.6-30
## [149] biovizBase_1.48.0 BiocManager_1.30.22
## [151] tzdb_0.4.0 httpuv_1.6.11
## [153] tidyr_1.3.0 openssl_2.1.0
## [155] purrr_1.0.2 reshape_0.8.9
## [157] xtable_1.8-4 restfulr_0.0.15
## [159] AnnotationFilter_1.24.0 rsvg_2.4.0
## [161] later_1.3.1 viridisLite_0.4.2
## [163] tibble_3.2.1 memoise_2.0.1
## [165] AnnotationDbi_1.62.2 GenomicAlignments_1.36.0
## [167] cluster_2.1.4