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_blood_ADOS.csv

  • limma_guthrie_ADOS.csv

  • limma_buccal_ADOS.csv

suppressPackageStartupMessages({
  library("limma")
  library("parallel")
  library("mitch")
  library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
  source("meth_functions.R")
  library("data.table")
  library("kableExtra")
  library("eulerr")
  library("GenomicRanges")
  data("dualmap850kEID")
})
## Warning in data("dualmap850kEID"): data set 'dualmap850kEID' not found
source("meth_functions.R")

Load data

bl_mvals <- readRDS(file="bl_mvals.rds")
gu_mvals <- readRDS(file="gu_mvals.rds")
bu_mvals <- readRDS(file="buccal_mvals.rds")

bl_design <- readRDS(file="bl_design_ados.rds")
gu_design <- readRDS(file="gu_design_ados.rds")
bu_design <- readRDS(file="buccal_design_ados.rds")

bl_lim <- read.csv("limma_blood_ADOS.csv")
gu_lim <- read.csv("limma_guthrie_ADOS.csv")
bu_lim <- read.csv("limma_buccal_ADOS.csv")

Probe sets

For each gene, extract out the probes.

anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","Regulatory_Feature_Group","Islands_Name","Relation_to_Island")])
promoters <- grep("Prom",myann$Regulatory_Feature_Group)
promoters <- myann[promoters,]
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
sets <- split(rep(names(gp2), lengths(gp2)), unlist(gp2))
summary(unlist(lapply(sets,length)))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    9.00   24.00   49.68   55.00 6778.00

Gene methylation enrichment for blood

For each gene, calculate the mean probe t-statistic.

bl <- pmea(mval=bl_mvals,design=bl_design,sets,cores=10)
str(bl)
## List of 2
##  $ df      :'data.frame':    27364 obs. of  6 variables:
##   ..$ nprobes: num [1:27364] 776 504 394 212 416 282 403 226 148 296 ...
##   ..$ mean   : num [1:27364] -0.426 -0.462 -0.51 -0.596 -0.474 ...
##   ..$ median : num [1:27364] -0.436 -0.575 -0.58 -0.673 -0.544 ...
##   ..$ P.Value: num [1:27364] 5.89e-35 2.87e-28 3.36e-27 4.15e-25 9.03e-24 ...
##   ..$ sig    : num [1:27364] 34.2 27.5 26.5 24.4 23 ...
##   ..$ FDR    : num [1:27364] 1.61e-30 7.84e-24 9.18e-23 1.13e-20 2.47e-19 ...
##  $ toptable:'data.frame':    802647 obs. of  6 variables:
##   ..$ logFC    : num [1:802647] -0.0855 0.1005 -0.0706 0.1229 -0.0612 ...
##   ..$ AveExpr  : num [1:802647] -4.336 4.685 -4.504 0.666 0.655 ...
##   ..$ t        : num [1:802647] -6.71 5.98 -5.8 5.77 -5.61 ...
##   ..$ P.Value  : num [1:802647] 2.28e-05 6.70e-05 8.77e-05 9.21e-05 1.19e-04 ...
##   ..$ adj.P.Val: num [1:802647] 0.876 0.876 0.876 0.876 0.876 ...
##   ..$ B        : num [1:802647] 3.05 2.04 1.78 1.74 1.49 ...
blm <- bl$df[,"mean",drop=FALSE]
blm <- blm[!is.na(blm$mean),,drop=FALSE] # remove NA vals
head(blm)
##               mean
## MAD1L1  -0.4256926
## HDAC4   -0.4618768
## ATP11A  -0.5100916
## GNAS    -0.5959431
## TRAPPC9 -0.4741987
## PDE4D   -0.6894219
nrow(blm)
## [1] 25739
par(mfrow=c(2,1))
hist(bl$toptable$t,xlab="probe t-stat",main="blood at assessment", breaks=30) ; grid() ; abline(v=0,lty=1,lwd=3)
hist(bl$df$mean,xlab="gene mean t-stat",main="blood at assessment", breaks=30) ; grid() ; abline(v=0,lty=1,lwd=3)

pdf("hist_bl_ados.pdf")
par(mfrow=c(2,1))
hist(bl$toptable$t,xlab="probe t-stat",main="blood at assessment",breaks=30) ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
hist(bl$df$mean,xlab="gene mean t-stat",main="blood at assessment", breaks=30) ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
dev.off()
## png 
##   2
par(mfrow=c(1,1))

gmea_boxplot(bl,sets)

Gene methylation enrichment for guthrie cards

gu <- pmea(mval=gu_mvals,design=gu_design,sets,cores=10)
str(gu)
## List of 2
##  $ df      :'data.frame':    27364 obs. of  6 variables:
##   ..$ nprobes: num [1:27364] 1395 153 68 81 88 ...
##   ..$ mean   : num [1:27364] -0.252 -0.458 -0.627 -0.608 -0.574 ...
##   ..$ median : num [1:27364] -0.235 -0.457 -0.51 -0.625 -0.624 ...
##   ..$ P.Value: num [1:27364] 1.72e-22 1.24e-12 3.40e-12 4.72e-12 4.91e-12 ...
##   ..$ sig    : num [1:27364] 21.8 11.9 11.5 11.3 11.3 ...
##   ..$ FDR    : num [1:27364] 4.72e-18 3.39e-08 9.30e-08 1.29e-07 1.34e-07 ...
##  $ toptable:'data.frame':    790658 obs. of  6 variables:
##   ..$ logFC    : num [1:790658] 0.0827 0.1035 0.1318 0.0916 -0.1333 ...
##   ..$ AveExpr  : num [1:790658] 3.05 1.8 1.6 -3.38 2.27 ...
##   ..$ t        : num [1:790658] 7.39 6.87 6.82 6.72 -6.34 ...
##   ..$ P.Value  : num [1:790658] 4.01e-06 8.92e-06 9.54e-06 1.12e-05 2.08e-05 ...
##   ..$ adj.P.Val: num [1:790658] 1 1 1 1 1 ...
##   ..$ B        : num [1:790658] 4.58 3.85 3.79 3.64 3.07 ...
gum <- gu$df[,"mean",drop=FALSE]
gum <- gum[!is.na(gum$mean),,drop=FALSE] # remove NA vals
head(gum)
##                mean
## PTPRN2   -0.2524635
## ANK3     -0.4580590
## MIR548F3 -0.6267218
## AKAP6    -0.6084892
## LRP1B    -0.5737076
## RNU6-71P -0.5286560
nrow(gum)
## [1] 25716
par(mfrow=c(2,1))
hist(gu$toptable$t,xlab="probe t-stat",main="neonatal Guthrie card",breaks=30) ; grid() ; abline(v=0,lty=1,lwd=3)
hist(gu$df$mean,xlab="gene mean t-stat",main="neonatal Guthrie card",breaks=30) ; grid() ; abline(v=0,lty=1,lwd=3)

