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("HGNChelper")
  library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
  source("meth_functions.R")
  library("data.table")
  library("kableExtra")
  library("eulerr")
  library("GenomicRanges")
#  data("dualmap850kEID")
  library("RhpcBLASctl")
})

source("meth_functions.R")

RhpcBLASctl::blas_set_num_threads(1)

Annotation

anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","UCSC_RefGene_Group","Islands_Name","Relation_to_Island")])
head(myann)
##            UCSC_RefGene_Name UCSC_RefGene_Group              Islands_Name
## cg18478105            YTHDF1             TSS200   chr20:61846843-61848103
## cg09835024            EIF2S3            TSS1500    chrX:24072558-24073135
## cg14361672              PKN3            TSS1500  chr9:131464843-131465830
## cg01763666            CCDC57               Body                          
## cg12950382         INF2;INF2          Body;Body                          
## cg02115394       CDC16;CDC16      TSS200;TSS200 chr13:115000148-115000874
##            Relation_to_Island
## cg18478105             Island
## cg09835024             Island
## cg14361672            N_Shore
## cg01763666            OpenSea
## cg12950382            OpenSea
## cg02115394             Island
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
gp2 <- lapply(gp2,unique)
gt <- stack(gp2)
colnames(gt) <- c("gene","probe")
gt$probe <- as.character(gt$probe)
dim(gt)
## [1] 684970      2
head(gt)
##     gene      probe
## 1 YTHDF1 cg18478105
## 2 EIF2S3 cg09835024
## 3   PKN3 cg14361672
## 4 CCDC57 cg01763666
## 5   INF2 cg12950382
## 6  CDC16 cg02115394
#update gene symbols
new.hgnc.table <- readRDS("new.hgnc.table.rds")
fix <- checkGeneSymbols(gt$gene,map=new.hgnc.table)
## Warning in checkGeneSymbols(gt$gene, map = new.hgnc.table): Human gene symbols
## should be all upper-case except for the 'orf' in open reading frames. The case
## of some letters was corrected.
## Warning in checkGeneSymbols(gt$gene, map = new.hgnc.table): x contains
## non-approved gene symbols
fix2 <- fix[which(fix$x != fix$Suggested.Symbol),]
length(unique(fix2$x))
## [1] 3254
gt$gene <- fix$Suggested.Symbol

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")

Mitch import

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

blm <- mitch_import(bl_lim,DEtype="limma",geneTable=gt,geneIDcol="Name")
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 802647
## Note: no. genes in output = 22285
## Warning in mitch_import(bl_lim, DEtype = "limma", geneTable = gt, geneIDcol = "Name"): Warning: less than half of the input genes are also in the
##         output
gum <- mitch_import(gu_lim,DEtype="limma",geneTable=gt,geneIDcol="Name")
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 790658
## Note: no. genes in output = 22269
## Warning in mitch_import(gu_lim, DEtype = "limma", geneTable = gt, geneIDcol = "Name"): Warning: less than half of the input genes are also in the
##         output
bum <- mitch_import(bu_lim,DEtype="limma",geneTable=gt,geneIDcol="Name")
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 801260
## Note: no. genes in output = 22286
## Warning in mitch_import(bu_lim, DEtype = "limma", geneTable = gt, geneIDcol = "Name"): Warning: less than half of the input genes are also in the
##         output

Probe level t-stats

dim(bl_lim)
## [1] 802647     14
par(mfrow=c(2,1))

hist(bl_lim$t,xlab="probe t-stat",main="blood at assessment", breaks=30, xlim=c(-4,4))
grid() ; abline(v=0,lty=1,lwd=3)
hist(blm$x,xlab="gene mean t-stat",main="blood at assessment", breaks=30, xlim=c(-2,2))
grid() ; abline(v=0,lty=1,lwd=3)

hist(gu_lim$t,xlab="probe t-stat",main="neonatal Guthrie card", breaks=30, xlim=c(-4,4))
grid() ; abline(v=0,lty=1,lwd=3)
hist(gum$x,xlab="gene mean t-stat",main="neonatal Guthrie card", breaks=30, xlim=c(-2,2))
grid() ; abline(v=0,lty=1,lwd=3)

hist(bu_lim$t,xlab="probe t-stat",main="buccal swab", breaks=30, xlim=c(-4,4))
grid() ; abline(v=0,lty=1,lwd=3)
hist(bum$x,xlab="gene mean t-stat",main="buccal swab", breaks=30, xlim=c(-2,2))
grid() ; abline(v=0,lty=1,lwd=3)

pdf("hist_bl_ados.pdf")
par(mfrow=c(2,1))
hist(bl_lim$t,xlab="probe t-stat",main="blood at assessment", breaks=30, xlim=c(-4,4))
grid() ; abline(v=0,lty=1,lwd=3)
hist(blm$x,xlab="gene mean t-stat",main="blood at assessment", breaks=30, xlim=c(-2,2))
grid() ; abline(v=0,lty=1,lwd=3)
dev.off()
## png 
##   2
pdf("hist_gu_ados.pdf")
hist(gu_lim$t,xlab="probe t-stat",main="neonatal Guthrie card", breaks=30, xlim=c(-4,4))
grid() ; abline(v=0,lty=1,lwd=3)
hist(gum$x,xlab="gene mean t-stat",main="neonatal Guthrie card", breaks=30, xlim=c(-2,2))
grid() ; abline(v=0,lty=1,lwd=3)
dev.off()
## png 
##   2
par(mfrow=c(1,1))

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, h=0, lty=2, lwd=2, col="blue")
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

Gene set level enrichment

par(mar = c(5.1, 26.1, 4.1, 2.1))

msigdb <- gmt_import("msigdb.v2023.2.Hs.symbols.gmt")

gobp <- msigdb[grep("^GOBP",names(msigdb))]
gomf <- msigdb[grep("^GOMF",names(msigdb))]
gocc <- msigdb[grep("^GOCC",names(msigdb))]
go <- c(gobp,gomf,gocc)
names(go) <- gsub("_"," ",names(go))

reactome <- gmt_import("ReactomePathways_2024-09-18.gmt")

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

cores = 8

## BLOOD
# go
bl_mres1 <- mitch_calc(blm,go,priority="effect",cores=cores,minsetsize=5)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
bl_mres1tbl <- subset(bl_mres1$enrichment_result,p.adjustANOVA<0.01)
top <- head(bl_mres1tbl,30)

