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("RIdeogram")
  library("GenomicRanges")
  library("tictoc")
  library("globaltest")
  library("ebGSEA")
  data("dualmap850kEID")
})

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
380 Eicosanoids 12 0.0001518 0.6314443 0.0015541
1163 SARS-CoV-1 modulates host translation machinery 34 0.0000001 0.5217548 0.0000050
212 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000001 -0.4893724 0.0000046
224 Competing endogenous RNAs (ceRNAs) regulate PTEN translation 16 0.0008579 -0.4812570 0.0078946
903 Peptide chain elongation 84 0.0000000 0.4663117 0.0000000
1215 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000003 -0.4646788 0.0000077
397 Eukaryotic Translation Termination 87 0.0000000 0.4644754 0.0000000
395 Eukaryotic Translation Elongation 88 0.0000000 0.4582574 0.0000000
1208 Selenocysteine synthesis 87 0.0000000 0.4487558 0.0000000
1496 Viral mRNA Translation 84 0.0000000 0.4451212 0.0000000
107 Autodegradation of the E3 ubiquitin ligase COP1 50 0.0000001 0.4331659 0.0000044
829 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 89 0.0000000 0.4277184 0.0000000
1145 Response of EIF2AK4 (GCN2) to amino acid deficiency 95 0.0000000 0.4275271 0.0000000
1482 Ubiquitin-dependent degradation of Cyclin D 50 0.0000002 0.4252715 0.0000063
1481 Ubiquitin Mediated Degradation of Phosphorylated Cdc25A 50 0.0000002 0.4252061 0.0000063
1533 p53-Independent DNA Damage Response 50 0.0000002 0.4252061 0.0000063
1534 p53-Independent G1/S DNA damage checkpoint 50 0.0000002 0.4252061 0.0000063
242 Cross-presentation of soluble exogenous antigens (endosomes) 46 0.0000007 0.4233906 0.0000152
449 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.0000040
1187 SRP-dependent cotranslational protein targeting to membrane 105 0.0000000 0.4187931 0.0000000
1172 SCF(Skp2)-mediated degradation of p27/p21 58 0.0000001 0.4125927 0.0000023
403 FBXL7 down-regulates AURKA during mitotic entry and in early mitosis 53 0.0000002 0.4119269 0.0000064
1111 Regulation of activated PAK-2p34 by proteasome mediated degradation 48 0.0000008 0.4117493 0.0000165
492 GSK3B and BTRC:CUL1-mediated-degradation of NFE2L2 50 0.0000005 0.4115645 0.0000120
1323 Stabilization of p53 55 0.0000001 0.4098243 0.0000051
565 Hh mutants are degraded by ERAD 54 0.0000002 0.4081817 0.0000064
493 GTP hydrolysis and joining of the 60S ribosomal subunit 104 0.0000000 0.4065254 0.0000000
1502 Vpu mediated degradation of CD4 50 0.0000007 0.4053813 0.0000157
663 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
212 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000004 -0.4676449 0.0000116
882 PINK1-PRKN Mediated Mitophagy 21 0.0003534 0.4502581 0.0028311
13 APC/C:Cdc20 mediated degradation of Cyclin B 24 0.0001441 0.4481226 0.0014011
753 Mitophagy 27 0.0001137 0.4290077 0.0011647
95 Aspirin ADME 44 0.0000009 -0.4283638 0.0000217
923 Phosphorylation of the APC/C 20 0.0013935 0.4127413 0.0083826
399 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.0031162
853 Olfactory Signaling Pathway 363 0.0000000 -0.4003579 0.0000000
1345 Synthesis of active ubiquitin: roles of E1 and E2 enzymes 27 0.0003176 0.4002390 0.0026182
25 Aberrant regulation of mitotic cell cycle due to RB1 defects 36 0.0000330 0.3997404 0.0004322
516 Glucuronidation 25 0.0005866 -0.3971243 0.0041983
1162 SARS-CoV-1 activates/modulates innate immune responses 37 0.0000397 0.3902516 0.0004993
1440 Transcriptional regulation of granulopoiesis 31 0.0002269 0.3825190 0.0020164
334 Diseases of mitotic cell cycle 37 0.0000601 0.3810843 0.0006953
1216 Sensory perception of taste 47 0.0000078 -0.3767662 0.0001326
1213 Senescence-Associated Secretory Phenotype (SASP) 53 0.0000031 0.3702452 0.0000613
106 Autodegradation of Cdh1 by Cdh1:APC/C 62 0.0000005 0.3700687 0.0000123
1215 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000417 -0.3697219 0.0005204
1172 SCF(Skp2)-mediated degradation of p27/p21 58 0.0000023 0.3588126 0.0000479
492 GSK3B and BTRC:CUL1-mediated-degradation of NFE2L2 50 0.0000143 0.3545438 0.0002195
773 N-glycan trimming in the ER and Calnexin/Calreticulin cycle 34 0.0003527 0.3540062 0.0028311
1022 RIPK1-mediated regulated necrosis 27 0.0016249 0.3503229 0.0096992
1126 Regulation of necroptotic cell death 27 0.0016249 0.3503229 0.0096992
335 Diseases of programmed cell death 43 0.0000781 0.3480437 0.0008688
1142 Respiratory electron transport 90 0.0000000 0.3474232 0.0000006
101 Assembly of the pre-replicative complex 82 0.0000001 0.3450819 0.0000023
840 Nuclear import of Rev protein 32 0.0007631 0.3437184 0.0050850
683 Lysosome Vesicle Biogenesis 33 0.0006818 0.3415878 0.0047268
1079 Regulation of APC/C activators between G1/S and early anaphase 79 0.0000002 0.3396609 0.0000053 
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
1163 SARS-CoV-1 modulates host translation machinery 34 2.00e-07 0.5177226 0.0000095
1384 TRAF6 mediated IRF7 activation 26 3.55e-05 -0.4684102 0.0012593
212 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 6.00e-07 -0.4614606 0.0000324
399 Expression and translocation of olfactory receptors 358 0.00e+00 -0.4178164 0.0000000
903 Peptide chain elongation 84 0.00e+00 0.4117980 0.0000000
853 Olfactory Signaling Pathway 365 0.00e+00 -0.4085329 0.0000000
1496 Viral mRNA Translation 84 0.00e+00 0.4052533 0.0000000
397 Eukaryotic Translation Termination 87 0.00e+00 0.3980434 0.0000000
395 Eukaryotic Translation Elongation 88 0.00e+00 0.3949691 0.0000000
1208 Selenocysteine synthesis 87 0.00e+00 0.3784774 0.0000001
829 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 89 0.00e+00 0.3728020 0.0000001
1145 Response of EIF2AK4 (GCN2) to amino acid deficiency 95 0.00e+00 0.3609384 0.0000001
449 Formation of a pool of free 40S subunits 94 0.00e+00 0.3592418 0.0000002
1215 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 7.68e-05 -0.3567897 0.0024224
455 Formation of the ternary complex, and subsequently, the 43S complex 46 4.47e-05 0.3477262 0.0015028
1169 SARS-CoV-2 modulates host translation machinery 46 6.31e-05 0.3408465 0.0020325
1207 Selenoamino acid metabolism 103 0.00e+00 0.3372552 0.0000003
663 L13a-mediated translational silencing of Ceruloplasmin expression 103 0.00e+00 0.3352722 0.0000003
493 GTP hydrolysis and joining of the 60S ribosomal subunit 104 0.00e+00 0.3338761 0.0000003
1156 Ribosomal scanning and start codon recognition 53 2.82e-05 0.3323725 0.0011193
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.0009975
1446 Translation initiation complex formation 53 3.12e-05 0.3305846 0.0011755
165 Cap-dependent Translation Initiation 111 0.00e+00 0.3255501 0.0000003
396 Eukaryotic Translation Initiation 111 0.00e+00 0.3255501 0.0000003
828 Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC) 108 0.00e+00 0.3241781 0.0000004
830 Nonsense-Mediated Decay (NMD) 108 0.00e+00 0.3241781 0.0000004
703 Major pathway of rRNA processing in the nucleolus and cytosol 171 0.00e+00 0.3154260 0.0000000
1187 SRP-dependent cotranslational protein targeting to membrane 105 1.00e-07 0.2992844 0.0000066
1540 rRNA processing in the nucleus and cytosol 180 0.00e+00 0.2985439 0.0000000
191 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")

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
1384 TRAF6 mediated IRF7 activation 26 0.0000355 -0.4684102 0.0012593
212 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000006 -0.4614606 0.0000324
399 Expression and translocation of olfactory receptors 358 0.0000000 -0.4178164 0.0000000
853 Olfactory Signaling Pathway 365 0.0000000 -0.4085329 0.0000000
1215 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000768 -0.3567897 0.0024224
197 Chemokine receptors bind chemokines 56 0.0001806 -0.2892053 0.0053680
1214 Sensory Perception 569 0.0000000 -0.2715433 0.0000000
659 Keratinization 213 0.0000266 -0.1668486 0.0010826
460 G alpha (i) signalling events 305 0.0000624 -0.1331282 0.0020325
584 Immune System 1870 0.0001565 -0.0524093 0.0047443

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/RtmpJkCVLl/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/RtmpJkCVLl/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/RtmpJkCVLl/rmarkdown-str59af858f7bed0.html
## 
## Output created: /tmp/RtmpJkCVLl/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
212 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000000 -0.4901753 -0.4676449 0.0000001 0.0000004 0.6774685 0.0159314 0.0000000
1365 TICAM1-dependent activation of IRF3/IRF7 12 0.0012354 0.4502023 0.4919727 0.0069205 0.0031649 0.6668728 0.0295361 0.0075492
380 Eicosanoids 12 0.0002957 0.6313349 -0.1018713 0.0001522 0.5411614 0.6395010 0.5184551 0.0021354
1163 SARS-CoV-1 modulates host translation machinery 34 0.0000001 0.5215509 0.3046011 0.0000001 0.0021091 0.6039844 0.1534067 0.0000014
1215 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
224 Competing endogenous RNAs (ceRNAs) regulate PTEN translation 16 0.0001272 -0.4822519 0.2740224 0.0008369 0.0577137 0.5546667 0.5347667 0.0009980
516 Glucuronidation 25 0.0000806 -0.3782663 -0.3971243 0.0010586 0.0005866 0.5484461 0.0133346 0.0006701
1172 SCF(Skp2)-mediated degradation of p27/p21 58 0.0000000 0.4123054 0.3588126 0.0000001 0.0000023 0.5465731 0.0378251 0.0000000
492 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.0010509
903 Peptide chain elongation 84 0.0000000 0.4660659 0.2690810 0.0000000 0.0000200 0.5381654 0.1392894 0.0000000
882 PINK1-PRKN Mediated Mitophagy 21 0.0004631 0.2872506 0.4502581 0.0226638 0.0003534 0.5340836 0.1152637 0.0031822
397 Eukaryotic Translation Termination 87 0.0000000 0.4642281 0.2632838 0.0000000 0.0000217 0.5336910 0.1420890 0.0000000
1481 Ubiquitin Mediated Degradation of Phosphorylated Cdc25A 50 0.0000000 0.4249295 0.3125505 0.0000002 0.0001310 0.5274968 0.0794640 0.0000004
1533 p53-Independent DNA Damage Response 50 0.0000000 0.4249295 0.3125505 0.0000002 0.0001310 0.5274968 0.0794640 0.0000004
1534 p53-Independent G1/S DNA damage checkpoint 50 0.0000000 0.4249295 0.3125505 0.0000002 0.0001310 0.5274968 0.0794640 0.0000004
1216 Sensory perception of taste 47 0.0000001 -0.3649348 -0.3767662 0.0000149 0.0000078 0.5245285 0.0083661 0.0000013
106 Autodegradation of Cdh1 by Cdh1:APC/C 62 0.0000000 0.3716103 0.3700687 0.0000004 0.0000005 0.5244474 0.0010901 0.0000000
565 Hh mutants are degraded by ERAD 54 0.0000000 0.4078909 0.3275877 0.0000002 0.0000310 0.5231526 0.0567829 0.0000002
395 Eukaryotic Translation Elongation 88 0.0000000 0.4580045 0.2508669 0.0000000 0.0000471 0.5222090 0.1464684 0.0000000
107 Autodegradation of the E3 ubiquitin ligase COP1 50 0.0000000 0.4328964 0.2891483 0.0000001 0.0004030 0.5205824 0.1016453 0.0000005
1404 The role of GTSE1 in G2/M progression after G2 checkpoint 58 0.0000000 0.3997866 0.3330399 0.0000001 0.0000114 0.5203315 0.0471970 0.0000001
1482 Ubiquitin-dependent degradation of Cyclin D 50 0.0000000 0.4249949 0.2999003 0.0000002 0.0002428 0.5201546 0.0884553 0.0000005
1496 Viral mRNA Translation 84 0.0000000 0.4448583 0.2671498 0.0000000 0.0000229 0.5189103 0.1256589 0.0000000
1208 Selenocysteine synthesis 87 0.0000000 0.4484943 0.2576795 0.0000000 0.0000323 0.5172484 0.1349264 0.0000000
1502 Vpu mediated degradation of CD4 50 0.0000001 0.4050884 0.3180831 0.0000007 0.0000993 0.5150471 0.0615221 0.0000008
1323 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
807 Negative regulation of NOTCH4 signaling 52 0.0000000 0.3942742 0.3241131 0.0000009 0.0000524 0.5103934 0.0496114 0.0000006
1173 SCF-beta-TrCP mediated degradation of Emi1 53 0.0000000 0.3917813 0.3249161 0.0000008 0.0000425 0.5089823 0.0472809 0.0000005

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] 13
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

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

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.2.1                                    
##  [2] GGally_2.1.2                                       
##  [3] ggplot2_3.4.3                                      
##  [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] ebGSEA_0.1.0                                       
## [12] globaltest_5.54.0                                  
## [13] survival_3.5-7                                     
## [14] tictoc_1.2                                         
## [15] RIdeogram_0.2.2                                    
## [16] eulerr_7.0.0                                       
## [17] kableExtra_1.3.4                                   
## [18] data.table_1.14.8                                  
## [19] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [20] minfi_1.46.0                                       
## [21] bumphunter_1.42.0                                  
## [22] locfit_1.5-9.8                                     
## [23] iterators_1.0.14                                   
## [24] foreach_1.5.2                                      
## [25] Biostrings_2.68.1                                  
## [26] XVector_0.40.0                                     
## [27] SummarizedExperiment_1.30.2                        
## [28] Biobase_2.60.0                                     
## [29] MatrixGenerics_1.12.3                              
## [30] matrixStats_1.0.0                                  
## [31] GenomicRanges_1.52.0                               
## [32] GenomeInfoDb_1.36.3                                
## [33] IRanges_2.34.1                                     
## [34] S4Vectors_0.38.1                                   
## [35] BiocGenerics_0.46.0                                
## [36] mitch_1.12.0                                       
## [37] limma_3.56.2                                       
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.2             later_1.3.1              
##   [3] BiocIO_1.10.0             bitops_1.0-7             
##   [5] filelock_1.0.2            preprocessCore_1.62.1    
##   [7] XML_3.99-0.14             lifecycle_1.0.3          
##   [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.24            jquerylib_0.1.4          
##  [17] yaml_2.3.7                httpuv_1.6.11            
##  [19] grImport2_0.2-0           doRNG_1.8.6              
##  [21] askpass_1.2.0             DBI_1.1.3                
##  [23] RColorBrewer_1.1-3        abind_1.4-5              
##  [25] zlibbioc_1.46.0           rvest_1.0.3              
##  [27] quadprog_1.5-8            purrr_1.0.2              
##  [29] RCurl_1.98-1.12           rappdirs_0.3.3           
##  [31] GenomeInfoDbData_1.2.10   genefilter_1.82.1        
##  [33] annotate_1.78.0           svglite_2.1.1            
##  [35] DelayedMatrixStats_1.22.6 codetools_0.2-19         
##  [37] DelayedArray_0.26.7       xml2_1.3.5               
##  [39] tidyselect_1.2.0          farver_2.1.1             
##  [41] beanplot_1.3.1            base64enc_0.1-3          
##  [43] BiocFileCache_2.8.0       illuminaio_0.42.0        
##  [45] webshot_0.5.5             GenomicAlignments_1.36.0 
##  [47] jsonlite_1.8.7            multtest_2.56.0          
##  [49] ellipsis_0.3.2            systemfonts_1.0.4        
##  [51] tools_4.3.2               progress_1.2.2           
##  [53] Rcpp_1.0.11               glue_1.6.2               
##  [55] gridExtra_2.3             xfun_0.40                
##  [57] HDF5Array_1.28.1          withr_2.5.0              
##  [59] fastmap_1.1.1             rhdf5filters_1.12.1      
##  [61] fansi_1.0.4               openssl_2.1.0            
##  [63] caTools_1.18.2            digest_0.6.33            
##  [65] R6_2.5.1                  mime_0.12                
##  [67] colorspace_2.1-0          rsvg_2.4.0               
##  [69] jpeg_0.1-10               biomaRt_2.56.1           
##  [71] RSQLite_2.3.1             tidyr_1.3.0              
##  [73] utf8_1.2.3                generics_0.1.3           
##  [75] rtracklayer_1.60.1        prettyunits_1.1.1        
##  [77] httr_1.4.7                htmlwidgets_1.6.2        
##  [79] S4Arrays_1.0.6            pkgconfig_2.0.3          
##  [81] gtable_0.3.4              blob_1.2.4               
##  [83] siggenes_1.74.0           htmltools_0.5.6          
##  [85] scales_1.2.1              png_0.1-8                
##  [87] knitr_1.44                rstudioapi_0.15.0        
##  [89] tzdb_0.4.0                rjson_0.2.21             
##  [91] nlme_3.1-163              curl_5.0.2               
##  [93] org.Hs.eg.db_3.17.0       cachem_1.0.8             
##  [95] rhdf5_2.44.0              stringr_1.5.0            
##  [97] KernSmooth_2.23-22        AnnotationDbi_1.62.2     
##  [99] restfulr_0.0.15           GEOquery_2.68.0          
## [101] pillar_1.9.0              grid_4.3.2               
## [103] reshape_0.8.9             vctrs_0.6.3              
## [105] promises_1.2.1            dbplyr_2.3.3             
## [107] xtable_1.8-4              evaluate_0.21            
## [109] readr_2.1.4               GenomicFeatures_1.52.2   
## [111] cli_3.6.1                 compiler_4.3.2           
## [113] Rsamtools_2.16.0          rlang_1.1.1              
## [115] crayon_1.5.2              rngtools_1.5.2           
## [117] kpmt_0.1.0                labeling_0.4.3           
## [119] nor1mix_1.3-0             mclust_6.0.0             
## [121] plyr_1.8.8                stringi_1.7.12           
## [123] viridisLite_0.4.2         BiocParallel_1.34.2      
## [125] munsell_0.5.0             Matrix_1.6-1.1           
## [127] hms_1.1.3                 sparseMatrixStats_1.12.2 
## [129] bit64_4.0.5               Rhdf5lib_1.22.1          
## [131] KEGGREST_1.40.0           shiny_1.7.5              
## [133] highr_0.10                memoise_2.0.1            
## [135] bslib_0.5.1               bit_4.0.5