pdf("hist_gu_ados.pdf")
par(mfrow=c(2,1))
hist(gu$toptable$t,xlab="probe t-stat",main="neonatal Guthrie card", breaks=30) ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
hist(gu$df$mean,xlab="gene mean t-stat",main="neonatal Guthrie card", breaks=30) ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
dev.off()
## png 
##   2
par(mfrow=c(1,1))

gmea_boxplot(gu,sets)

Gene methylation enrichment for buccal samples

bu <- pmea(mval=bu_mvals,design=bu_design,sets,cores=10)
str(bu)
## List of 2
##  $ df      :'data.frame':    27364 obs. of  6 variables:
##   ..$ nprobes: num [1:27364] 656 382 257 376 212 475 114 121 84 48 ...
##   ..$ mean   : num [1:27364] 0.65 0.7 0.643 0.559 0.654 ...
##   ..$ median : num [1:27364] 0.715 0.724 0.646 0.598 0.83 ...
##   ..$ P.Value: num [1:27364] 2.02e-53 2.70e-39 2.99e-25 2.66e-21 2.02e-19 ...
##   ..$ sig    : num [1:27364] 52.7 38.6 24.5 20.6 18.7 ...
##   ..$ FDR    : num [1:27364] 5.53e-49 7.39e-35 8.18e-21 7.28e-17 5.53e-15 ...
##  $ toptable:'data.frame':    801260 obs. of  6 variables:
##   ..$ logFC    : num [1:801260] 0.276 -0.187 0.171 -0.178 -0.177 ...
##   ..$ AveExpr  : num [1:801260] -3.79 5.52 -3.96 1.34 1.6 ...
##   ..$ t        : num [1:801260] 12.16 -10.05 9.5 -9.41 -9.01 ...
##   ..$ P.Value  : num [1:801260] 8.38e-06 2.82e-05 4.01e-05 4.27e-05 5.60e-05 ...
##   ..$ adj.P.Val: num [1:801260] 0.905 0.905 0.905 0.905 0.905 ...
##   ..$ B        : num [1:801260] 3.1 2.41 2.19 2.15 1.97 ...
bum <- bu$df[,"mean",drop=FALSE]
bum <- bum[!is.na(bum$mean),,drop=FALSE] # remove NA vals
head(bum)
##             mean
## PRDM16 0.6501046
## CDH4   0.6997851
## SEPT9  0.6428529
## KCNQ1  0.5587846
## GNAS   0.6544056
## SHANK2 0.4526450
nrow(bum)
## [1] 25739
par(mfrow=c(2,1))
hist(bu$toptable$t,xlab="probe t-stat",main="buccal swab") ; grid() ; abline(v=0,lty=1,lwd=3)
hist(bu$df$mean,xlab="gene mean t-stat",main="buccal swab") ; grid() ; abline(v=0,lty=1,lwd=3)

pdf("hist_bu_ados.pdf")
par(mfrow=c(2,1))
hist(bu$toptable$t,xlab="probe t-stat",main="buccal swab") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
hist(bu$df$mean,xlab="gene mean t-stat",main="buccal swab") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
dev.off()
## png 
##   2
par(mfrow=c(1,1))

gmea_boxplot(bu,sets)

Comparison of blood and guthrie card

Now compare gene methylation scores for blood and guthrie cards. The mean limma t-statistic is used for a scatterplot, then rank-rank plot.

m1 <- merge(blm,gum,by=0)
rownames(m1) <- m1$Row.names
m1$Row.names=NULL
colnames(m1) <- c("bl","gu")

plot(m1$bl,m1$gu,xlab="blood at assessment",ylab="neonatal Guthrie card",
  col = rgb(red = 0.6, green = 0.6, blue = 0.6, alpha = 0.5) ,
  pch=19,cex=0.9)
grid()
abline(v=0,lty=1,lwd=3,col="black") ; abline(h=0,lty=1,lwd=3,col="black")
mtext("aggregated t-statistic for each gene (mean)")

rm1 <- apply(m1,2,rank)
mydiff <- apply(m1,2,function(x) { length(which(x<0)) } )
rm2 <- rm1
rm2[,1] <- rm2[,1] - mydiff[1]
rm2[,2] <- rm2[,2] - mydiff[2]
rnk <- rm2

palette <- colorRampPalette(c("white", "yellow", "orange", "red","darkred", "black"))
xmin = min(rnk[, 1])
xmax = max(rnk[, 1])
ymin = min(rnk[, 2])
ymax = max(rnk[, 2])
k <- MASS::kde2d(rnk[, 1], rnk[, 2])
X_AXIS = "blood at assessment"
Y_AXIS = "neonatal Guthrie card"

filled.contour(k, xlim = c(xmin, xmax), ylim = c(ymin, ymax),
  color.palette = palette, plot.title = {
  abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
  title(main="gene rank",xlab = X_AXIS, ylab = Y_AXIS)
})

pdf("genedens_ados.pdf")
par(mfrow=c(2,2))
plot(1)
plot(1)

plot(m1$bl,m1$gu,xlab="blood at assessment",ylab="neonatal Guthrie card",
  col = rgb(red = 0.6, green = 0.6, blue = 0.6, alpha = 0.5) ,
  pch=19,cex=0.9)
grid()
abline(v=0, h=0, lty=2, lwd=2, col="blue")
mtext("aggregated t-statistic for each gene (mean)")

filled.contour(k, xlim = c(xmin, xmax), ylim = c(ymin, ymax),
  color.palette = palette, plot.title = {
  abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
  title(main = "gene rank", xlab = X_AXIS, ylab = Y_AXIS)
})

dev.off()
## png 
##   2

Now we need a table of genes with top scores.