head(top[order(-abs(top$s.dist)),],30) %>%
  kbl(caption="BLOOD: Top effect size GOs found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)
BLOOD: Top effect size GOs found with mitch after 1% FDR filtering
set setSize pANOVA s.dist p.adjustANOVA
7737 GOMF BITTER TASTE RECEPTOR ACTIVITY 24 0.00e+00 -0.8235629 0.0000000
9098 GOMF TASTE RECEPTOR ACTIVITY 30 0.00e+00 -0.6360009 0.0000004
9058 GOMF STRUCTURAL CONSTITUENT OF CHROMATIN 80 0.00e+00 0.5544292 0.0000000
9513 GOCC CYTOSOLIC SMALL RIBOSOMAL SUBUNIT 37 0.00e+00 0.5192766 0.0000088
10005 GOCC PROTEASOME CORE COMPLEX 20 9.28e-05 0.5047743 0.0084929
9510 GOCC CYTOSOLIC LARGE RIBOSOMAL SUBUNIT 54 0.00e+00 0.4579300 0.0000013
9887 GOCC NUCLEOSOME 116 0.00e+00 0.4562063 0.0000000
1452 GOBP EXCRETION 29 6.37e-05 0.4289211 0.0061609
8600 GOMF ODORANT BINDING 104 0.00e+00 -0.4274379 0.0000000
8835 GOMF PROTEIN CARRIER CHAPERONE 33 2.43e-05 0.4245420 0.0027021
10111 GOCC SMALL RIBOSOMAL SUBUNIT 71 0.00e+00 0.4030221 0.0000010
9512 GOCC CYTOSOLIC RIBOSOME 108 0.00e+00 0.3914337 0.0000000
9067 GOMF STRUCTURAL CONSTITUENT OF RIBOSOME 156 0.00e+00 0.3866463 0.0000000
10052 GOCC RIBOSOMAL SUBUNIT 178 0.00e+00 0.3820720 0.0000000
6859 GOBP SENSORY PERCEPTION OF BITTER TASTE 44 1.41e-05 -0.3782674 0.0016998
5067 GOBP PROTON MOTIVE FORCE DRIVEN ATP SYNTHESIS 65 2.00e-07 0.3764149 0.0000264
8521 GOMF NAD P H DEHYDROGENASE QUINONE ACTIVITY 37 7.42e-05 0.3763764 0.0070369
9021 GOMF SNRNA BINDING 48 7.70e-06 0.3731843 0.0009540
9730 GOCC LARGE RIBOSOMAL SUBUNIT 111 0.00e+00 0.3678348 0.0000000
1008 GOBP CYTOPLASMIC TRANSLATION 141 0.00e+00 0.3444921 0.0000000
8662 GOMF OXIDOREDUCTION DRIVEN ACTIVE TRANSMEMBRANE TRANSPORTER ACTIVITY 53 1.55e-05 0.3431073 0.0018432
6775 GOBP RIBOSOMAL SMALL SUBUNIT BIOGENESIS 101 0.00e+00 0.3351064 0.0000013
2347 GOBP MATURATION OF SSU RRNA 51 5.71e-05 0.3256768 0.0056317
10044 GOCC RESPIRASOME 85 4.00e-07 0.3176767 0.0000656
3629 GOBP NUCLEOSOME ORGANIZATION 119 0.00e+00 0.3125481 0.0000009
10004 GOCC PROTEASOME COMPLEX 59 4.77e-05 0.3060288 0.0049390
10053 GOCC RIBOSOME 225 0.00e+00 0.3052046 0.0000000
278 GOBP ATP SYNTHESIS COUPLED ELECTRON TRANSPORT 85 2.30e-06 0.2961706 0.0003245
8971 GOMF RRNA BINDING 62 5.99e-05 0.2946323 0.0058480
4916 GOBP PROTEIN DNA COMPLEX ASSEMBLY 203 0.00e+00 0.2921723 0.0000000 
# Barplot
up <- head(subset(bl_mres1$enrichment_result,p.adjustANOVA<0.01 & s.dist>0),15)
vup <- up$s.dist
names(vup) <- up$set
dn <- head(subset(bl_mres1$enrichment_result,p.adjustANOVA<0.01 & s.dist<0),15)
vdn <- dn$s.dist
names(vdn) <- dn$set
barplot(sort(c(vup,vdn)),horiz=TRUE,las=1,cex.names=0.8,xlab="ES",main="Blood at assessment")
mtext("GO")

pdf("ados_barplot_bl_go.pdf")
par(mar = c(5.1, 29.1, 4.1, 2.1))
barplot(sort(c(vup,vdn)),horiz=TRUE,las=1,cex.names=0.9,xlab="ES",main="Blood at assessment")
mtext("GO")
dev.off()
## png 
##   2
# reactome
bl_mres2 <- mitch_calc(blm,reactome,priority="effect",cores=cores,minsetsize=5)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
bl_mres2tbl <- subset(bl_mres2$enrichment_result,p.adjustANOVA<0.01)
top <- head(bl_mres2tbl,30)

head(top[order(-abs(top$s.dist)),],30) %>%
  kbl(caption="BLOOD: Top effect size Reactome pathways found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)
BLOOD: Top effect size Reactome pathways found with mitch after 1% FDR filtering
set setSize pANOVA s.dist p.adjustANOVA
501 Eicosanoids 12 0.0001260 0.6391146 0.0011508
1338 RNA Polymerase I Promoter Opening 26 0.0000000 0.6321954 0.0000008
1201 Pexophagy 10 0.0009548 0.6032323 0.0071598
1551 SIRT1 negatively regulates rRNA expression 30 0.0000000 0.5845503 0.0000010
1530 SARS-CoV-1 modulates host translation machinery 34 0.0000000 0.5511160 0.0000009
1185 Packaging Of Telomere Ends 29 0.0000003 0.5496077 0.0000067
787 Inhibition of DNA recombination at telomere 43 0.0000000 0.5350176 0.0000001
286 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000000 -0.5064904 0.0000013
359 DNA methylation 28 0.0000038 0.5046647 0.0000504
527 Eukaryotic Translation Termination 87 0.0000000 0.4901141 0.0000000
1194 Peptide chain elongation 84 0.0000000 0.4886089 0.0000000
273 Chromatin modifications during the maternal to zygotic transition (MZT) 32 0.0000022 0.4830458 0.0000323
525 Eukaryotic Translation Elongation 88 0.0000000 0.4821736 0.0000000
1600 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000001 -0.4809650 0.0000027
138 Assembly of the ORC complex at the origin of replication 32 0.0000030 0.4768332 0.0000411
416 Deposition of new CENPA-containing nucleosomes at the centromere 46 0.0000000 0.4750810 0.0000008
1116 Nucleosome assembly 46 0.0000000 0.4750810 0.0000008
1495 Replacement of protamines by nucleosomes in the male pronucleus 22 0.0001237 0.4726779 0.0011350
1590 Selenocysteine synthesis 87 0.0000000 0.4706539 0.0000000
1961 Viral mRNA Translation 84 0.0000000 0.4696297 0.0000000
49 Activated PKN1 stimulates transcription of AR (androgen receptor) regulated genes KLK2 and KLK3 29 0.0000154 0.4636193 0.0001736
1468 Regulation of endogenous retroelements by the Human Silencing Hub (HUSH) complex 32 0.0000078 0.4564651 0.0000973
1102 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 89 0.0000000 0.4558332 0.0000000
308 Condensation of Prophase Chromosomes 36 0.0000026 0.4525010 0.0000364
1509 Response of EIF2AK4 (GCN2) to amino acid deficiency 95 0.0000000 0.4512108 0.0000000
292 Cleavage of the damaged pyrimidine 38 0.0000015 0.4506817 0.0000233
418 Depyrimidination 38 0.0000015 0.4506817 0.0000233
1397 Recognition and association of DNA glycosylase with site containing an affected pyrimidine 38 0.0000015 0.4506817 0.0000233
401 Defective pyroptosis 34 0.0000068 0.4457156 0.0000855
152 B-WICH complex positively regulates rRNA expression 54 0.0000000 0.4444203 0.0000006 
# Barplot
up <- head(subset(bl_mres2$enrichment_result,p.adjustANOVA<0.01 & s.dist>0),15)
vup <- up$s.dist
names(vup) <- up$set
dn <- head(subset(bl_mres2$enrichment_result,p.adjustANOVA<0.01 & s.dist<0),15)
vdn <- dn$s.dist
names(vdn) <- dn$set
barplot(sort(c(vup,vdn)),horiz=TRUE,las=1,cex.names=0.9,xlab="ES",main="Blood at assessment")
mtext("Reactome")

pdf("ados_barplot_bl_reactome.pdf")
par(mar = c(5.1, 29.1, 4.1, 2.1))
barplot(sort(c(vup,vdn)),horiz=TRUE,las=1,cex.names=0.9,xlab="ES",main="Blood at assessment")
mtext("Reactome")
dev.off()
## png 
##   2
## GUTHRIE
# GO
gu_mres1 <- mitch_calc(gum,go,priority="effect",cores=cores,minsetsize=5)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
gu_mres1tbl <- subset(gu_mres1$enrichment_result,p.adjustANOVA<0.01)
top <- head(gu_mres1tbl,30)

head(top[order(-abs(top$s.dist)),],30) %>%
  kbl(caption="GUTHRIE: Top effect size GOs found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)
GUTHRIE: Top effect size GOs found with mitch after 1% FDR filtering
set setSize pANOVA s.dist p.adjustANOVA
8455 GOMF MITOCHONDRIAL RIBOSOME BINDING 9 0.0003252 0.6918139 0.0085677
7737 GOMF BITTER TASTE RECEPTOR ACTIVITY 24 0.0000000 -0.6653780 0.0000014
6809 GOBP ROSTROCAUDAL NEURAL TUBE PATTERNING 10 0.0003625 0.6511434 0.0093353
5079 GOBP PTERIDINE CONTAINING COMPOUND BIOSYNTHETIC PROCESS 12 0.0003691 0.5936484 0.0094339
557 GOBP CELLULAR GLUCURONIDATION 21 0.0000100 -0.5565829 0.0004970
7253 GOBP TRANSLATIONAL TERMINATION 15 0.0002359 0.5482879 0.0066048
4537 GOBP POSITIVE REGULATION OF PEPTIDYL SERINE PHOSPHORYLATION OF STAT PROTEIN 16 0.0002912 -0.5230699 0.0078760
6859 GOBP SENSORY PERCEPTION OF BITTER TASTE 44 0.0000000 -0.5023111 0.0000007
9098 GOMF TASTE RECEPTOR ACTIVITY 30 0.0000026 -0.4960025 0.0001396
8956 GOMF RNA POLYMERASE II CORE PROMOTER SEQUENCE SPECIFIC DNA BINDING 18 0.0003136 0.4905647 0.0083277
8600 GOMF ODORANT BINDING 103 0.0000000 -0.4855928 0.0000000
4321 GOBP POSITIVE REGULATION OF INTERFERON BETA PRODUCTION 37 0.0000004 0.4836005 0.0000236
5025 GOBP PROTEIN PEPTIDYL PROLYL ISOMERIZATION 19 0.0003651 0.4722413 0.0093773
3942 GOBP POLYPRENOL METABOLIC PROCESS 23 0.0001179 0.4637238 0.0038489
9058 GOMF STRUCTURAL CONSTITUENT OF CHROMATIN 80 0.0000000 0.4620634 0.0000000
7852 GOMF CIS TRANS ISOMERASE ACTIVITY 39 0.0000009 0.4549454 0.0000522
1080 GOBP DETECTION OF CHEMICAL STIMULUS INVOLVED IN SENSORY PERCEPTION OF TASTE 46 0.0000002 -0.4396757 0.0000170
176 GOBP ANAPHASE PROMOTING COMPLEX DEPENDENT CATABOLIC PROCESS 24 0.0002601 0.4305949 0.0071472
8602 GOMF OLFACTORY RECEPTOR ACTIVITY 354 0.0000000 -0.4190990 0.0000000
7468 GOBP VIRAL BUDDING 25 0.0003568 0.4124402 0.0092109
7952 GOMF DISULFIDE OXIDOREDUCTASE ACTIVITY 34 0.0000387 0.4077382 0.0015848
8521 GOMF NAD P H DEHYDROGENASE QUINONE ACTIVITY 37 0.0000192 0.4059598 0.0008939
7387 GOBP URONIC ACID METABOLIC PROCESS 26 0.0003909 -0.4017236 0.0097714
1961 GOBP INTERFERON BETA PRODUCTION 56 0.0000002 0.3990215 0.0000166
8662 GOMF OXIDOREDUCTION DRIVEN ACTIVE TRANSMEMBRANE TRANSPORTER ACTIVITY 53 0.0000006 0.3961958 0.0000373
6864 GOBP SENSORY PERCEPTION OF SMELL 380 0.0000000 -0.3861774 0.0000000
6860 GOBP SENSORY PERCEPTION OF CHEMICAL STIMULUS 455 0.0000000 -0.3796892 0.0000000
8835 GOMF PROTEIN CARRIER CHAPERONE 33 0.0001714 0.3779157 0.0050436
10044 GOCC RESPIRASOME 85 0.0000000 0.3707781 0.0000003
1093 GOBP DETECTION OF STIMULUS INVOLVED IN SENSORY PERCEPTION 469 0.0000000 -0.3696058 0.0000000 
# Barplot
up <- head(subset(gu_mres1$enrichment_result,p.adjustANOVA<0.01 & s.dist>0),15)
vup <- up$s.dist
names(vup) <- up$set
dn <- head(subset(gu_mres1$enrichment_result,p.adjustANOVA<0.01 & s.dist<0),15)
vdn <- dn$s.dist
names(vdn) <- dn$set
barplot(sort(c(vup,vdn)),horiz=TRUE,las=1,cex.names=0.8,xlab="ES",main="Neonatal Guthrie card")
mtext("GO")

pdf("ados_barplot_gu_go.pdf")
par(mar = c(5.1, 29.1, 4.1, 2.1))
barplot(sort(c(vup,vdn)),horiz=TRUE,las=1,cex.names=0.7,xlab="ES",main="Neonatal Guthrie card")
mtext("GO")
dev.off()
## png 
##   2
# reactome
gu_mres2 <- mitch_calc(gum,reactome,priority="effect",cores=cores,minsetsize=5)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
gu_mres2tbl <- subset(gu_mres2$enrichment_result,p.adjustANOVA<0.01)
top <- head(gu_mres2tbl,30)

head(top[order(-abs(top$s.dist)),],30) %>%
  kbl(caption="GUTHRIE: Top effect size Reactome pathways found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)
GUTHRIE: Top effect size Reactome pathways found with mitch after 1% FDR filtering
set setSize pANOVA s.dist p.adjustANOVA
305 Complex III assembly 24 0.0000134 0.5132689 0.0001885
138 Assembly of the ORC complex at the origin of replication 32 0.0000022 0.4836941 0.0000445
16 APC/C:Cdc20 mediated degradation of Cyclin B 24 0.0000423 0.4827752 0.0004654
286 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000002 -0.4815622 0.0000056
1338 RNA Polymerase I Promoter Opening 26 0.0000213 0.4814739 0.0002687
401 Defective pyroptosis 34 0.0000012 0.4811347 0.0000275
1468 Regulation of endogenous retroelements by the Human Silencing Hub (HUSH) complex 32 0.0000031 0.4760633 0.0000582
1185 Packaging Of Telomere Ends 29 0.0000125 0.4686554 0.0001793
1178 PRC2 methylates histones and DNA 35 0.0000021 0.4633226 0.0000440
134 Assembly Of The HIV Virion 16 0.0014511 0.4597526 0.0082389
49 Activated PKN1 stimulates transcription of AR (androgen receptor) regulated genes KLK2 and KLK3 29 0.0000208 0.4565694 0.0002631
1220 Phosphorylation of the APC/C 20 0.0005181 0.4482988 0.0034890
1886 Transcriptional regulation of granulopoiesis 53 0.0000000 0.4469233 0.0000007
416 Deposition of new CENPA-containing nucleosomes at the centromere 46 0.0000002 0.4466201 0.0000047
1116 Nucleosome assembly 46 0.0000002 0.4466201 0.0000047
158 Base-Excision Repair, AP Site Formation 40 0.0000011 0.4453349 0.0000259
292 Cleavage of the damaged pyrimidine 38 0.0000021 0.4443108 0.0000440
418 Depyrimidination 38 0.0000021 0.4443108 0.0000440
1397 Recognition and association of DNA glycosylase with site containing an affected pyrimidine 38 0.0000021 0.4443108 0.0000440
1551 SIRT1 negatively regulates rRNA expression 30 0.0000275 0.4422756 0.0003255
34 Aberrant regulation of mitotic exit in cancer due to RB1 defects 20 0.0006341 0.4412468 0.0041594
14 APC-Cdc20 mediated degradation of Nek2A 26 0.0001052 0.4393258 0.0009876
493 ERCC6 (CSB) and EHMT2 (G9a) positively regulate rRNA expression 38 0.0000029 0.4384690 0.0000558
424 Digestion 18 0.0013510 -0.4362850 0.0077797
359 DNA methylation 28 0.0000692 0.4343507 0.0007198
1164 PINK1-PRKN Mediated Mitophagy 31 0.0000287 0.4340405 0.0003365
321 Conversion from APC/C:Cdc20 to APC/C:Cdh1 in late anaphase 20 0.0008135 0.4324104 0.0051209
291 Cleavage of the damaged purine 33 0.0000176 0.4317132 0.0002336
417 Depurination 33 0.0000176 0.4317132 0.0002336
1396 Recognition and association of DNA glycosylase with site containing an affected purine 33 0.0000176 0.4317132 0.0002336 
# Barplot
up <- head(subset(gu_mres2$enrichment_result,p.adjustANOVA<0.01 & s.dist>0),15)
vup <- up$s.dist
names(vup) <- up$set
dn <- head(subset(gu_mres2$enrichment_result,p.adjustANOVA<0.01 & s.dist<0),15)
vdn <- dn$s.dist
names(vdn) <- dn$set
barplot(sort(c(vup,vdn)),horiz=TRUE,las=1,cex.names=0.9,xlab="ES",main="Neonatal Guthrie card")
mtext("Reactome")

pdf("ados_barplot_gu_reactome.pdf")
par(mar = c(5.1, 29.1, 4.1, 2.1))
barplot(sort(c(vup,vdn)),horiz=TRUE,las=1,cex.names=0.75,xlab="ES",main="Neonatal Guthrie card")
mtext("Reactome")
dev.off()
## png 
##   2
par(mar = c(5.1, 4.1, 4.1, 2.1))

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
res1 <- mitch_calc(m, go, priority="effect",cores=cores)
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
mitch_report(res1,"blgu_mitch_go.html",overwrite=TRUE)
## Dataset saved as " /tmp/RtmpGyw2UW/blgu_mitch_go.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: /asd_meth/mitch.knit.md
## /usr/bin/pandoc +RTS -K512m -RTS /asd_meth/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpGyw2UW/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 --embed-resources --standalone --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/RtmpGyw2UW/rmarkdown-str1ce332e6d2.html
## 
## Output created: /tmp/RtmpGyw2UW/mitch_report.html
## [1] TRUE
mitch_plots(res=res1,outfile="blgu_mitch_go_ados.pdf")
## png 
##   2
sig1 <- subset(res1$enrichment_result,p.adjustMANOVA<0.01)

head(sig1[order(-abs(sig1$s.dist)),],30) %>%
  kbl(caption="Top GOs found with mitch (effect size after 1% FDR filtering)") %>%
  kable_paper("hover", full_width = F)
Top GOs 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
5599 GOMF BITTER TASTE RECEPTOR ACTIVITY 24 0.0000000 -0.8243351 -0.6653780 0.0000000 0.0000000 1.0593659 0.1123996 0.0000000
6498 GOMF TASTE RECEPTOR ACTIVITY 30 0.0000000 -0.6367073 -0.4960025 0.0000000 0.0000026 0.8071026 0.0994933 0.0000000
6561 GOMF TYPE I INTERFERON RECEPTOR BINDING 13 0.0000911 -0.5293093 -0.5440676 0.0009505 0.0006814 0.7590638 0.0104357 0.0022587
1821 GOBP MITOCHONDRIAL ELECTRON TRANSPORT UBIQUINOL TO CYTOCHROME C 12 0.0003427 0.5657995 0.4608138 0.0006885 0.0057059 0.7297112 0.0742361 0.0065738
406 GOBP CELLULAR GLUCURONIDATION 21 0.0000012 -0.4603731 -0.5565829 0.0002595 0.0000100 0.7223074 0.0680306 0.0000456
6470 GOMF STRUCTURAL CONSTITUENT OF CHROMATIN 80 0.0000000 0.5542960 0.4620634 0.0000000 0.0000000 0.7216278 0.0652183 0.0000000
3604 GOBP PROTEIN LOCALIZATION TO CENP A CONTAINING CHROMATIN 16 0.0000837 0.5487856 0.4076304 0.0001440 0.0047543 0.6836140 0.0998117 0.0020883
3301 GOBP POSITIVE REGULATION OF PEPTIDYL SERINE PHOSPHORYLATION OF STAT PROTEIN 16 0.0001008 -0.4329192 -0.5230699 0.0027139 0.0002912 0.6789854 0.0637462 0.0024747
5294 GOBP TRANSLATIONAL TERMINATION 15 0.0003030 0.3525898 0.5482879 0.0180562 0.0002359 0.6518736 0.1383795 0.0060339
6709 GOCC CHROMOSOME CENTROMERIC CORE DOMAIN 17 0.0001418 0.5052183 0.4021952 0.0003098 0.0040884 0.6457604 0.0728484 0.0032705
6149 GOMF ODORANT BINDING 103 0.0000000 -0.4225346 -0.4855928 0.0000000 0.0000000 0.6436893 0.0445889 0.0000000
5008 GOBP SENSORY PERCEPTION OF BITTER TASTE 44 0.0000000 -0.3788506 -0.5023111 0.0000137 0.0000000 0.6291615 0.0872997 0.0000000
6385 GOMF RIBOSOMAL SMALL SUBUNIT BINDING 16 0.0003917 0.3763987 0.4992417 0.0091371 0.0005445 0.6252345 0.0868631 0.0073782
2268 GOBP NEGATIVE REGULATION OF MEGAKARYOCYTE DIFFERENTIATION 18 0.0002681 0.5142091 0.3031974 0.0001584 0.0259383 0.5969419 0.1492078 0.0055055
6781 GOCC CYTOSOLIC SMALL RIBOSOMAL SUBUNIT 37 0.0000001 0.5191056 0.2739659 0.0000000 0.0039251 0.5869651 0.1733399 0.0000033
6310 GOMF PROTEIN CARRIER CHAPERONE 33 0.0000018 0.4243106 0.3779157 0.0000245 0.0001714 0.5682075 0.0328061 0.0000668
7031 GOCC NUCLEOSOME 116 0.0000000 0.4560035 0.3167104 0.0000000 0.0000000 0.5551979 0.0984951 0.0000000
5386 GOBP URONIC ACID METABOLIC PROCESS 26 0.0000482 -0.3827548 -0.4017236 0.0007284 0.0003909 0.5548722 0.0134130 0.0013170
7119 GOCC PROTEASOME CORE COMPLEX 20 0.0002973 0.5045980 0.2293856 0.0000933 0.0757365 0.5542895 0.1946045 0.0059643
6095 GOMF NAD P H DEHYDROGENASE QUINONE ACTIVITY 37 0.0000008 0.3761099 0.4059598 0.0000750 0.0000192 0.5534094 0.0211070 0.0000319
6778 GOCC CYTOSOLIC LARGE RIBOSOMAL SUBUNIT 54 0.0000000 0.4577263 0.3009887 0.0000000 0.0001299 0.5478207 0.1108302 0.0000001
788 GOBP DETECTION OF CHEMICAL STIMULUS INVOLVED IN SENSORY PERCEPTION OF TASTE 46 0.0000001 -0.3074664 -0.4396757 0.0003080 0.0000002 0.5365168 0.0934861 0.0000029
3135 GOBP POSITIVE REGULATION OF INTERFERON BETA PRODUCTION 37 0.0000012 0.2097537 0.4836005 0.0272406 0.0000004 0.5271300 0.1936389 0.0000456
3678 GOBP PROTON MOTIVE FORCE DRIVEN ATP SYNTHESIS 65 0.0000000 0.3761484 0.3676510 0.0000002 0.0000003 0.5259800 0.0060086 0.0000000
6193 GOMF OXIDOREDUCTION DRIVEN ACTIVE TRANSMEMBRANE TRANSPORTER ACTIVITY 53 0.0000000 0.3428194 0.3961958 0.0000157 0.0000006 0.5239239 0.0377429 0.0000008
6150 GOMF OLFACTORY RECEPTOR ACTIVITY 354 0.0000000 -0.2855009 -0.4190990 0.0000000 0.0000000 0.5071042 0.0944681 0.0000000
1815 GOBP MITOCHONDRIAL CYTOCHROME C OXIDASE ASSEMBLY 25 0.0004709 0.3891458 0.3060133 0.0007559 0.0080796 0.4950542 0.0587836 0.0085812
6992 GOCC NADH DEHYDROGENASE COMPLEX 40 0.0000054 0.3296010 0.3679540 0.0003089 0.0000563 0.4939909 0.0271197 0.0001893
725 GOBP CYTOCHROME COMPLEX ASSEMBLY 38 0.0000123 0.3231239 0.3671758 0.0005661 0.0000895 0.4891085 0.0311493 0.0003929
7146 GOCC RESPIRASOME 85 0.0000000 0.3173745 0.3707781 0.0000004 0.0000000 0.4880604 0.0377621 0.0000000 
res2 <- mitch_calc(m, reactome, priority="effect",cores=cores)
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
mitch_report(res2,"blgu_mitch_reactome.html",overwrite=TRUE)
## Dataset saved as " /tmp/RtmpGyw2UW/blgu_mitch_reactome.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: /asd_meth/mitch.knit.md
## /usr/bin/pandoc +RTS -K512m -RTS /asd_meth/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpGyw2UW/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 --embed-resources --standalone --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/RtmpGyw2UW/rmarkdown-str1ce8a4fa71.html
## 
## Output created: /tmp/RtmpGyw2UW/mitch_report.html
## [1] TRUE
mitch_plots(res=res2,outfile="blgu_mitch_reactome_ados.pdf")
## Warning in title(main = ll$set, xlab = X_AXIS, ylab = Y_AXIS): conversion
## failure on 'Regulation of TBK1, IKKε-mediated activation of IRF3, IRF7 upon
## TLR3 ligation' in 'mbcsToSbcs': for ε (U+03B5)
## Warning in title(...): conversion failure on 'Regulation of TBK1, IKKε-mediated
## activation of IRF3, IRF7 upon TLR3 ligation' in 'mbcsToSbcs': for ε (U+03B5)
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Regulation of TBK1, IKKε-mediated activation of IRF3,
## IRF7 upon TLR3 ligation' in 'mbcsToSbcs': for ε (U+03B5)
## png 
##   2
sig2 <- subset(res2$enrichment_result,p.adjustMANOVA<0.01)

head(sig2[order(-abs(sig2$s.dist)),],30) %>%
  kbl(caption="Top Reactomes found with mitch (effect size after 1% FDR filtering)") %>%
  kable_paper("hover", full_width = F)
Top Reactomes 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
1072 RNA Polymerase I Promoter Opening 26 0.0000000 0.6320969 0.4814739 0.0000000 0.0000213 0.7945839 0.1065066 0.0000000
954 Pexophagy 10 0.0010451 0.6031268 0.4249337 0.0009568 0.0199681 0.7377877 0.1260015 0.0048149
1247 SIRT1 negatively regulates rRNA expression 30 0.0000000 0.5844193 0.4422756 0.0000000 0.0000275 0.7329076 0.1005108 0.0000000
938 Packaging Of Telomere Ends 29 0.0000000 0.5494542 0.4686554 0.0000003 0.0000125 0.7221758 0.0571334 0.0000001
220 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000000 -0.5071617 -0.4815622 0.0000000 0.0000002 0.6993677 0.0181016 0.0000000
105 Assembly of the ORC complex at the origin of replication 32 0.0000000 0.4766310 0.4836941 0.0000030 0.0000022 0.6790708 0.0049944 0.0000001
281 DNA methylation 28 0.0000002 0.5044753 0.4343507 0.0000038 0.0000692 0.6656995 0.0495856 0.0000018
1176 Regulation of endogenous retroelements by the Human Silencing Hub (HUSH) complex 32 0.0000000 0.4562483 0.4760633 0.0000079 0.0000031 0.6593927 0.0140113 0.0000004
307 Defective pyroptosis 34 0.0000000 0.4455019 0.4811347 0.0000069 0.0000012 0.6557153 0.0251962 0.0000002
320 Deposition of new CENPA-containing nucleosomes at the centromere 46 0.0000000 0.4748772 0.4466201 0.0000000 0.0000002 0.6519033 0.0199808 0.0000000
884 Nucleosome assembly 46 0.0000000 0.4748772 0.4466201 0.0000000 0.0000002 0.6519033 0.0199808 0.0000000
33 Activated PKN1 stimulates transcription of AR (androgen receptor) regulated genes KLK2 and KLK3 29 0.0000002 0.4634039 0.4565694 0.0000156 0.0000208 0.6505373 0.0048327 0.0000022
392 Eicosanoids 12 0.0001985 0.6390274 -0.1175660 0.0001263 0.4807042 0.6497521 0.5349923 0.0010513
627 Inhibition of DNA recombination at telomere 43 0.0000000 0.5348544 0.3517042 0.0000000 0.0000657 0.6401290 0.1295068 0.0000000
1228 SARS-CoV-1 modulates host translation machinery 34 0.0000000 0.5509650 0.3111311 0.0000000 0.0016892 0.6327440 0.1695881 0.0000003
225 Cleavage of the damaged pyrimidine 38 0.0000000 0.4504663 0.4443108 0.0000015 0.0000021 0.6327179 0.0043526 0.0000001
322 Depyrimidination 38 0.0000000 0.4504663 0.4443108 0.0000015 0.0000021 0.6327179 0.0043526 0.0000001
1114 Recognition and association of DNA glycosylase with site containing an affected pyrimidine 38 0.0000000 0.4504663 0.4443108 0.0000015 0.0000021 0.6327179 0.0043526 0.0000001
1441 TICAM1-dependent activation of IRF3/IRF7 13 0.0019033 0.4465589 0.4342615 0.0053029 0.0067032 0.6228948 0.0086956 0.0079803
209 Chromatin modifications during the maternal to zygotic transition (MZT) 32 0.0000002 0.4828453 0.3835682 0.0000023 0.0001727 0.6166556 0.0701995 0.0000022
224 Cleavage of the damaged purine 33 0.0000002 0.4374560 0.4317132 0.0000136 0.0000176 0.6146089 0.0040608 0.0000018
321 Depurination 33 0.0000002 0.4374560 0.4317132 0.0000136 0.0000176 0.6146089 0.0040608 0.0000018
1113 Recognition and association of DNA glycosylase with site containing an affected purine 33 0.0000002 0.4374560 0.4317132 0.0000136 0.0000176 0.6146089 0.0040608 0.0000018
1282 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000000 -0.4816470 -0.3788695 0.0000001 0.0000269 0.6128017 0.0726746 0.0000001
238 Condensation of Prophase Chromosomes 36 0.0000001 0.4522844 0.4064616 0.0000026 0.0000243 0.6080890 0.0324016 0.0000007
121 Base-Excision Repair, AP Site Formation 40 0.0000000 0.4117459 0.4453349 0.0000066 0.0000011 0.6065129 0.0237510 0.0000002
235 Complex III assembly 24 0.0000154 0.3131678 0.5132689 0.0079086 0.0000134 0.6012645 0.1414929 0.0000998
384 ERCC6 (CSB) and EHMT2 (G9a) positively regulate rRNA expression 38 0.0000001 0.3932595 0.4384690 0.0000272 0.0000029 0.5889891 0.0319680 0.0000008
935 PRC2 methylates histones and DNA 35 0.0000003 0.3630514 0.4633226 0.0002011 0.0000021 0.5886206 0.0709024 0.0000023
1280 Senescence-Associated Secretory Phenotype (SASP) 74 0.0000000 0.4004682 0.4268231 0.0000000 0.0000000 0.5852800 0.0186357 0.0000000

Top table of genes

rnk <- res1$ranked_profile

tbl <- data.frame(table(gt$gene))

toptab <- cbind(m,rnk)

toptab$Row.names = NULL

toptab$med <- apply(data.frame(toptab$bl,toptab$gu),1,median)

toptab2 <- merge(toptab,tbl,by.x=0,by.y="Var1")
## Warning in merge.data.frame(toptab, tbl, by.x = 0, by.y = "Var1"): column names
## 'bl', 'gu' are duplicated in the result
toptab2 <- toptab2[order(toptab2$med),]

toptab2 <- subset(toptab2,Freq>=5)

head(toptab2,100) %>%
  kbl(caption="Top genes with reducing methylation") %>%
  kable_paper("hover", full_width = F)
Top genes with reducing methylation
Row.names bl gu bl.1 gu.1 med Freq
3995 CXADRP3 -2.2202814 -1.1159814 -18756 -15299 -1.668131 6
4102 CYP4Z2P -2.0461284 -1.1477867 -18732 -15316 -1.596958 6
10870 MIR1251 -1.7447979 -1.4227722 -18654 -15436 -1.583785 6
14974 POM121L8P -1.4626322 -1.5979785 -18426 -15472 -1.530305 5
1850 BZW1P2 -1.7535132 -1.2517218 -18657 -15385 -1.502617 5
14846 PMCHL2 -2.0334253 -0.9552235 -18730 -15107 -1.494324 8
12615 NALF1-IT1 -2.4591520 -0.4664148 -18774 -12679 -1.462783 7
14136 PCAT1 -2.2841133 -0.5947558 -18761 -13830 -1.439435 6
11530 MIR521-1 -2.1697975 -0.7089284 -18751 -14443 -1.439363 5
7041 GUSBP2 -1.6851804 -1.1895332 -18615 -15344 -1.437357 10
11533 MIR523 -1.1520646 -1.6403193 -17941 -15480 -1.396192 5
9520 LINC01069 -2.0860694 -0.6893484 -18738 -14352 -1.387709 5
13595 OR2J2 -1.6794291 -1.0927888 -18606 -15278 -1.386109 5
13488 OR10A5 -2.0015621 -0.7063150 -18721 -14434 -1.353938 6
13664 OR4M1 -2.1233483 -0.5783751 -18744 -13709 -1.350862 7
5494 FABP2 -1.8938632 -0.7898329 -18700 -14740 -1.341848 5
9421 LINC00871 -1.7116509 -0.9581566 -18636 -15109 -1.334904 9
3304 CLCA3P -2.2461886 -0.4172496 -18758 -11985 -1.331719 6
9345 LINC00613 -1.9736218 -0.6840457 -18716 -14331 -1.328834 6
9897 LINC02903 -1.7815191 -0.8509152 -18666 -14897 -1.316217 6
13761 OR5M1 -1.6720546 -0.9440198 -18604 -15085 -1.308037 6
8477 KERA -1.7025120 -0.8885475 -18630 -14995 -1.295530 5
533 AJM1 -1.3580746 -1.2324007 -18310 -15367 -1.295238 5
13707 OR52E8 -1.3805243 -1.2052246 -18330 -15351 -1.292874 8
16189 RHAG -1.7088007 -0.8673835 -18635 -14948 -1.288092 10
12587 NAALADL2-AS2 -1.5367760 -1.0166810 -18500 -15197 -1.276728 7
13756 OR5K2 -1.7772626 -0.7748558 -18665 -14697 -1.276059 5
19168 TAAR9 -1.4843207 -1.0655353 -18447 -15256 -1.274928 6
20844 UGT2B4 -1.6452967 -0.8965763 -18589 -15004 -1.270936 5
11184 MIR340 -1.4893008 -1.0337634 -18458 -15218 -1.261532 5
5705 FAM99A -0.7217411 -1.7931826 -16008 -15501 -1.257462 5
15560 PTH -1.7006418 -0.8074612 -18628 -14791 -1.254052 6
6747 GPC5-AS2 -1.6988305 -0.8025703 -18627 -14773 -1.250700 5
17059 SELE -1.4490469 -1.0517365 -18410 -15235 -1.250392 10
11471 MIR490 -1.5650594 -0.9179919 -18523 -15037 -1.241526 5
19166 TAAR6 -2.0587643 -0.3887210 -18736 -11494 -1.223743 9
18261 SNORD116-6 -1.2049191 -1.2419042 -18063 -15370 -1.223412 8
20029 TMPRSS11F -1.5533554 -0.8862582 -18511 -14993 -1.219807 13
479 AGMO -1.4975357 -0.9394015 -18469 -15078 -1.218469 24
9360 LINC00644 -1.4304641 -0.9998308 -18388 -15175 -1.215148 5
3090 CFHR3 -1.9095133 -0.5031942 -18706 -13080 -1.206354 6
2877 CDRT15 -1.2639027 -1.1294407 -18182 -15308 -1.196672 8
13748 OR5F1 -1.3528213 -1.0326622 -18301 -15217 -1.192742 6
18174 SNORD114-24 -1.9158299 -0.4693366 -18708 -12723 -1.192583 8
13623 OR2W1 -1.6068012 -0.7675571 -18562 -14665 -1.187179 7
8829 KRTAP22-1 -1.4951714 -0.8714871 -18465 -14957 -1.183329 7
9559 LINC01128 -0.8066052 -1.5509223 -16676 -15463 -1.178764 16
7656 HTN3 -1.9559818 -0.3975766 -18714 -11645 -1.176779 5
3829 CRYGC -1.3967692 -0.9518646 -18350 -15097 -1.174317 11
15339 PRPS1L1 -1.1868217 -1.1435632 -18026 -15315 -1.165192 12
7128 H4C11 -1.9141119 -0.4156714 -18707 -11961 -1.164892 11
19819 TMEM14EP -1.4951674 -0.8270972 -18464 -14839 -1.161132 5
1886 C12orf50 -1.5898104 -0.7218940 -18544 -14502 -1.155852 7
12305 MSTN -1.5756431 -0.7262507 -18532 -14516 -1.150947 7
5878 FER1L6-AS1 -1.4109479 -0.8899857 -18359 -14999 -1.150467 7
492 AGR3 -1.5219150 -0.7772743 -18491 -14705 -1.149595 5
17470 SLC15A5 -1.7078951 -0.5891900 -18634 -13781 -1.148543 10
11715 MIR643 -1.2281898 -1.0582404 -18115 -15244 -1.143215 6
18285 SNORD18C -1.4378113 -0.8428891 -18397 -14880 -1.140350 5
449 AFM -1.1436373 -1.1289105 -17918 -15307 -1.136274 6
18160 SNORD114-10 -1.4592133 -0.8027844 -18418 -14776 -1.130999 6
18991 STX19 -1.7012398 -0.5602416 -18629 -13572 -1.130741 10
9228 LINC00375 -1.2615999 -0.9910004 -18178 -15160 -1.126300 6
8671 KLRC2 -1.6223704 -0.6127489 -18573 -13945 -1.117560 5
6467 GFRAL -1.4819990 -0.7488058 -18445 -14602 -1.115402 10
15023 PP2D1 -1.3738722 -0.8535397 -18322 -14904 -1.113706 8
13778 OR6C1 -1.6603029 -0.5655805 -18597 -13618 -1.112942 5
13981 PABPC1P2 -0.8542082 -1.3689877 -16933 -15422 -1.111598 7
9322 LINC00575 -1.4873242 -0.7347831 -18450 -14550 -1.111054 6
8494 KIAA0408 -1.4994918 -0.7191214 -18472 -14488 -1.109307 8
19258 TAS2R39 -1.4712488 -0.7472461 -18435 -14596 -1.109247 5
13837 OR8K3 -1.4747777 -0.7408440 -18440 -14570 -1.107811 5
9260 LINC00448 -1.6053980 -0.6089826 -18561 -13913 -1.107190 6
18175 SNORD114-25 -1.6291811 -0.5809166 -18579 -13732 -1.105049 7
19254 TAS2R3 -1.5152879 -0.6916815 -18482 -14363 -1.103485 7
9361 LINC00645 -1.8557007 -0.3482450 -18689 -10718 -1.101973 6
20839 UGT2B10 -1.6194379 -0.5814921 -18572 -13737 -1.100465 5
9234 LINC00381 -1.6460939 -0.5542307 -18592 -13521 -1.100162 6
19406 TDO2 -1.6123622 -0.5843650 -18567 -13749 -1.098364 6
4507 DLGAP1-AS3 -1.2469504 -0.9490306 -18153 -15090 -1.097991 6
6642 GNAT3 -1.4264669 -0.7674098 -18381 -14664 -1.096938 13
8863 KRTAP6-2 -1.6839797 -0.5072401 -18610 -13119 -1.095610 7
11577 MIR548H2 -1.4099618 -0.7788018 -18357 -14713 -1.094382 30
21841 ZNF404 -1.2972700 -0.8886111 -18230 -14996 -1.092941 6
17936 SMIM23 -1.2314945 -0.9505836 -18122 -15093 -1.091039 5
2803 CDH19 -1.4879501 -0.6922320 -18452 -14365 -1.090091 18
18184 SNORD114-5 -1.5821343 -0.5943741 -18536 -13822 -1.088254 5
13799 OR6T1 -1.2934939 -0.8756839 -18226 -14968 -1.084589 8
17965 SMR3B -1.4944858 -0.6742780 -18463 -14291 -1.084382 8
22202 ZNG1E -1.5979788 -0.5679598 -18553 -13642 -1.082969 5
9886 LINC02877 -1.0617309 -1.0972361 -17693 -15282 -1.079484 6
13620 OR2T8 -1.4301558 -0.7256101 -18387 -14513 -1.077883 6
10880 MIR125B2 -1.6301545 -0.5213534 -18580 -13254 -1.075754 5
19266 TAS2R50 -1.5281616 -0.6180385 -18497 -13978 -1.073100 7
16288 RMST -1.3685723 -0.7709529 -18317 -14677 -1.069763 19
18849 ST7-OT3 -1.6173440 -0.5168714 -18570 -13209 -1.067108 6
19253 TAS2R20 -1.4654574 -0.6665000 -18429 -14250 -1.065979 6
19896 TMEM225 -0.9008958 -1.2287856 -17156 -15363 -1.064841 6
20157 TONSL-AS1 -1.0277125 -1.1009145 -17596 -15287 -1.064313 6
15588 PTPN20 -0.8769013 -1.2506030 -17052 -15383 -1.063752 18

Now look at the pathways with hypomethylation.

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

par(mar=c(9.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), cex.names=0.8,
  xlab="enrichment score")

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

pdf("adosbarplot_hypo_go.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 bitter taste receptors.

genesets <- go

bitter <- genesets[grep("GOMF BITTER",names(genesets))]
bitter
## $`GOMF BITTER TASTE RECEPTOR ACTIVITY`
##  [1] "TAS2R1"  "TAS2R10" "TAS2R13" "TAS2R14" "TAS2R16" "TAS2R19" "TAS2R20"
##  [8] "TAS2R3"  "TAS2R30" "TAS2R31" "TAS2R38" "TAS2R39" "TAS2R4"  "TAS2R40"
## [15] "TAS2R41" "TAS2R42" "TAS2R43" "TAS2R45" "TAS2R46" "TAS2R5"  "TAS2R50"
## [22] "TAS2R60" "TAS2R7"  "TAS2R8"  "TAS2R9"
rnk2 <- rnk[which(rownames(rnk) %in% bitter[[1]]),]

plot(rnk2,main=names(bitter),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("bitter.pdf")

plot(rnk2,main=names(bitter),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

Heatmap

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

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

bitter_probes <- gt[which(gt$gene %in% bitter[[1]]),"probe"]

bitter_mvals <- bl_mvals[rownames(bl_mvals) %in% bitter_probes,]

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

heatmap.2(bitter_mvals,trace="none",col=colfunc(25), margins = c(5,8), scale="row",
  main="Bitter taste receptors",
  ColSideColors = colcols )

Now with reactome

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

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_hypo_reactome.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.

genesets <- reactome

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 <- gt[which(gt$gene %in% phero[[1]]),"probe"]

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 )

SFARI

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        "1061"        
## pMANOVA        "1.008455e-05"
## s.bl           "-0.08046307" 
## s.gu           "-0.04932763" 
## p.bl           "9.386797e-06"
## p.gu           "0.006608851" 
## s.dist         "0.09437967"  
## SD             "0.02201608"  
## p.adjustMANOVA "1.008455e-05"
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=1,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 1061 0.0000101 -0.0804631 -0.0493276 0.0000094 0.0066089 0.0943797 0.0220161 0.0104375
664 random 663 1000 0.0006176 -0.0618008 -0.0485970 0.0009386 0.0092852 0.0786194 0.0093365 0.3196333
211 random 210 1000 0.0024331 -0.0513177 -0.0493682 0.0060138 0.0082259 0.0712090 0.0013785 0.5460660
328 random 327 1000 0.0038395 -0.0618170 -0.0204526 0.0009357 0.2736267 0.0651126 0.0292490 0.5460660
929 random 928 1000 0.0082477 0.0417376 0.0478681 0.0254736 0.0103968 0.0635089 0.0043349 0.5460660
286 random 285 1000 0.0092805 0.0400016 0.0482284 0.0322570 0.0098334 0.0626586 0.0058172 0.5460660 
mitch_plots(res=sfarires2,outfile="sfari_mitch_ados.pdf")
## png 
##   2
mitch_report(res=sfarires2,outfile="sfari_mitch_ados.html",overwrite=TRUE)
## Dataset saved as " /tmp/RtmpGyw2UW/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: /asd_meth/mitch.knit.md
## /usr/bin/pandoc +RTS -K512m -RTS /asd_meth/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpGyw2UW/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 --embed-resources --standalone --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/RtmpGyw2UW/rmarkdown-str1ce511c332a.html
## 
## Output created: /tmp/RtmpGyw2UW/mitch_report.html
## [1] TRUE

Buccal pathway analysis

GO

bumg <- mitch_calc(bum, go, priority="effect",cores=8)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
bumgres <- bumg$enrichment_result

head(subset(bumgres,s.dist<0 ),10) %>%
  kbl(caption="Buccal: Top gene sets with reducing methylation") %>%
  kable_paper("hover", full_width = F)
Buccal: Top gene sets with reducing methylation
set setSize pANOVA s.dist p.adjustANOVA
6561 GOMF TYPE I INTERFERON RECEPTOR BINDING 14 0.0000000 -0.8664242 0.0000051
5599 GOMF BITTER TASTE RECEPTOR ACTIVITY 24 0.0000000 -0.7573668 0.0000000
3301 GOBP POSITIVE REGULATION OF PEPTIDYL SERINE PHOSPHORYLATION OF STAT PROTEIN 17 0.0000003 -0.7138465 0.0000698
5024 GOBP SERINE PHOSPHORYLATION OF STAT PROTEIN 21 0.0000014 -0.6078321 0.0002238
1993 GOBP NATURAL KILLER CELL CHEMOTAXIS 10 0.0011721 -0.5926558 0.0420861
6498 GOMF TASTE RECEPTOR ACTIVITY 30 0.0000000 -0.5788372 0.0000098
5367 GOBP T HELPER 2 CELL CYTOKINE PRODUCTION 13 0.0009465 -0.5294959 0.0360945
389 GOBP CD4 POSITIVE ALPHA BETA T CELL CYTOKINE PRODUCTION 19 0.0001190 -0.5098671 0.0073503
4317 GOBP REGULATION OF NATURAL KILLER CELL PROLIFERATION 12 0.0030836 -0.4933256 0.0851371
3153 GOBP POSITIVE REGULATION OF ISOTYPE SWITCHING TO IGG ISOTYPES 11 0.0071092 -0.4686665 0.1452383 
head(subset(bumgres,s.dist>0 ),10) %>%
  kbl(caption="Buccal: Top gene sets with higher methylation") %>%
  kable_paper("hover", full_width = F)
Buccal: Top gene sets with higher methylation
set setSize pANOVA s.dist p.adjustANOVA
3240 GOBP POSITIVE REGULATION OF MUSCLE ADAPTATION 13 0.0007494 0.5398844 0.0300134
6485 GOMF SULFURTRANSFERASE ACTIVITY 11 0.0020579 0.5365697 0.0625492
5872 GOMF HISTONE CHAPERONE ACTIVITY 12 0.0019201 0.5171650 0.0600680
3604 GOBP PROTEIN LOCALIZATION TO CENP A CONTAINING CHROMATIN 16 0.0004605 0.5057140 0.0205943
6781 GOCC CYTOSOLIC SMALL RIBOSOMAL SUBUNIT 37 0.0000001 0.4992402 0.0000324
6709 GOCC CHROMOSOME CENTROMERIC CORE DOMAIN 17 0.0005815 0.4818542 0.0249346
6470 GOMF STRUCTURAL CONSTITUENT OF CHROMATIN 80 0.0000000 0.4683464 0.0000000
6385 GOMF RIBOSOMAL SMALL SUBUNIT BINDING 16 0.0012137 0.4671643 0.0433652
3413 GOBP POSITIVE REGULATION OF STRIATED MUSCLE CONTRACTION 13 0.0042987 0.4573423 0.1065751
2802 GOBP PHOSPHATIDYLETHANOLAMINE BIOSYNTHETIC PROCESS 10 0.0137988 0.4496768 0.2220303 
mysets <- c("GOMF BITTER TASTE RECEPTOR ACTIVITY",
"GOMF TYPE I INTERFERON RECEPTOR BINDING",
"GOBP CELLULAR GLUCURONIDATION",
"GOBP POSITIVE REGULATION OF PEPTIDYL SERINE PHOSPHORYLATION OF STAT PROTEIN",
"GOMF ODORANT BINDING",
"GOMF OLFACTORY RECEPTOR ACTIVITY",
"GOBP URONIC ACID METABOLIC PROCESS")

bumgres[which(bumgres$set %in% mysets),] %>%
  kbl(caption="Selected gene sets in buccal samples") %>%
  kable_paper("hover", full_width = F)
Selected gene sets in buccal samples
set setSize pANOVA s.dist p.adjustANOVA
6561 GOMF TYPE I INTERFERON RECEPTOR BINDING 14 0.0000000 -0.8664242 0.0000051
5599 GOMF BITTER TASTE RECEPTOR ACTIVITY 24 0.0000000 -0.7573668 0.0000000
3301 GOBP POSITIVE REGULATION OF PEPTIDYL SERINE PHOSPHORYLATION OF STAT PROTEIN 17 0.0000003 -0.7138465 0.0000698
6149 GOMF ODORANT BINDING 104 0.0000000 -0.4472935 0.0000000
6150 GOMF OLFACTORY RECEPTOR ACTIVITY 356 0.0000000 -0.4388773 0.0000000
406 GOBP CELLULAR GLUCURONIDATION 21 0.0159259 -0.3038487 0.2422823
5386 GOBP URONIC ACID METABOLIC PROCESS 26 0.0672663 -0.2073191 0.4698194 
MAX=max(bumg$ranked_profile)
MIN=min(bumg$ranked_profile)

par(mfrow=c(3,1))
par(mar=c(4.1, 4.1, 5.1, 1.1))

beeswarm(bumg$detailed_sets["GOMF BITTER TASTE RECEPTOR ACTIVITY"],
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
  abline(v=c(MIN,MAX),lty=2)
mtext("GOMF BITTER TASTE RECEPTOR ACTIVITY (FDR=4.8e-8)",cex=0.7)

beeswarm(bumg$detailed_sets["GOMF TYPE I INTERFERON RECEPTOR BINDING"],
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
  abline(v=c(MIN,MAX),lty=2)
mtext("GOMF TYPE I INTERFERON RECEPTOR BINDING (FDR=5.1e-6)",cex=0.7)

mygenes <- go[[which(names(go)=="GOBP CELLULAR GLUCURONIDATION")]]
myranks <- bumg$ranked_profile[which(rownames(bumg$ranked_profile) %in% mygenes),]

beeswarm(myranks,
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
mtext("GOBP CELLULAR GLUCURONIDATION (FDR=0.24)",cex=0.7)
  abline(v=c(MIN,MAX),lty=2)

pdf("buccal_go.pdf")
par(mfrow=c(3,1))
par(mar=c(4.1, 4.1, 5.1, 1.1))

beeswarm(bumg$detailed_sets["GOMF BITTER TASTE RECEPTOR ACTIVITY"],
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
  abline(v=c(MIN,MAX),lty=2)
mtext("GOMF BITTER TASTE RECEPTOR ACTIVITY (FDR=4.8e-8)",cex=0.7)

beeswarm(bumg$detailed_sets["GOMF TYPE I INTERFERON RECEPTOR BINDING"],
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
  abline(v=c(MIN,MAX),lty=2)
mtext("GOMF TYPE I INTERFERON RECEPTOR BINDING (FDR=5.1e-6)",cex=0.7)

mygenes <- go[[which(names(go)=="GOBP CELLULAR GLUCURONIDATION")]]
myranks <- bumg$ranked_profile[which(rownames(bumg$ranked_profile) %in% mygenes),]

beeswarm(myranks,
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
mtext("GOBP CELLULAR GLUCURONIDATION (FDR=0.24)",cex=0.7)
  abline(v=c(MIN,MAX),lty=2)
dev.off()
## png 
##   2

REACTOME

bumr <- mitch_calc(bum, reactome, priority="effect",cores=8)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
bumrres <- bumr$enrichment_result

head(subset(bumrres,s.dist<0 ),10) %>%
  kbl(caption="Buccal: Top gene sets with reducing methylation") %>%
  kable_paper("hover", full_width = F)
Buccal: Top gene sets with reducing methylation
set setSize pANOVA s.dist p.adjustANOVA
1460 TRAF6 mediated IRF7 activation 26 0.0000226 -0.4800332 0.0008008
220 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000003 -0.4747295 0.0000142
642 Interaction With Cumulus Cells And The Zona Pellucida 11 0.0107594 -0.4440731 0.0913985
412 Expression and translocation of olfactory receptors 357 0.0000000 -0.4376575 0.0000000
895 Olfactory Signaling Pathway 364 0.0000000 -0.4280739 0.0000000
816 NF-kB activation through FADD/RIP-1 pathway mediated by caspase-8 and -10 12 0.0136817 -0.4110248 0.1057573
674 Intrinsic Pathway of Fibrin Clot Formation 21 0.0023778 -0.3829735 0.0313148
179 Caspase activation via Death Receptors in the presence of ligand 16 0.0104302 -0.3698080 0.0900090
1282 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000461 -0.3676606 0.0013685
656 Interleukin-15 signaling 13 0.0217792 -0.3674715 0.1455698 
head(subset(bumrres,s.dist>0 ),10) %>%
  kbl(caption="Buccal: Top gene sets with higher methylation") %>%
  kable_paper("hover", full_width = F)
Buccal: Top gene sets with higher methylation
set setSize pANOVA s.dist p.adjustANOVA
681 Josephin domain DUBs 10 0.0009048 0.6059795 0.0162170
954 Pexophagy 10 0.0018332 0.5689980 0.0282078
214 Chylomicron remodeling 10 0.0020745 0.5623092 0.0296795
1461 TRAF6 mediated IRF7 activation in TLR7/8 or 9 signaling 11 0.0013317 0.5587307 0.0221631
1228 SARS-CoV-1 modulates host translation machinery 34 0.0000001 0.5275375 0.0000055
1265 SUMOylation of immune response proteins 11 0.0073661 0.4666014 0.0737059
957 Phase 4 - resting membrane potential 19 0.0009019 0.4398319 0.0162170
1072 RNA Polymerase I Promoter Opening 26 0.0001435 0.4307122 0.0035451
1130 Regulation of CDH11 gene transcription 10 0.0192168 0.4275543 0.1357350
817 NF-kB is activated and signals survival 12 0.0112359 0.4226527 0.0930244 
mysets <- c("Class C/3 (Metabotropic glutamate/pheromone receptors)",
"Scavenging of heme from plasma",
"Sensory perception of sweet, bitter, and umami (glutamate) taste",
"Glucuronidation",
"Sensory perception of taste",
"Digestion",
"Aspirin ADME",
"Expression and translocation of olfactory receptors",
"Olfactory Signaling Pathway",
"Dectin-2 family")

bumrres[which(bumrres$set %in% mysets),] %>%
  kbl(caption="Selected gene sets in buccal samples") %>%
  kable_paper("hover", full_width = F)
Selected gene sets in buccal samples
set setSize pANOVA s.dist p.adjustANOVA
220 Class C/3 (Metabotropic glutamate/pheromone receptors) 39 0.0000003 -0.4747295 0.0000142
412 Expression and translocation of olfactory receptors 357 0.0000000 -0.4376575 0.0000000
895 Olfactory Signaling Pathway 364 0.0000000 -0.4280739 0.0000000
1282 Sensory perception of sweet, bitter, and umami (glutamate) taste 41 0.0000461 -0.3676606 0.0013685
1283 Sensory perception of taste 47 0.0014228 -0.2689132 0.0229762
540 Glucuronidation 24 0.0360649 -0.2471738 0.1966605
100 Aspirin ADME 43 0.0234052 -0.1997691 0.1514837
288 Dectin-2 family 28 0.1643760 -0.1518170 0.4865649
1271 Scavenging of heme from plasma 12 0.4337351 -0.1305109 0.7357329
328 Digestion 18 0.8381293 -0.0278127 0.9414524 
MAX=max(bumr$ranked_profile)
MIN=min(bumr$ranked_profile)

par(mfrow=c(3,1))
par(mar=c(4.1, 4.1, 5.1, 1.1))

beeswarm(bumr$detailed_sets["Class C/3 (Metabotropic glutamate/pheromone receptors)"],
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
  abline(v=c(MIN,MAX),lty=2)
mtext("Class C/3 (Metabotropic glutamate/pheromone receptors) (FDR=1.4e-5)",cex=0.7)

mygenes <- reactome[[which(names(reactome)=="Expression and translocation of olfactory receptors")]]
myranks <- bumr$ranked_profile[which(rownames(bumr$ranked_profile) %in% mygenes),]
beeswarm(myranks,
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
mtext("Expression and translocation of olfactory receptors (FDR=8.5e-43)",cex=0.7)
  abline(v=c(MIN,MAX),lty=2)

mygenes <- reactome[[which(names(reactome)=="Olfactory Signaling Pathway")]]
myranks <- bumr$ranked_profile[which(rownames(bumr$ranked_profile) %in% mygenes),]
beeswarm(myranks,
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
mtext("Olfactory Signaling Pathway (FDR=5.4e-42)",cex=0.7)
  abline(v=c(MIN,MAX),lty=2)

pdf("buccal_reactome.pdf")
par(mfrow=c(3,1))
par(mar=c(4.1, 4.1, 5.1, 1.1))

beeswarm(bumr$detailed_sets["Class C/3 (Metabotropic glutamate/pheromone receptors)"],
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
  abline(v=c(MIN,MAX),lty=2)
mtext("Class C/3 (Metabotropic glutamate/pheromone receptors) (FDR=1.4e-5)",cex=0.7)

mygenes <- reactome[[which(names(reactome)=="Expression and translocation of olfactory receptors")]]
myranks <- bumr$ranked_profile[which(rownames(bumr$ranked_profile) %in% mygenes),]
beeswarm(myranks,
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
mtext("Expression and translocation of olfactory receptors (FDR=8.5e-43)",cex=0.7)
  abline(v=c(MIN,MAX),lty=2)

mygenes <- reactome[[which(names(reactome)=="Olfactory Signaling Pathway")]]
myranks <- bumr$ranked_profile[which(rownames(bumr$ranked_profile) %in% mygenes),]
beeswarm(myranks,
  horizontal=TRUE,labels="", xlim=c(MIN,MAX), bty="n",
  xlab="ranked score in buccal swab")
mtext("Olfactory Signaling Pathway (FDR=5.4e-42)",cex=0.7)
  abline(v=c(MIN,MAX),lty=2)

dev.off()
## png 
##   2

Session information

For reproducibility

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 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] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] pkgload_1.3.4                                      
##  [2] GGally_2.2.1                                       
##  [3] ggplot2_3.5.1                                      
##  [4] reshape2_1.4.4                                     
##  [5] beeswarm_0.4.0                                     
##  [6] gplots_3.1.3.1                                     
##  [7] gtools_3.9.5                                       
##  [8] tibble_3.2.1                                       
##  [9] dplyr_1.1.4                                        
## [10] echarts4r_0.4.5                                    
## [11] RhpcBLASctl_0.23-42                                
## [12] eulerr_7.0.2                                       
## [13] kableExtra_1.4.0                                   
## [14] data.table_1.15.4                                  
## [15] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [16] minfi_1.50.0                                       
## [17] bumphunter_1.46.0                                  
## [18] locfit_1.5-9.10                                    
## [19] iterators_1.0.14                                   
## [20] foreach_1.5.2                                      
## [21] Biostrings_2.72.1                                  
## [22] XVector_0.44.0                                     
## [23] SummarizedExperiment_1.34.0                        
## [24] Biobase_2.64.0                                     
## [25] MatrixGenerics_1.16.0                              
## [26] matrixStats_1.3.0                                  
## [27] GenomicRanges_1.56.1                               
## [28] GenomeInfoDb_1.40.1                                
## [29] IRanges_2.38.1                                     
## [30] S4Vectors_0.42.1                                   
## [31] BiocGenerics_0.50.0                                
## [32] HGNChelper_0.8.14                                  
## [33] mitch_1.16.0                                       
## [34] limma_3.60.4                                       
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.4.1             later_1.3.2              
##   [3] BiocIO_1.14.0             bitops_1.0-8             
##   [5] preprocessCore_1.61.0     XML_3.99-0.17            
##   [7] lifecycle_1.0.4           lattice_0.22-6           
##   [9] MASS_7.3-61               base64_2.0.1             
##  [11] scrime_1.3.5              magrittr_2.0.3           
##  [13] sass_0.4.9                rmarkdown_2.27           
##  [15] jquerylib_0.1.4           yaml_2.3.8               
##  [17] httpuv_1.6.15             doRNG_1.8.6              
##  [19] askpass_1.2.0             DBI_1.2.3                
##  [21] RColorBrewer_1.1-3        abind_1.4-5              
##  [23] zlibbioc_1.50.0           quadprog_1.5-8           
##  [25] purrr_1.0.2               RCurl_1.98-1.16          
##  [27] GenomeInfoDbData_1.2.12   genefilter_1.86.0        
##  [29] annotate_1.82.0           svglite_2.1.3            
##  [31] DelayedMatrixStats_1.26.0 codetools_0.2-20         
##  [33] DelayedArray_0.30.1       xml2_1.3.6               
##  [35] tidyselect_1.2.1          farver_2.1.2             
##  [37] UCSC.utils_1.0.0          beanplot_1.3.1           
##  [39] illuminaio_0.46.0         GenomicAlignments_1.40.0 
##  [41] jsonlite_1.8.8            multtest_2.60.0          
##  [43] survival_3.7-0            systemfonts_1.1.0        
##  [45] tools_4.4.1               Rcpp_1.0.12              
##  [47] glue_1.7.0                gridExtra_2.3            
##  [49] SparseArray_1.4.8         xfun_0.45                
##  [51] HDF5Array_1.32.1          withr_3.0.0              
##  [53] fastmap_1.2.0             rhdf5filters_1.16.0      
##  [55] fansi_1.0.6               openssl_2.2.0            
##  [57] caTools_1.18.2            digest_0.6.36            
##  [59] R6_2.5.1                  mime_0.12                
##  [61] colorspace_2.1-1          RSQLite_2.3.7            
##  [63] utf8_1.2.4                tidyr_1.3.1              
##  [65] generics_0.1.3            rtracklayer_1.64.0       
##  [67] httr_1.4.7                htmlwidgets_1.6.4        
##  [69] S4Arrays_1.4.1            ggstats_0.6.0            
##  [71] pkgconfig_2.0.3           gtable_0.3.5             
##  [73] blob_1.2.4                siggenes_1.78.0          
##  [75] htmltools_0.5.8.1         scales_1.3.0             
##  [77] png_0.1-8                 knitr_1.47               
##  [79] rstudioapi_0.16.0         tzdb_0.4.0               
##  [81] rjson_0.2.22              nlme_3.1-165             
##  [83] curl_5.2.1                cachem_1.1.0             
##  [85] rhdf5_2.48.0              stringr_1.5.1            
##  [87] KernSmooth_2.23-24        AnnotationDbi_1.66.0     
##  [89] restfulr_0.0.15           GEOquery_2.72.0          
##  [91] pillar_1.9.0              grid_4.4.1               
##  [93] reshape_0.8.9             vctrs_0.6.5              
##  [95] promises_1.3.0            xtable_1.8-4             
##  [97] evaluate_0.24.0           readr_2.1.5              
##  [99] GenomicFeatures_1.56.0    cli_3.6.3                
## [101] compiler_4.4.1            Rsamtools_2.20.0         
## [103] rlang_1.1.4               crayon_1.5.3             
## [105] rngtools_1.5.2            labeling_0.4.3           
## [107] nor1mix_1.3-3             mclust_6.1.1             
## [109] plyr_1.8.9                stringi_1.8.4            
## [111] viridisLite_0.4.2         BiocParallel_1.38.0      
## [113] munsell_0.5.1             Matrix_1.7-0             
## [115] hms_1.1.3                 sparseMatrixStats_1.16.0 
## [117] bit64_4.0.5               Rhdf5lib_1.26.0          
## [119] KEGGREST_1.44.1           statmod_1.5.0            
## [121] shiny_1.8.1.1             highr_0.11               
## [123] memoise_2.0.1             bslib_0.7.0              
## [125] bit_4.0.5                 splitstackshape_1.4.8
save.image("downstream2_ados.Rdata")