Correlates of ADOS and methylation - compartment analysis
Mark Ziemann
2023-11-14
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("meth_functions.R")
library("data.table")
library("kableExtra")
library("eulerr")
library("GenomicRanges")
library("qqman")
library("forestplot")
})
source("meth_functions.R")Load Limma 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 51 0.9301106 0.67479681 0.68667348 1.2349086
## 3'UTR 22433 34 1.2129434 0.25108786 0.83420682 1.7089364
## 5'UTR 102937 114 0.8621083 0.15307788 0.70173200 1.0508847
## Body 333689 425 1.0234530 0.73585600 0.89455836 1.1707154
## ExonBnd 6713 2 0.2347955 0.02208455 0.02838447 0.8510027
## TSS1500 115845 168 1.1886017 0.04591205 0.99853936 1.4083218
## TSS200 75104 87 0.9117204 0.44569459 0.72200706 1.1389184
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Intergenic 43381 32 0.5948092 2.631735e-03 0.4040653 0.8466315
## 3'UTR 22433 31 1.1492426 4.336872e-01 0.7754905 1.6447055
## 5'UTR 102937 105 0.8218606 6.471350e-02 0.6633950 1.0095343
## Body 333689 476 1.4132284 5.625484e-07 1.2310194 1.6234086
## ExonBnd 6713 9 1.1107978 7.216921e-01 0.5059310 2.1199885
## TSS1500 115845 134 0.9491165 6.106978e-01 0.7831326 1.1431354
## TSS200 75104 59 0.6235731 2.353519e-04 0.4701899 0.8131389make_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 71 1.4275748 0.005873763 1.10296399 1.8227643
## 3'UTR 22758 18 0.6667923 0.093588036 0.39309639 1.0608189
## 5'UTR 104442 132 1.0921342 0.352381927 0.89972425 1.3172212
## Body 339559 366 0.8551137 0.026309339 0.74365053 0.9826902
## ExonBnd 6862 3 0.3701634 0.075077642 0.07615547 1.0859788
## TSS1500 117509 149 1.0987143 0.304641328 0.91392779 1.3136635
## TSS200 75457 94 1.0699632 0.535910255 0.85376183 1.3277874
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Intergenic 43528 33 0.6594701 1.787720e-02 0.4506115 0.9343429
## 3'UTR 22758 30 1.1784990 3.648760e-01 0.7894626 1.6968371
## 5'UTR 104442 117 0.9948600 1.000000e+00 0.8104788 1.2121831
## Body 339559 439 1.3311840 6.077279e-05 1.1552821 1.5346520
## ExonBnd 6862 12 1.5636266 1.419278e-01 0.8040955 2.7487368
## TSS1500 117509 101 0.7294968 2.652237e-03 0.5860767 0.8999640
## TSS200 75457 67 0.7696615 3.862801e-02 0.5900269 0.9894609make_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 159 0.8060765 0.012332080 0.6760051 0.9561060
## OpenSea 442427 570 1.0434181 0.523897750 0.9189347 1.1854343
## Shelf 55491 53 0.7412019 0.034887781 0.5510427 0.9779716
## Shore 142558 218 1.2677690 0.002311228 1.0859120 1.4751269
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 150182 162 0.8242457 0.0238509853 0.692205 0.9764387
## OpenSea 442427 518 0.8457016 0.0089327487 0.745533 0.9594641
## Shelf 55491 95 1.3913745 0.0035060520 1.113978 1.7207079
## Shore 142558 225 1.3203745 0.0003369766 1.133018 1.5337616make_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 177 0.9342229 0.44077756 0.7895856 1.100408
## OpenSea 452360 580 1.0694516 0.30733742 0.9415295 1.215564
## Shelf 56174 55 0.7731764 0.06278458 0.5779977 1.015328
## Shore 143927 188 1.0597333 0.48313995 0.8993294 1.243516
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 150186 92 0.4398291 7.197931e-17 0.3509793 0.5455411
## OpenSea 452360 698 1.7909408 3.305622e-18 1.5623637 2.0568003
## Shelf 56174 63 0.8933467 4.201049e-01 0.6807654 1.1538172
## Shore 143927 147 0.7885126 7.296286e-03 0.6572722 0.9404946make_forest_plots_up(comp)
make_forest_plots_dn(comp)
Session information
For reproducibility
sessionInfo()## R version 4.3.1 (2023-06-16)
## 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/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] forestplot_3.1.3
## [2] abind_1.4-5
## [3] checkmate_2.3.0
## [4] qqman_0.1.9
## [5] eulerr_7.0.0
## [6] data.table_1.14.8
## [7] DMRcatedata_2.18.0
## [8] IlluminaHumanMethylation450kmanifest_0.4.0
## [9] IlluminaHumanMethylationEPICmanifest_0.3.0
## [10] vioplot_0.4.0
## [11] zoo_1.8-12
## [12] sm_2.2-5.7.1
## [13] mitch_1.12.0
## [14] FlowSorted.Blood.EPIC_2.4.2
## [15] ExperimentHub_2.8.1
## [16] AnnotationHub_3.8.0
## [17] BiocFileCache_2.11.1
## [18] dbplyr_2.4.0
## [19] DMRcate_2.14.1
## [20] reshape2_1.4.4
## [21] FlowSorted.Blood.450k_1.38.0
## [22] WGCNA_1.72-1
## [23] fastcluster_1.2.3
## [24] dynamicTreeCut_1.63-1
## [25] limma_3.56.2
## [26] missMethyl_1.34.0
## [27] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [28] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [29] minfi_1.46.0
## [30] bumphunter_1.42.0
## [31] locfit_1.5-9.8
## [32] iterators_1.0.14
## [33] foreach_1.5.2
## [34] Biostrings_2.68.1
## [35] XVector_0.40.0
## [36] SummarizedExperiment_1.30.2
## [37] Biobase_2.60.0
## [38] MatrixGenerics_1.12.3
## [39] matrixStats_1.1.0
## [40] GenomicRanges_1.52.1
## [41] GenomeInfoDb_1.36.4
## [42] IRanges_2.34.1
## [43] S4Vectors_0.38.2
## [44] BiocGenerics_0.46.0
## [45] beeswarm_0.4.0
## [46] ggplot2_3.4.4
## [47] gplots_3.1.3
## [48] RColorBrewer_1.1-3
## [49] dplyr_1.1.3
## [50] kableExtra_1.3.4
##
## 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 doParallel_1.0.17
## [7] tools_4.3.1 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.2
## [15] rhdf5filters_1.12.1 permute_0.9-7
## [17] withr_2.5.0 GGally_2.1.2
## [19] prettyunits_1.1.1 gridExtra_2.3
## [21] base64_2.0.1 preprocessCore_1.61.0
## [23] cli_3.6.1 labeling_0.4.3
## [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.2
## [35] R.utils_2.12.2 dichromat_2.0-0.1
## [37] scrime_1.3.5 BSgenome_1.68.0
## [39] readxl_1.4.3 rstudioapi_0.15.0
## [41] impute_1.74.1 RSQLite_2.3.3
## [43] generics_0.1.3 BiocIO_1.10.0
## [45] gtools_3.9.4 GO.db_3.17.0
## [47] Matrix_1.6-1 interp_1.1-4
## [49] fansi_1.0.4 R.methodsS3_1.8.2
## [51] lifecycle_1.0.3 edgeR_3.42.4
## [53] yaml_2.3.7 rhdf5_2.44.0
## [55] blob_1.2.4 promises_1.2.1
## [57] crayon_1.5.2 lattice_0.21-8
## [59] echarts4r_0.4.5 GenomicFeatures_1.52.2
## [61] annotate_1.78.0 KEGGREST_1.40.1
## [63] pillar_1.9.0 knitr_1.43
## [65] beanplot_1.3.1 rjson_0.2.21
## [67] codetools_0.2-19 glue_1.6.2
## [69] vctrs_0.6.3 png_0.1-8
## [71] cellranger_1.1.0 gtable_0.3.4
## [73] cachem_1.0.8 xfun_0.40
## [75] mime_0.12 S4Arrays_1.0.6
## [77] survival_3.5-7 statmod_1.5.0
## [79] ellipsis_0.3.2 interactiveDisplayBase_1.38.0
## [81] nlme_3.1-163 bit64_4.0.5
## [83] bsseq_1.36.0 progress_1.2.2
## [85] filelock_1.0.2 bslib_0.5.1
## [87] nor1mix_1.3-0 KernSmooth_2.23-22
## [89] rpart_4.1.19 colorspace_2.1-0
## [91] DBI_1.1.3 Hmisc_5.1-1
## [93] nnet_7.3-19 tidyselect_1.2.0
## [95] bit_4.0.5 compiler_4.3.1
## [97] curl_5.0.2 rvest_1.0.3
## [99] htmlTable_2.4.2 BiasedUrn_2.0.11
## [101] xml2_1.3.5 DelayedArray_0.26.7
## [103] rtracklayer_1.60.1 scales_1.2.1
## [105] caTools_1.18.2 quadprog_1.5-8
## [107] rappdirs_0.3.3 stringr_1.5.0
## [109] digest_0.6.33 rmarkdown_2.24
## [111] GEOquery_2.68.0 htmltools_0.5.6
## [113] pkgconfig_2.0.3 jpeg_0.1-10
## [115] base64enc_0.1-3 sparseMatrixStats_1.12.2
## [117] highr_0.10 fastmap_1.1.1
## [119] ensembldb_2.24.1 rlang_1.1.1
## [121] htmlwidgets_1.6.2 shiny_1.7.5
## [123] DelayedMatrixStats_1.22.6 farver_2.1.1
## [125] jquerylib_0.1.4 jsonlite_1.8.7
## [127] BiocParallel_1.34.2 mclust_6.0.0
## [129] R.oo_1.25.0 VariantAnnotation_1.46.0
## [131] RCurl_1.98-1.13 magrittr_2.0.3
## [133] Formula_1.2-5 GenomeInfoDbData_1.2.10
## [135] Rhdf5lib_1.22.1 munsell_0.5.0
## [137] Rcpp_1.0.11 stringi_1.7.12
## [139] zlibbioc_1.46.0 MASS_7.3-60
## [141] plyr_1.8.9 org.Hs.eg.db_3.17.0
## [143] deldir_1.0-9 splines_4.3.1
## [145] multtest_2.56.0 hms_1.1.3
## [147] rngtools_1.5.2 biomaRt_2.56.1
## [149] BiocVersion_3.17.1 XML_3.99-0.15
## [151] evaluate_0.21 calibrate_1.7.7
## [153] latticeExtra_0.6-30 biovizBase_1.48.0
## [155] BiocManager_1.30.22 httpuv_1.6.11
## [157] tzdb_0.4.0 tidyr_1.3.0
## [159] openssl_2.1.0 purrr_1.0.2
## [161] reshape_0.8.9 xtable_1.8-4
## [163] restfulr_0.0.15 AnnotationFilter_1.24.0
## [165] later_1.3.1 viridisLite_0.4.2
## [167] tibble_3.2.1 memoise_2.0.1
## [169] AnnotationDbi_1.62.2 GenomicAlignments_1.36.0
## [171] cluster_2.1.4