m3 <- m1
m3$med <- apply(rnk,1,median)
mg <- merge(gu$df,bl$df,by=0)
rownames(mg) <- mg$Row.names
mg <- mg[rownames(mg) %in% rownames(m3),]
dim(mg)
## [1] 25716    13
mg$med <- m3$med
mg <- mg[order(mg$med),]
mg <-subset(mg,nprobes.x>=5)
mg <-subset(mg,nprobes.x>=5)
dim(mg)
## [1] 21872    14
mg$Row.names = mg$sig.x = mg$sig.y = mg$FDR.x = mg$FDR.y = mg$median.x = mg$median.y = mg$nprobes.y = NULL

head(mg,20) %>%
  kbl(caption="Genes with coordinated hypomethylation in neonatal Guthrie card and blood at assessment") %>%
  kable_paper("hover", full_width = F)
Genes with coordinated hypomethylation in neonatal Guthrie card and blood at assessment
nprobes.x mean.x P.Value.x mean.y P.Value.y med
MIR1251 6 -1.4255875 0.0000323 -1.749671 0.0003086 -19468.5
BZW1L1 5 -1.2532992 0.0060798 -1.754408 0.0025244 -19435.0
LINC00871 7 -0.9595183 0.0061097 -1.713879 0.0001170 -19234.0
TAAR9 5 -1.0688897 0.0374636 -1.482612 0.0011922 -19206.5
LOC101926944 5 -0.9957280 0.0045594 -1.576021 0.0159759 -19203.0
NAALADL2-AS2 7 -1.0180771 0.0103643 -1.538871 0.0070632 -19200.0
OR52E8 8 -1.2018857 0.0004812 -1.384415 0.0189218 -19194.5
MIR340 5 -1.0338178 0.0274113 -1.491043 0.0038681 -19184.5
KERA 5 -0.8913769 0.1288030 -1.703936 0.0100214 -19146.0
LOC400655 9 -0.9597790 0.0000456 -1.527574 0.0191155 -19140.5
UGT2B4 5 -0.8995623 0.0232803 -1.651494 0.0002245 -19131.0
MIR490 5 -0.9209369 0.0340287 -1.566863 0.0087053 -19116.5
RHAG 8 -0.8685327 0.0006638 -1.706321 0.0000769 -19115.5
LINC00644 5 -1.0006948 0.0277128 -1.430676 0.0406665 -19108.5
SELE 9 -1.1272002 0.0383219 -1.289182 0.0007408 -19084.0
TMPRSS11F 11 -0.8876613 0.0001080 -1.551513 0.0007112 -19064.5
HISPPD1 5 -0.8494605 0.0259018 -1.618099 0.0001746 -19041.0
CRYGC 9 -0.9482129 0.0003966 -1.392217 0.0000031 -19020.5
FABP2 5 -0.7913966 0.0215867 -1.895935 0.0069802 -19020.5
PTH 5 -0.8111156 0.0339078 -1.703538 0.0184793 -19012.0

Gene set level enrichment

genesets <- gmt_import("ReactomePathways.gmt")

sfari <- readLines("sfari.txt")
lsfari <- list(sfari)
names(lsfari) <- "sfari"

cores = 8

Mitch enrichment analysis

Here I’m using the mean t-statistic for downstream enrichment. These have already been saved to objects called blm, gum and bum. Pathways are from REACTOME and analysis using mitch.

mitch_bl <- mitch_calc(blm,genesets,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
sig <- subset(mitch_bl$enrichment_result,`p.adjustANOVA`<0.01)

head(sig[order(-abs(sig$s.dist)),],30) %>%
  kbl(caption="BLOOD: Top effect size pathways found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)
BLOOD: Top effect size pathways found with mitch after 1% FDR filtering
set setSize pANOVA s.dist p.adjustANOVA
384 Eicosanoids 12 0.0001518 0.6314443 0.0015531
1194 SARS-CoV-1 modulates host translation machinery 34 0.0000001 0.5217548 0.0000052
216 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000001 -0.4893724 0.0000047
228 Competing endogenous RNAs (ceRNAs) regulate PTEN translation 16 0.0008579 -0.4812570 0.0078648
927 Peptide chain elongation 84 0.0000000 0.4663117 0.0000000
1246 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000003 -0.4646788 0.0000079
401 Eukaryotic Translation Termination 87 0.0000000 0.4644754 0.0000000
399 Eukaryotic Translation Elongation 88 0.0000000 0.4582574 0.0000000
1239 Selenocysteine synthesis 87 0.0000000 0.4487558 0.0000000
1534 Viral mRNA Translation 84 0.0000000 0.4451212 0.0000000
108 Autodegradation of the E3 ubiquitin ligase COP1 50 0.0000001 0.4331659 0.0000046
851 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 89 0.0000000 0.4277184 0.0000000
1176 Response of EIF2AK4 (GCN2) to amino acid deficiency 95 0.0000000 0.4275271 0.0000000
1519 Ubiquitin-dependent degradation of Cyclin D 50 0.0000002 0.4252715 0.0000065
1518 Ubiquitin Mediated Degradation of Phosphorylated Cdc25A 50 0.0000002 0.4252061 0.0000065
1573 p53-Independent DNA Damage Response 50 0.0000002 0.4252061 0.0000065
1574 p53-Independent G1/S DNA damage checkpoint 50 0.0000002 0.4252061 0.0000065
246 Cross-presentation of soluble exogenous antigens (endosomes) 46 0.0000007 0.4233906 0.0000154
454 Formation of a pool of free 40S subunits 94 0.0000000 0.4230774 0.0000000
21 AUF1 (hnRNP D0) binds and destabilizes mRNA 53 0.0000001 0.4225398 0.0000041
1218 SRP-dependent cotranslational protein targeting to membrane 105 0.0000000 0.4187931 0.0000000
1203 SCF(Skp2)-mediated degradation of p27/p21 58 0.0000001 0.4125927 0.0000024
407 FBXL7 down-regulates AURKA during mitotic entry and in early mitosis 53 0.0000002 0.4119269 0.0000065
1142 Regulation of activated PAK-2p34 by proteasome mediated degradation 48 0.0000008 0.4117493 0.0000167
504 GSK3B and BTRC:CUL1-mediated-degradation of NFE2L2 50 0.0000005 0.4115645 0.0000122
1357 Stabilization of p53 55 0.0000001 0.4098243 0.0000053
579 Hh mutants are degraded by ERAD 54 0.0000002 0.4081817 0.0000065
505 GTP hydrolysis and joining of the 60S ribosomal subunit 104 0.0000000 0.4065254 0.0000000
1540 Vpu mediated degradation of CD4 50 0.0000007 0.4053813 0.0000159
679 L13a-mediated translational silencing of Ceruloplasmin expression 103 0.0000000 0.4052466 0.0000000 
mitch_gu <- mitch_calc(gum,genesets,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
sig <- subset(mitch_gu$enrichment_result,`p.adjustANOVA`<0.01)

head(sig[order(-abs(sig$s.dist)),],30) %>%
  kbl(caption="GUTHRIE: Top effect size pathways found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)
GUTHRIE: Top effect size pathways found with mitch after 1% FDR filtering
set setSize pANOVA s.dist p.adjustANOVA
216 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000004 -0.4676449 0.0000122
904 PINK1-PRKN Mediated Mitophagy 21 0.0003534 0.4502581 0.0028893
13 APC/C:Cdc20 mediated degradation of Cyclin B 24 0.0001441 0.4481226 0.0014374
773 Mitophagy 27 0.0001137 0.4290077 0.0012108
96 Aspirin ADME 44 0.0000009 -0.4283638 0.0000228
947 Phosphorylation of the APC/C 20 0.0013935 0.4127413 0.0083714
403 Expression and translocation of olfactory receptors 356 0.0000000 -0.4043606 0.0000000
11 APC-Cdc20 mediated degradation of Nek2A 26 0.0004031 0.4007755 0.0031652
875 Olfactory Signaling Pathway 363 0.0000000 -0.4003579 0.0000000
1379 Synthesis of active ubiquitin: roles of E1 and E2 enzymes 27 0.0003176 0.4002390 0.0026717
25 Aberrant regulation of mitotic cell cycle due to RB1 defects 36 0.0000330 0.3997404 0.0004589
528 Glucuronidation 25 0.0005866 -0.3971243 0.0043069
1475 Transcriptional regulation of granulopoiesis 31 0.0002269 0.3825190 0.0020451
338 Diseases of mitotic cell cycle 37 0.0000601 0.3810843 0.0007329
1193 SARS-CoV-1 activates/modulates innate immune responses 36 0.0000841 0.3786107 0.0009601
1247 Sensory perception of taste 47 0.0000078 -0.3767662 0.0001361
1244 Senescence-Associated Secretory Phenotype (SASP) 53 0.0000031 0.3702452 0.0000614
107 Autodegradation of Cdh1 by Cdh1:APC/C 62 0.0000005 0.3700687 0.0000128
1246 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000417 -0.3697219 0.0005563
1203 SCF(Skp2)-mediated degradation of p27/p21 58 0.0000023 0.3588126 0.0000478
504 GSK3B and BTRC:CUL1-mediated-degradation of NFE2L2 50 0.0000143 0.3545438 0.0002252
793 N-glycan trimming in the ER and Calnexin/Calreticulin cycle 34 0.0003527 0.3540062 0.0028893
339 Diseases of programmed cell death 43 0.0000781 0.3480437 0.0009109
1173 Respiratory electron transport 90 0.0000000 0.3474232 0.0000006
102 Assembly of the pre-replicative complex 82 0.0000001 0.3450819 0.0000023
862 Nuclear import of Rev protein 32 0.0007631 0.3437184 0.0052620
699 Lysosome Vesicle Biogenesis 33 0.0006818 0.3415878 0.0048491
1104 Regulation of APC/C activators between G1/S and early anaphase 79 0.0000002 0.3396609 0.0000055
143 CDK-mediated phosphorylation and removal of Cdc6 71 0.0000010 0.3360274 0.0000247
224 Collagen chain trimerization 42 0.0001738 -0.3347096 0.0016603 
mitch_bu <- mitch_calc(bum,genesets,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
sig <- subset(mitch_bu$enrichment_result,`p.adjustANOVA`<0.01)

head(sig[order(-abs(sig$s.dist)),],30) %>%
  kbl(caption="BUCCAL: Top effect size pathways found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)
BUCCAL: Top effect size pathways found with mitch after 1% FDR filtering
set setSize pANOVA s.dist p.adjustANOVA
1194 SARS-CoV-1 modulates host translation machinery 34 2.00e-07 0.5177226 0.0000098
1418 TRAF6 mediated IRF7 activation 25 5.09e-05 -0.4679910 0.0017180
216 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 6.00e-07 -0.4614606 0.0000333
403 Expression and translocation of olfactory receptors 358 0.00e+00 -0.4178164 0.0000000
927 Peptide chain elongation 84 0.00e+00 0.4117980 0.0000000
875 Olfactory Signaling Pathway 365 0.00e+00 -0.4085329 0.0000000
1534 Viral mRNA Translation 84 0.00e+00 0.4052533 0.0000000
401 Eukaryotic Translation Termination 87 0.00e+00 0.3980434 0.0000000
399 Eukaryotic Translation Elongation 88 0.00e+00 0.3949691 0.0000000
1239 Selenocysteine synthesis 87 0.00e+00 0.3784774 0.0000001
851 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 89 0.00e+00 0.3728020 0.0000001
1176 Response of EIF2AK4 (GCN2) to amino acid deficiency 95 0.00e+00 0.3609384 0.0000001
454 Formation of a pool of free 40S subunits 94 0.00e+00 0.3592418 0.0000002
1246 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 7.68e-05 -0.3567897 0.0024354
466 Formation of the ternary complex, and subsequently, the 43S complex 46 4.47e-05 0.3477262 0.0015759
1200 SARS-CoV-2 modulates host translation machinery 46 6.31e-05 0.3408465 0.0020425
1238 Selenoamino acid metabolism 103 0.00e+00 0.3372552 0.0000003
679 L13a-mediated translational silencing of Ceruloplasmin expression 103 0.00e+00 0.3352722 0.0000003
505 GTP hydrolysis and joining of the 60S ribosomal subunit 104 0.00e+00 0.3338761 0.0000003
1187 Ribosomal scanning and start codon recognition 53 2.82e-05 0.3323725 0.0010922
54 Activation of the mRNA upon binding of the cap-binding complex and eIFs, and subsequent binding to 43S 54 2.39e-05 0.3322706 0.0009721
1481 Translation initiation complex formation 53 3.12e-05 0.3305846 0.0011714
167 Cap-dependent Translation Initiation 111 0.00e+00 0.3255501 0.0000003
400 Eukaryotic Translation Initiation 111 0.00e+00 0.3255501 0.0000003
850 Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC) 108 0.00e+00 0.3241781 0.0000004
852 Nonsense-Mediated Decay (NMD) 108 0.00e+00 0.3241781 0.0000004
720 Major pathway of rRNA processing in the nucleolus and cytosol 171 0.00e+00 0.3154260 0.0000000
1218 SRP-dependent cotranslational protein targeting to membrane 105 1.00e-07 0.2992844 0.0000067
1580 rRNA processing in the nucleus and cytosol 180 0.00e+00 0.2985439 0.0000000
194 Cellular response to starvation 147 0.00e+00 0.2901196 0.0000001 
mitch_plots(res=mitch_bu,outfile="buccal_mitch_ados.pdf")
## png 
##   2
mitch_report(res=mitch_bu,outfile="buccal_mitch_ados.html")
## Dataset saved as " /tmp/RtmpupMn4p/buccal_mitch_ados.rds ".
## 
## 
## processing file: mitch.Rmd
## 1/34                             
## 2/34 [checklibraries]            
## 3/34                             
## 4/34 [peek]                      
## 5/34                             
## 6/34 [metrics]                   
## 7/34                             
## 8/34 [scatterplot]
## 9/34                             
## 10/34 [contourplot]               
## 11/34                             
## 12/34 [input_geneset_metrics1]    
## 13/34                             
## 14/34 [input_geneset_metrics2]
## 15/34                             
## 16/34 [input_geneset_metrics3]
## 17/34                             
## 18/34 [echart1d]                  
## 19/34 [echart2d]                  
## 20/34                             
## 21/34 [heatmap]                   
## 22/34                             
## 23/34 [effectsize]                
## 24/34                             
## 25/34 [results_table]             
## 26/34                             
## 27/34 [results_table_complete]    
## 28/34                             
## 29/34 [detailed_geneset_reports1d]
## 30/34                             
## 31/34 [detailed_geneset_reports2d]
## 32/34                             
## 33/34 [session_info]              
## 34/34
## output file: /home/mdz/projects/asd_meth/mitch.knit.md
## /usr/bin/pandoc +RTS -K512m -RTS /home/mdz/projects/asd_meth/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpupMn4p/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /tmp/RtmpupMn4p/rmarkdown-str9bd7a7f09c9c0.html
## 
## Output created: /tmp/RtmpupMn4p/mitch_report.html
## [1] TRUE
subset(mitch_bu$enrichment_result,p.adjustANOVA<0.01 & s.dist<0) %>%
  kbl(caption="Reactome pathways with an inverse methylation association with ADOS in buccal swab samples.") %>%
  kable_paper("hover", full_width = F)
Reactome pathways with an inverse methylation association with ADOS in buccal swab samples.
set setSize pANOVA s.dist p.adjustANOVA
1418 TRAF6 mediated IRF7 activation 25 0.0000509 -0.4679910 0.0017180
216 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000006 -0.4614606 0.0000333
403 Expression and translocation of olfactory receptors 358 0.0000000 -0.4178164 0.0000000
875 Olfactory Signaling Pathway 365 0.0000000 -0.4085329 0.0000000
1246 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000768 -0.3567897 0.0024354
200 Chemokine receptors bind chemokines 56 0.0001806 -0.2892053 0.0056149
1245 Sensory Perception 571 0.0000000 -0.2683827 0.0000000
675 Keratinization 213 0.0000266 -0.1668486 0.0010550
471 G alpha (i) signalling events 305 0.0000624 -0.1331282 0.0020425

Mitch joint analysis of blood and guthrie card. As using the median value is apparently the most sensitive, we’ll go with that.

m <- merge(blm,gum,by=0)
rownames(m) <- m$Row.names
m$Row.names = NULL
colnames(m) <- c("bl","gu")
m[is.na(m)] <- 0
res <- mitch_calc(m, genesets, priority="effect")
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
mitch_report(res,"blgu_mitch.html",overwrite=TRUE)
## Note: overwriting existing report
## Dataset saved as " /tmp/RtmpupMn4p/blgu_mitch.rds ".
## 
## 
## processing file: mitch.Rmd
## 1/34                             
## 2/34 [checklibraries]            
## 3/34                             
## 4/34 [peek]                      
## 5/34                             
## 6/34 [metrics]                   
## 7/34                             
## 8/34 [scatterplot]
## 9/34                             
## 10/34 [contourplot]
## 11/34                             
## 12/34 [input_geneset_metrics1]    
## 13/34                             
## 14/34 [input_geneset_metrics2]
## 15/34                             
## 16/34 [input_geneset_metrics3]
## 17/34                             
## 18/34 [echart1d]                  
## 19/34 [echart2d]                  
## 20/34                             
## 21/34 [heatmap]
## 22/34                             
## 23/34 [effectsize]                
## 24/34                             
## 25/34 [results_table]             
## 26/34                             
## 27/34 [results_table_complete]    
## 28/34                             
## 29/34 [detailed_geneset_reports1d]
## 30/34                             
## 31/34 [detailed_geneset_reports2d]
## 32/34                             
## 33/34 [session_info]              
## 34/34
## output file: /home/mdz/projects/asd_meth/mitch.knit.md
## /usr/bin/pandoc +RTS -K512m -RTS /home/mdz/projects/asd_meth/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpupMn4p/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /tmp/RtmpupMn4p/rmarkdown-str9bd7a792c17b9.html
## 
## Output created: /tmp/RtmpupMn4p/mitch_report.html
## [1] TRUE
mitch_plots(res=res,outfile="blgu_mitch_ados.pdf")
## png 
##   2
sig <- subset(res$enrichment_result,p.adjustMANOVA<0.01)

head(sig[order(-abs(sig$s.dist)),],30) %>%
  kbl(caption="Top pathways found with mitch (effect size after 1% FDR filtering)") %>%
  kable_paper("hover", full_width = F)
Top pathways found with mitch (effect size after 1% FDR filtering)
set setSize pMANOVA s.bl s.gu p.bl p.gu s.dist SD p.adjustMANOVA
216 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000000 -0.4901753 -0.4676449 0.0000001 0.0000004 0.6774685 0.0159314 0.0000000
1399 TICAM1-dependent activation of IRF3/IRF7 12 0.0012354 0.4502023 0.4919727 0.0069205 0.0031649 0.6668728 0.0295361 0.0076240
384 Eicosanoids 12 0.0002957 0.6313349 -0.1018713 0.0001522 0.5411614 0.6395010 0.5184551 0.0021605
1194 SARS-CoV-1 modulates host translation machinery 34 0.0000001 0.5215509 0.3046011 0.0000001 0.0021091 0.6039844 0.1534067 0.0000015
1246 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000000 -0.4654969 -0.3697219 0.0000002 0.0000417 0.5944592 0.0677232 0.0000002
228 Competing endogenous RNAs (ceRNAs) regulate PTEN translation 16 0.0001272 -0.4822519 0.2740224 0.0008369 0.0577137 0.5546667 0.5347667 0.0010187
528 Glucuronidation 25 0.0000806 -0.3782663 -0.3971243 0.0010586 0.0005866 0.5484461 0.0133346 0.0006625
1203 SCF(Skp2)-mediated degradation of p27/p21 58 0.0000000 0.4123054 0.3588126 0.0000001 0.0000023 0.5465731 0.0378251 0.0000000
504 GSK3B and BTRC:CUL1-mediated-degradation of NFE2L2 50 0.0000000 0.4112772 0.3545438 0.0000005 0.0000143 0.5430011 0.0401166 0.0000002
13 APC/C:Cdc20 mediated degradation of Cyclin B 24 0.0001346 0.3003821 0.4481226 0.0108444 0.0001441 0.5394843 0.1044683 0.0010727
927 Peptide chain elongation 84 0.0000000 0.4660659 0.2690810 0.0000000 0.0000200 0.5381654 0.1392894 0.0000000
904 PINK1-PRKN Mediated Mitophagy 21 0.0004631 0.2872506 0.4502581 0.0226638 0.0003534 0.5340836 0.1152637 0.0032216
401 Eukaryotic Translation Termination 87 0.0000000 0.4642281 0.2632838 0.0000000 0.0000217 0.5336910 0.1420890 0.0000000
1518 Ubiquitin Mediated Degradation of Phosphorylated Cdc25A 50 0.0000000 0.4249295 0.3125505 0.0000002 0.0001310 0.5274968 0.0794640 0.0000004
1573 p53-Independent DNA Damage Response 50 0.0000000 0.4249295 0.3125505 0.0000002 0.0001310 0.5274968 0.0794640 0.0000004
1574 p53-Independent G1/S DNA damage checkpoint 50 0.0000000 0.4249295 0.3125505 0.0000002 0.0001310 0.5274968 0.0794640 0.0000004
1247 Sensory perception of taste 47 0.0000001 -0.3649348 -0.3767662 0.0000149 0.0000078 0.5245285 0.0083661 0.0000013
107 Autodegradation of Cdh1 by Cdh1:APC/C 62 0.0000000 0.3716103 0.3700687 0.0000004 0.0000005 0.5244474 0.0010901 0.0000000
579 Hh mutants are degraded by ERAD 54 0.0000000 0.4078909 0.3275877 0.0000002 0.0000310 0.5231526 0.0567829 0.0000002
399 Eukaryotic Translation Elongation 88 0.0000000 0.4580045 0.2508669 0.0000000 0.0000471 0.5222090 0.1464684 0.0000000
108 Autodegradation of the E3 ubiquitin ligase COP1 50 0.0000000 0.4328964 0.2891483 0.0000001 0.0004030 0.5205824 0.1016453 0.0000005
1519 Ubiquitin-dependent degradation of Cyclin D 50 0.0000000 0.4249949 0.2999003 0.0000002 0.0002428 0.5201546 0.0884553 0.0000005
1534 Viral mRNA Translation 84 0.0000000 0.4448583 0.2671498 0.0000000 0.0000229 0.5189103 0.1256589 0.0000000
1239 Selenocysteine synthesis 87 0.0000000 0.4484943 0.2576795 0.0000000 0.0000323 0.5172484 0.1349264 0.0000000
1540 Vpu mediated degradation of CD4 50 0.0000001 0.4050884 0.3180831 0.0000007 0.0000993 0.5150471 0.0615221 0.0000008
1357 Stabilization of p53 55 0.0000000 0.4095348 0.3100843 0.0000001 0.0000691 0.5136838 0.0703221 0.0000002
21 AUF1 (hnRNP D0) binds and destabilizes mRNA 53 0.0000000 0.4222620 0.2887859 0.0000001 0.0002748 0.5115687 0.0943818 0.0000004
829 Negative regulation of NOTCH4 signaling 52 0.0000000 0.3942742 0.3241131 0.0000009 0.0000524 0.5103934 0.0496114 0.0000006
1204 SCF-beta-TrCP mediated degradation of Emi1 53 0.0000000 0.3917813 0.3249161 0.0000008 0.0000425 0.5089823 0.0472809 0.0000005
1142 Regulation of activated PAK-2p34 by proteasome mediated degradation 48 0.0000002 0.4114624 0.2964797 0.0000008 0.0003784 0.5071504 0.0813050 0.0000020

Now look at the pathways with hypomethylation.

hypo <- subset(res$enrichment_result,p.adjustMANOVA<0.01 & s.bl <0 & s.gu < 0 )
nrow(hypo)
## [1] 15
hypo <- hypo[nrow(hypo):1,]

par(mar=c(4.1, 25.1, 5.1, 1.1))

barplot(t(abs(hypo[,c("s.bl","s.gu")])),
  beside=TRUE,horiz=TRUE,las=1,names.arg=hypo$set,
  col=gray.colors(2),
  xlab="enrichment score")

legend("topright", inset=0, title="tissue",
  c("blood at assessment","neonatal Guthrie card"),
  fill=gray.colors(2), horiz=FALSE, cex=1)

pdf("adosbarplot.pdf")

par(mar=c(4.1, 25.1, 5.1, 1.1))

barplot(t(abs(hypo[,c("s.bl","s.gu")])),
  beside=TRUE,horiz=TRUE,las=1,names.arg=hypo$set, 
  col=gray.colors(2),
  xlab="enrichment score")

legend("topright", inset=0, title="tissue",
  c("blood at assessment","neonatal Guthrie card"), 
  fill=gray.colors(2), horiz=FALSE, cex=1)

dev.off()
## png 
##   2
mar = c(5.1, 4.1, 4.1, 2.1)

Now look at metatropic glutamate/pheromone receptors.

phero <- genesets[grep("pheromone",names(genesets))]
phero
## $`Class C/3 (Metabotropic glutamate/pheromone receptors)`
##  [1] "CASR"    "GABBR1"  "GABBR2"  "GPRC6A"  "GRM1"    "GRM2"    "GRM3"   
##  [8] "GRM4"    "GRM5"    "GRM6"    "GRM7"    "GRM8"    "TAS1R1"  "TAS1R2" 
## [15] "TAS1R3"  "TAS2R1"  "TAS2R10" "TAS2R13" "TAS2R14" "TAS2R16" "TAS2R19"
## [22] "TAS2R20" "TAS2R3"  "TAS2R30" "TAS2R31" "TAS2R38" "TAS2R39" "TAS2R4" 
## [29] "TAS2R40" "TAS2R41" "TAS2R42" "TAS2R43" "TAS2R45" "TAS2R46" "TAS2R5" 
## [36] "TAS2R50" "TAS2R60" "TAS2R7"  "TAS2R8"  "TAS2R9"
rnk2 <- rnk[which(rownames(rnk) %in% phero[[1]]),]

plot(rnk2,main=names(phero),xlim=c(-21569,4146),ylim=c(-17678,8037),
  xlab="Rank in blood at assessment",
  ylab="Rank in neonatal Guthrie card",
  cex=1.5,pch=19,col="lightgray")

text(rnk2,labels=rownames(rnk2))

abline(v=0,h=0,lwd=2,lty=2)

pdf("phero.pdf")

plot(rnk2,main=names(phero),xlim=c(-21569,4146),ylim=c(-17678,8037),
  xlab="Rank in blood at assessment",
  ylab="Rank in neonatal Guthrie card",
  cex=1.5,pch=19,col="lightgray")

text(rnk2,labels=rownames(rnk2))

abline(v=0,h=0,lwd=2,lty=2)

dev.off()
## png 
##   2

Make a heatmap of top TAS2R gene probes.

colcolfunc <- colorRampPalette(c("yellow", "orange"))

colcols <- colcolfunc(20)[bl_design[,"ADOS"]]

phero_probes <- unlist(sets[which(names(sets) %in% phero[[1]])])

phero_mvals <- bl_mvals[rownames(bl_mvals) %in% phero_probes,]

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(phero_mvals,trace="none",col=colfunc(25), margins = c(5,8), scale="row",
  main="Class C/3 (Metabotropic glutamate/pheromone receptors) probes",
  ColSideColors = colcols )

Now let’s take a look at the SFARI gene list.

str(lsfari)
## List of 1
##  $ sfari: chr [1:1152] "ABAT" "ABCA10" "ABCA13" "ABCA7" ...
sfarires <- mitch_calc(m1, lsfari, priority="effect")
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
t(sfarires$enrichment_result)
##                1             
## set            "sfari"       
## setSize        "1034"        
## pMANOVA        "7.495978e-06"
## s.bl           "-0.08355524" 
## s.gu           "-0.04650798" 
## p.bl           "5.119228e-06"
## p.gu           "0.01115786"  
## s.dist         "0.09562673"  
## SD             "0.02619637"  
## p.adjustMANOVA "7.495978e-06"
lrand <-lapply(1:1034,function(i) { set.seed(i*1000) ; sample(rownames(m1),1000) })
names(lrand) <- paste("random",1:1034)
lsets <- c(lsfari,lrand)

sfarires2 <- mitch_calc(m1, lsets, priority="effect",cores=16,resrows=10)
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
head(sfarires2$enrichment_result) %>%
  kbl(caption = "Enrichment analysis of SFARI genes") %>%
  kable_paper("hover", full_width = F)
Enrichment analysis of SFARI genes
set setSize pMANOVA s.bl s.gu p.bl p.gu s.dist SD p.adjustMANOVA
1 sfari 1034 0.0000075 -0.0835552 -0.0465080 0.0000051 0.0111579 0.0956267 0.0261964 0.0077583
380 random 379 1000 0.0009309 0.0685287 0.0016664 0.0002332 0.9287050 0.0685490 0.0472788 0.3211440
620 random 619 1000 0.0003263 0.0332501 -0.0589756 0.0741981 0.0015405 0.0677030 0.0652134 0.1688350
68 random 67 1000 0.0034017 0.0251146 0.0613512 0.1774881 0.0009858 0.0662926 0.0256232 0.6312853
886 random 885 1000 0.0044702 0.0315138 0.0576910 0.0906163 0.0019486 0.0657371 0.0185100 0.6312853
77 random 76 1000 0.0029505 0.0635025 0.0096897 0.0006494 0.6028683 0.0642375 0.0380514 0.6312853 
mitch_plots(res=sfarires2,outfile="sfari_mitch_ados.pdf")
## png 
##   2
mitch_report(res=sfarires2,outfile="sfari_mitch_ados.html",overwrite=TRUE)
## Note: overwriting existing report
## Dataset saved as " /tmp/RtmpupMn4p/sfari_mitch_ados.rds ".
## 
## 
## processing file: mitch.Rmd
## 1/34                             
## 2/34 [checklibraries]            
## 3/34                             
## 4/34 [peek]                      
## 5/34                             
## 6/34 [metrics]                   
## 7/34                             
## 8/34 [scatterplot]
## 9/34                             
## 10/34 [contourplot]
## 11/34                             
## 12/34 [input_geneset_metrics1]    
## 13/34                             
## 14/34 [input_geneset_metrics2]
## 15/34                             
## 16/34 [input_geneset_metrics3]
## 17/34                             
## 18/34 [echart1d]                  
## 19/34 [echart2d]                  
## 20/34                             
## 21/34 [heatmap]
## 22/34                             
## 23/34 [effectsize]                
## 24/34                             
## 25/34 [results_table]             
## 26/34                             
## 27/34 [results_table_complete]    
## 28/34                             
## 29/34 [detailed_geneset_reports1d]
## 30/34                             
## 31/34 [detailed_geneset_reports2d]
## 32/34                             
## 33/34 [session_info]              
## 34/34
## output file: /home/mdz/projects/asd_meth/mitch.knit.md
## /usr/bin/pandoc +RTS -K512m -RTS /home/mdz/projects/asd_meth/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpupMn4p/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /tmp/RtmpupMn4p/rmarkdown-str9bd7a3f3aed9b.html
## 
## Output created: /tmp/RtmpupMn4p/mitch_report.html
## [1] TRUE

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] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] pkgload_1.3.3                                      
##  [2] GGally_2.1.2                                       
##  [3] ggplot2_3.4.4                                      
##  [4] reshape2_1.4.4                                     
##  [5] beeswarm_0.4.0                                     
##  [6] gplots_3.1.3                                       
##  [7] gtools_3.9.4                                       
##  [8] tibble_3.2.1                                       
##  [9] dplyr_1.1.3                                        
## [10] echarts4r_0.4.5                                    
## [11] eulerr_7.0.0                                       
## [12] kableExtra_1.3.4                                   
## [13] data.table_1.14.8                                  
## [14] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [15] minfi_1.48.0                                       
## [16] bumphunter_1.44.0                                  
## [17] locfit_1.5-9.8                                     
## [18] iterators_1.0.14                                   
## [19] foreach_1.5.2                                      
## [20] Biostrings_2.70.1                                  
## [21] XVector_0.42.0                                     
## [22] SummarizedExperiment_1.32.0                        
## [23] Biobase_2.62.0                                     
## [24] MatrixGenerics_1.14.0                              
## [25] matrixStats_1.1.0                                  
## [26] GenomicRanges_1.54.1                               
## [27] GenomeInfoDb_1.38.1                                
## [28] IRanges_2.36.0                                     
## [29] S4Vectors_0.40.1                                   
## [30] BiocGenerics_0.48.1                                
## [31] mitch_1.14.0                                       
## [32] limma_3.58.1                                       
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.2             later_1.3.1              
##   [3] BiocIO_1.12.0             bitops_1.0-7             
##   [5] filelock_1.0.2            preprocessCore_1.64.0    
##   [7] XML_3.99-0.15             lifecycle_1.0.4          
##   [9] lattice_0.22-5            MASS_7.3-60              
##  [11] base64_2.0.1              scrime_1.3.5             
##  [13] magrittr_2.0.3            sass_0.4.7               
##  [15] rmarkdown_2.25            jquerylib_0.1.4          
##  [17] yaml_2.3.7                httpuv_1.6.12            
##  [19] doRNG_1.8.6               askpass_1.2.0            
##  [21] DBI_1.1.3                 RColorBrewer_1.1-3       
##  [23] abind_1.4-5               zlibbioc_1.48.0          
##  [25] rvest_1.0.3               quadprog_1.5-8           
##  [27] purrr_1.0.2               RCurl_1.98-1.13          
##  [29] rappdirs_0.3.3            GenomeInfoDbData_1.2.11  
##  [31] genefilter_1.84.0         annotate_1.80.0          
##  [33] svglite_2.1.2             DelayedMatrixStats_1.24.0
##  [35] codetools_0.2-19          DelayedArray_0.28.0      
##  [37] xml2_1.3.5                tidyselect_1.2.0         
##  [39] farver_2.1.1              beanplot_1.3.1           
##  [41] BiocFileCache_2.10.1      illuminaio_0.44.0        
##  [43] webshot_0.5.5             GenomicAlignments_1.38.0 
##  [45] jsonlite_1.8.7            multtest_2.58.0          
##  [47] ellipsis_0.3.2            survival_3.5-7           
##  [49] systemfonts_1.0.5         tools_4.3.2              
##  [51] progress_1.2.2            Rcpp_1.0.11              
##  [53] glue_1.6.2                gridExtra_2.3            
##  [55] SparseArray_1.2.2         xfun_0.41                
##  [57] HDF5Array_1.30.0          withr_2.5.2              
##  [59] fastmap_1.1.1             rhdf5filters_1.14.1      
##  [61] fansi_1.0.5               openssl_2.1.1            
##  [63] caTools_1.18.2            digest_0.6.33            
##  [65] R6_2.5.1                  mime_0.12                
##  [67] colorspace_2.1-0          biomaRt_2.58.0           
##  [69] RSQLite_2.3.3             tidyr_1.3.0              
##  [71] utf8_1.2.4                generics_0.1.3           
##  [73] rtracklayer_1.62.0        prettyunits_1.2.0        
##  [75] httr_1.4.7                htmlwidgets_1.6.2        
##  [77] S4Arrays_1.2.0            pkgconfig_2.0.3          
##  [79] gtable_0.3.4              blob_1.2.4               
##  [81] siggenes_1.76.0           htmltools_0.5.7          
##  [83] scales_1.2.1              png_0.1-8                
##  [85] knitr_1.45                rstudioapi_0.15.0        
##  [87] tzdb_0.4.0                rjson_0.2.21             
##  [89] nlme_3.1-163              curl_5.1.0               
##  [91] cachem_1.0.8              rhdf5_2.46.0             
##  [93] stringr_1.5.0             KernSmooth_2.23-22       
##  [95] AnnotationDbi_1.64.1      restfulr_0.0.15          
##  [97] GEOquery_2.70.0           pillar_1.9.0             
##  [99] grid_4.3.2                reshape_0.8.9            
## [101] vctrs_0.6.4               promises_1.2.1           
## [103] dbplyr_2.4.0              xtable_1.8-4             
## [105] evaluate_0.23             readr_2.1.4              
## [107] GenomicFeatures_1.54.1    cli_3.6.1                
## [109] compiler_4.3.2            Rsamtools_2.18.0         
## [111] rlang_1.1.2               crayon_1.5.2             
## [113] rngtools_1.5.2            labeling_0.4.3           
## [115] nor1mix_1.3-0             mclust_6.0.0             
## [117] plyr_1.8.9                stringi_1.7.12           
## [119] viridisLite_0.4.2         BiocParallel_1.36.0      
## [121] munsell_0.5.0             Matrix_1.6-2             
## [123] hms_1.1.3                 sparseMatrixStats_1.14.0 
## [125] bit64_4.0.5               Rhdf5lib_1.24.0          
## [127] KEGGREST_1.42.0           statmod_1.5.0            
## [129] shiny_1.7.5.1             highr_0.10               
## [131] memoise_2.0.1             bslib_0.5.1              
## [133] bit_4.0.5