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:

  • ASD_blood_top_dmps_onADOS_withintw_Nov27_limma.csv

  • ASD_guthrie_top_dmps_onADOS_withintw_Nov27_limma.csv

suppressPackageStartupMessages({
  library("parallel")
  library("mitch")
  library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
  source("https://raw.githubusercontent.com/markziemann/ART_methylation/master/meth_functions.R")
  library("data.table")
  library("kableExtra")
  library("eulerr")
  library("tictoc")
})

Load data

bl_design <- readRDS(file="bl_design_ados.rds")
bl_mvals <- readRDS(file="bl_mvals.rds")

gu_design <- readRDS(file="gu_design_ados.rds")
gu_mvals <- readRDS(file="gu_mvals.rds")

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 <- promoters[,"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    8.00   15.00   21.27   26.00  536.00

GMEA functions

This is an enrichment technique.

# get the median, mean and 1-sample t-test result for each gene
pmea <- function(mval,design,sets,cores=2) {
  fit <- lmFit(mval, design)
  fit <- eBayes(fit)
  top <- topTable(fit,coef=ncol(design),num=Inf, sort.by = "P")
  l <- mclapply(seq(1,length(sets)), function(i) {
    g <- names(sets[i])
    tstats <- top[rownames(top) %in% sets[[i]],"t"]
    myn <- length(tstats)
    mymean <- mean(tstats)
    mymedian <- median(tstats)
    if ( length(tstats) < 2 ) {
      pval=1
    } else {
      wtselfcont <- t.test(tstats)
      pval=wtselfcont$p.value
    }
    c("gene"=g,"nprobes"=myn,"mean"=mymean,"median"=mymedian,
      "P.Value"=pval)
  } , mc.cores=cores)

  df <- do.call(rbind, l)
  rownames(df) <- df[,1]
  df <- df[,-1]
  tmp <- apply(df,2,as.numeric)
  rownames(tmp) <- rownames(df)
  df <- as.data.frame(tmp)
  df$sig <- -log10(df[,4])
  df <- df[order(-df$sig),]
  df$FDR <- p.adjust(df$P.Value)
  out <- list("df"=df,"toptable"=top)
  return(out)
}
# pmea_res <- pmea(mval=mval,design=design,sets=head(sets,20),cores=detectCores()/2)

# Run the fry test for each gene. This is a more conservative test.
run_fry <- function(mval,design,sets,cores=2) {
  split_sets <- split(sets, ceiling(seq_along(sets)/200))
  fry_l <- mclapply(split_sets,function(l) {
    fry(y=mval, index = l, design = design,
      contrast = ncol(design) )
  } , mc.cores=cores )
  fry_res <- do.call(rbind,fry_l)
  rownames(fry_res) <- sub("\\.","@",rownames(fry_res))
  rownames(fry_res) <- sapply(strsplit(rownames(fry_res),"@"),"[[",2)
  fry_res[is.na(fry_res$PValue),"PValue"] <- 1
  fry_res <- fry_res[order(fry_res$PValue),]
  fry_res$FDR <- p.adjust(fry_res$PValue,method="fdr")
  return(fry_res)
}
#fry_res <- run_fry(mval=mval,design=design,sets=sets,cores=cores)

# main function to perform 1-sample t-test and fry test and merge the results.
main <- function(mval,design,sets,cores=2){
  pmea_res <- pmea(mval=mval,design=design,sets=sets,cores=cores)
  pmea_df <- pmea_res[[1]]
  limma_df <- pmea_res[[2]]
  fry_res <- run_fry(mval=mval,design=design,sets=sets,cores=cores)
  m <- merge(pmea_df,fry_res,by=0)
  rownames(m) <- m$Row.names
  m$Row.names = NULL
  m <- m[,c("nprobes","median","PValue","FDR.y")]
  colnames(m) <- c("nprobes","median","PValue","FDR")
  m <- m[order(m$PValue),]
  out <- list("res"=m,"limma_df"=limma_df)
  return(out)
}
#res <- main(mval,design,sets,cores=detectCores()/2)

probe_bias <- function(res) {
  res$sig <- -log10(res$PValue)
  sig <- subset(res,FDR < 0.05)
  plot(res$nprobes,res$sig,log="x",
    pch=19,cex=0.6,
    xlab="no. probes",ylab="-log10(p-value)")
  points(sig$nprobes,sig$sig,col="red",pch=19,cex=0.62)
  SIG = nrow(sig)
  TOT = nrow(res)
  HEADER <- paste(TOT, "genes in total.", SIG, "with FDR<0.05.")
  mtext(HEADER)
}

volcano_plot <- function(res) {
  res$sig <- -log10(res$PValue)
  sig <- subset(res,FDR < 0.05)
  plot(res$median,res$sig,
    pch=19,cex=0.6,
    xlab="median t-statistic",ylab="-log10(p-value)")
  points(sig$median,sig$sig,col="red",pch=19,cex=0.62)
  SIG = nrow(sig)
  UP = nrow(subset(sig,median>0))
  DN = nrow(subset(sig,median<0))
  TOT = nrow(res)
  HEADER <- paste(TOT, "genes in total.", SIG, "with FDR<0.05;",DN,"down,",UP,"up")
  mtext(HEADER)
}

gmea_boxplot <- function(res,sets,n=50) {
  df <- res[[1]]
  limma_df <- res[[2]]
  # smallest pval
  par(mfrow=c(1,2))
  gs <- head(rownames(df),n)
  mysets <- sets[names(sets) %in% gs]
  tstats <- lapply(mysets, function(set) {
    limma_df[rownames(limma_df) %in% set,"t"]
  })
  tstats <- tstats[order(unlist(lapply(tstats,median)))]
  boxplot(tstats,horizontal=TRUE,las=1,
    main="smallest p-val",cex.axis=0.6,
    xlab="t-statistic")
  grid()
  # biggest effect size (median)
  sig <- subset(df,FDR < 0.05)
  gs <- head(rownames(sig[order(-abs(sig$median)),]),n)
  if ( length(gs) >2 ) {
    tstats <- lapply(gs, function(g) {
      df[which(df$genes==g),"tvals"]
    })
    names(tstats) <- gs
    tstats <- tstats[order(unlist(lapply(tstats,median)))]
    boxplot(tstats,horizontal=TRUE,las=1,
      main="biggest effect size(median)",cex.axis=0.6,
      xlab="t-statistic")
    grid()
  } else {
    plot(1)
    mtext("too few significant genes found")
  }
    par(mfrow=c(1,1))
}

GMEA for blood

Effect size is estimated with the median t-statistic for each gene. Significance is estimated using the fry test. Results are merged for the user.

  • PMEA step 196.063 sec elapsed on 8 cores

  • Fry step 293.569 sec elapsed on 8 cores

  • Fry single core 1205.519 sec

bl_design <- readRDS(file="bl_design_ados.rds")
bl_mvals <- readRDS(file="bl_mvals.rds")

tic()
res_bl <- main(mval=bl_mvals,design=bl_design,sets,cores=detectCores()/2)
res_bl_df <- res_bl[[1]] # enrichment results
res_bl_top <- res_bl[[2]] # limma results
toc() ## 502-522 sec elapsed on 8 cores (8.3-8.7 mins)
## 238.275 sec elapsed
head(res_bl_df,30) %>%
  kbl(caption="Top DM genes in blood identified with GMEA") %>%
  kable_paper("hover", full_width = F)
Top DM genes in blood identified with GMEA
nprobes median PValue FDR
DALRD3 25 -0.6269258 0.0009725 1
CIZ1 18 0.8751525 0.0011199 1
PDIA6 12 -0.9445400 0.0015166 1
TMEM49 16 0.8937461 0.0015360 1
ATP2B4 24 0.3807179 0.0019190 1
SLC19A1 6 -1.2469666 0.0022313 1
C22orf25 7 -0.7343741 0.0023912 1
CMTR2 2 -1.3024494 0.0024033 1
SRSF11 3 1.5974439 0.0028130 1
ZNF713 3 2.9556156 0.0028865 1
ZNF169 7 0.8634934 0.0029543 1
JAM2 13 -0.6277891 0.0029771 1
MBD3 7 -0.9816297 0.0029825 1
ST3GAL4-AS1 5 -1.0745998 0.0032239 1
C19orf66 11 -0.6774133 0.0034851 1
BAIAP2 13 -0.7621523 0.0037113 1
TBC1D9B 6 -0.8557415 0.0038188 1
NLE1 10 0.5536044 0.0040251 1
CYHR1 22 -0.7675003 0.0041759 1
PITPNA 8 -0.5708743 0.0044339 1
C4orf37 8 -1.1385454 0.0044778 1
ZNF589 13 -0.3344240 0.0045633 1
UFC1 12 1.0438509 0.0045835 1
PBXIP1 11 -0.4188172 0.0047662 1
NDUFAF3 18 -0.5845658 0.0049528 1
RANBP1 10 -0.3142372 0.0052860 1
TRMT2A 10 -0.3142372 0.0052860 1
FAM194A 8 1.4169134 0.0053030 1
FLJ90757 12 -0.7390261 0.0053785 1
ZNF506 5 -1.6560745 0.0058725 1

Now some visualisations.

volcano_plot(res_bl_df)

probe_bias(res_bl_df)
## Warning in xy.coords(x, y, xlabel, ylabel, log): 412 x values <= 0 omitted from
## logarithmic plot

gmea_boxplot(res_bl,sets)

GMEA for guthrie

gu_design <- readRDS(file="gu_design_ados.rds")
gu_mvals <- readRDS(file="gu_mvals.rds")

tic()
res_gu <- main(mval=gu_mvals,design=gu_design,sets,cores=detectCores()/2)
res_gu_df <- res_gu[[1]] # enrichment results
res_gu_top <- res_gu[[2]] # limma results
toc() ## 502-522 sec elapsed on 8 cores (8.3-8.7 mins)
## 220.497 sec elapsed
head(res_gu_df,30) %>%
  kbl(caption="Top DM genes in Guthrie card identified with GMEA") %>%
  kable_paper("hover", full_width = F)
Top DM genes in Guthrie card identified with GMEA
nprobes median PValue FDR
SRGAP3 16 0.9290154 0.0001047 1
GPAM 11 0.6168445 0.0003158 1
ATMIN 10 1.1732044 0.0004538 1
TMEM63A 9 0.5676956 0.0008588 1
ARID1A 1 -4.2765420 0.0009078 1
MIPEP 15 0.4143746 0.0009778 1
BCR 6 0.8488818 0.0011635 1
ERLIN1 14 0.9973830 0.0013343 1
RAD23B 12 0.9244582 0.0013689 1
ZNF600 4 1.0510153 0.0015567 1
BEST4 6 0.5148220 0.0018416 1
SNORD44 1 -4.0993344 0.0018628 1
C1orf124 13 0.7695296 0.0020870 1
FOXP4-AS1 1 1.9057650 0.0021012 1
GNAI3 16 0.9329586 0.0022403 1
CPD 4 1.4684642 0.0022877 1
TBCCD1 29 0.4670442 0.0023720 1
SNX20 12 -0.8902897 0.0024846 1
ZNF337 12 0.6011552 0.0025426 1
GOLGA2 9 1.0107937 0.0027118 1
ACAD8 15 0.2787049 0.0028382 1
THYN1 15 0.2787049 0.0028382 1
ZNRF1 13 0.5628678 0.0028526 1
SDHAF4 1 -3.9549955 0.0028672 1
NAF1 2 1.2461706 0.0032179 1
TPK1 17 -0.2004702 0.0032548 1
SDE2 3 1.4131508 0.0034158 1
LRRC58 12 0.4902156 0.0034537 1
HNRPLL 7 1.1035547 0.0035194 1
LOC101927478 2 1.5358760 0.0035301 1

Now some visualisations.

volcano_plot(res_gu_df)

probe_bias(res_gu_df)
## Warning in xy.coords(x, y, xlabel, ylabel, log): 421 x values <= 0 omitted from
## logarithmic plot

gmea_boxplot(res_gu,sets)

Mitch enrichment analysis

Reactome was downloaded December 2022.

genesets <- gmt_import("ReactomePathways.gmt")

res_bl_df$bl <- res_bl_df$median * -log10(res_bl_df$PValue)
res_gu_df$gu <- res_gu_df$median * -log10(res_gu_df$PValue)

m <- merge(res_bl_df,res_gu_df,by=0)
rownames(m) <- m$Row.names
m2 <- m[,c("bl","gu")]
m2[is.na(m2)] <- 0

res <- mitch_calc(m2, genesets, priority="effect")
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
#mitch_report(res,"blgu_mitch.html",overwrite=TRUE)

head(res$enrichment_result,30)  %>%
  kbl(caption="Top pathways found with GMEA/fry+mitch (effect size)") %>%
  kable_paper("hover", full_width = F)
Top pathways found with GMEA/fry+mitch (effect size)
set setSize pMANOVA s.bl s.gu p.bl p.gu s.dist SD p.adjustMANOVA
925 SARS-CoV-1 targets host intracellular signalling and regulatory pathways 15 0.0009625 -0.4728207 -0.3389279 0.0015231 0.0230699 0.5817487 0.0946765 0.1196350
703 PTK6 Regulates RHO GTPases, RAS GTPase and MAP kinases 10 0.0179034 -0.3974593 -0.3710357 0.0295459 0.0422113 0.5437291 0.0186843 0.2298180
931 SARS-CoV-2 targets host intracellular signalling and regulatory pathways 11 0.0213227 -0.4151958 -0.2880242 0.0171192 0.0981689 0.5053172 0.0899239 0.2398571
982 Signaling by Activin 12 0.0156837 -0.3149874 -0.3932787 0.0588931 0.0183438 0.5038702 0.0553604 0.2190594
598 Molecules associated with elastic fibres 11 0.0243138 -0.4597960 0.0709336 0.0082830 0.6838097 0.4652353 0.3752825 0.2398571
1078 Synthesis of PIPs at the early endosome membrane 13 0.0221373 0.4011543 0.2261147 0.0122796 0.1581596 0.4604918 0.1237717 0.2398571
561 Metabolism of nitric oxide: NOS3 activation and regulation 12 0.0164172 -0.4204784 0.1833185 0.0116768 0.2716299 0.4587023 0.4269489 0.2247639
375 Gap junction trafficking and regulation 12 0.0327984 -0.2652200 -0.3711816 0.1117221 0.0260106 0.4561988 0.0749262 0.2578836
861 Regulation of IFNA/IFNB signaling 12 0.0334829 -0.3092933 -0.3352426 0.0636194 0.0443801 0.4561249 0.0183489 0.2585049
374 Gap junction trafficking 11 0.0443257 -0.2149567 -0.3978455 0.2171206 0.0223428 0.4522029 0.1293219 0.2932274
597 Mitotic Telophase/Cytokinesis 11 0.0466857 0.3517420 0.2838517 0.0434137 0.1031362 0.4519892 0.0480057 0.3001804
1093 TNF receptor superfamily (TNFSF) members mediating non-canonical NF-kB pathway 10 0.0518836 0.4437490 0.0672367 0.0151146 0.7128039 0.4488139 0.2662344 0.3071015
340 Formation of tubulin folding intermediates by CCT/TriC 14 0.0131745 0.4359022 -0.0828050 0.0047488 0.5917549 0.4436974 0.3667814 0.2126735
1069 Synaptic adhesion-like molecules 10 0.0699489 -0.2628166 -0.3543152 0.1501882 0.0523953 0.4411483 0.0646993 0.3547503
681 Other interleukin signaling 15 0.0266821 -0.2367143 -0.3467872 0.1125314 0.0200732 0.4198750 0.0778333 0.2418578
1019 Signaling by NODAL 12 0.0582911 -0.2385514 -0.3406946 0.1525571 0.0410342 0.4159081 0.0722261 0.3278546
719 Phosphorylation of CD3 and TCR zeta chains 10 0.0580476 0.3467656 -0.2272302 0.0576264 0.2134896 0.4145841 0.4058763 0.3278546
290 Elastic fibre formation 16 0.0127777 -0.3994462 0.1079073 0.0056773 0.4550217 0.4137648 0.3587531 0.2126735
1092 TICAM1-dependent activation of IRF3/IRF7 12 0.0469723 0.4114432 0.0142990 0.0136036 0.9316688 0.4116916 0.2808234 0.3001804
335 Formation of apoptosome 10 0.0638975 0.3650266 -0.1858284 0.0456599 0.3089859 0.4096054 0.3895133 0.3438293
898 Regulation of the apoptosome activity 10 0.0638975 0.3650266 -0.1858284 0.0456599 0.3089859 0.4096054 0.3895133 0.3438293
161 Cobalamin (Cbl, vitamin B12) transport and metabolism 11 0.0797357 -0.1395866 -0.3782807 0.4228697 0.0298463 0.4032129 0.1687822 0.3753221
412 HDR through MMEJ (alt-NHEJ) 11 0.0878474 0.2842990 0.2859123 0.1025948 0.1006587 0.4032019 0.0011408 0.3926979
495 Interleukin-7 signaling 15 0.0389323 -0.2406089 -0.3171629 0.1067381 0.0334759 0.3981017 0.0541318 0.2829994
745 Prefoldin mediated transfer of substrate to CCT/TriC 18 0.0096476 0.3579903 -0.1720743 0.0085656 0.2064085 0.3971985 0.3748123 0.1898616
1073 Synthesis of IP3 and IP4 in the cytosol 16 0.0298942 -0.3617757 -0.1609204 0.0122469 0.2652332 0.3959508 0.1420262 0.2493860
897 Regulation of signaling by CBL 14 0.0320092 0.0915535 -0.3831808 0.5532112 0.0130647 0.3939665 0.3356879 0.2550479
213 Defective HDR through Homologous Recombination (HRR) due to BRCA1 loss-of-function 22 0.0052857 0.3916445 -0.0354065 0.0014763 0.7738320 0.3932417 0.3019707 0.1712989
214 Defective HDR through Homologous Recombination (HRR) due to PALB2 loss of function 22 0.0052857 0.3916445 -0.0354065 0.0014763 0.7738320 0.3932417 0.3019707 0.1712989
215 Defective HDR through Homologous Recombination Repair (HRR) due to PALB2 loss of BRCA1 binding function 22 0.0052857 0.3916445 -0.0354065 0.0014763 0.7738320 0.3932417 0.3019707 0.1712989
head(res$enrichment_result[order(res$enrichment_result$pMANOVA),],30) %>%
  kbl(caption="Top pathways found with GMEA/fry+mitch (significance)") %>%
  kable_paper("hover", full_width = F)
Top pathways found with GMEA/fry+mitch (significance)
set setSize pMANOVA s.bl s.gu p.bl p.gu s.dist SD p.adjustMANOVA
554 Metabolism of RNA 576 0.0000049 0.1217234 0.0153399 0.0000008 0.5335176 0.1226861 0.0752245 0.0061049
126 Cell Cycle 549 0.0001037 0.1009704 0.0482310 0.0000615 0.0556761 0.1118984 0.0372924 0.0644588
360 G2/M Checkpoints 122 0.0003345 0.1934542 0.1012741 0.0002301 0.0538641 0.2183598 0.0651812 0.1085590
414 HIV Infection 195 0.0004380 0.1436326 0.0930916 0.0005665 0.0254918 0.1711618 0.0357378 0.1085590
203 DNA Repair 254 0.0004668 0.1334974 0.0654974 0.0002644 0.0735534 0.1486992 0.0480833 0.1085590
1063 Stabilization of p53 50 0.0005751 0.3128257 0.0761644 0.0001311 0.3519161 0.3219642 0.1673448 0.1085590
358 G1/S Transition 117 0.0006940 0.1698568 0.1306021 0.0015370 0.0148770 0.2142621 0.0277572 0.1085590
128 Cell Cycle, Mitotic 439 0.0006987 0.0979748 0.0526103 0.0004782 0.0607715 0.1112066 0.0320776 0.1085590
307 Extracellular matrix organization 93 0.0008468 -0.2241475 0.0058509 0.0001910 0.9224450 0.2242238 0.1626334 0.1169543
925 SARS-CoV-1 targets host intracellular signalling and regulatory pathways 15 0.0009625 -0.4728207 -0.3389279 0.0015231 0.0230699 0.5817487 0.0946765 0.1196350
127 Cell Cycle Checkpoints 233 0.0011147 0.1313005 0.0639092 0.0005835 0.0941917 0.1460281 0.0476529 0.1259561
920 S Phase 146 0.0013803 0.1689349 0.0603479 0.0004393 0.2093068 0.1793903 0.0767826 0.1429752
212 Defective CFTR causes cystic fibrosis 53 0.0015533 0.2822645 0.0732634 0.0003820 0.3566420 0.2916175 0.1477860 0.1485187
517 Late Phase of HIV Life Cycle 119 0.0016786 0.1419593 0.1402224 0.0075962 0.0083681 0.1995364 0.0012282 0.1490379
908 Respiratory electron transport 85 0.0020601 0.0704994 0.2155224 0.2619223 0.0006028 0.2267599 0.1025467 0.1707106
205 DNA Replication Pre-Initiation 92 0.0025096 0.1923876 0.1009600 0.0014490 0.0947211 0.2172692 0.0646490 0.1712989
73 Autodegradation of the E3 ubiquitin ligase COP1 45 0.0026323 0.2950075 0.0654534 0.0006217 0.4478129 0.3021814 0.1623193 0.1712989
431 Hh mutants are degraded by ERAD 48 0.0033385 0.2728914 0.0982639 0.0010801 0.2392670 0.2900440 0.1234803 0.1712989
530 M Phase 309 0.0037508 0.1017881 0.0548100 0.0022089 0.0993882 0.1156069 0.0332186 0.1712989
591 Mitotic G1 phase and G1/S transition 131 0.0037717 0.1429260 0.1050124 0.0048231 0.0383793 0.1773568 0.0268090 0.1712989
646 Neuronal System 142 0.0041875 -0.1457039 -0.0835494 0.0027868 0.0864357 0.1679587 0.0439499 0.1712989
679 Orc1 removal from chromatin 63 0.0042882 0.2152378 0.1292895 0.0031594 0.0762473 0.2510839 0.0607746 0.1712989
430 Hh mutants abrogate ligand secretion 50 0.0042932 0.2660770 0.0733128 0.0011435 0.3702389 0.2759923 0.1363049 0.1712989
909 Respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins. 87 0.0045071 0.0585021 0.2005342 0.3463204 0.0012433 0.2088934 0.1004318 0.1712989
593 Mitotic Metaphase and Anaphase 203 0.0046313 0.1246772 0.0614087 0.0022746 0.1328553 0.1389799 0.0447376 0.1712989
855 Regulation of Apoptosis 46 0.0046609 0.2786414 0.0491636 0.0010844 0.5643263 0.2829454 0.1622653 0.1712989
1192 Ubiquitin Mediated Degradation of Phosphorylated Cdc25A 45 0.0048035 0.2802916 0.0566044 0.0011492 0.5115396 0.2859500 0.1581707 0.1712989
1230 p53-Independent DNA Damage Response 45 0.0048035 0.2802916 0.0566044 0.0011492 0.5115396 0.2859500 0.1581707 0.1712989
1231 p53-Independent G1/S DNA damage checkpoint 45 0.0048035 0.2802916 0.0566044 0.0011492 0.5115396 0.2859500 0.1581707 0.1712989
415 HIV Life Cycle 131 0.0050268 0.1084247 0.1347483 0.0325103 0.0078768 0.1729538 0.0186136 0.1712989
sig <- subset(res$enrichment_result,p.adjustMANOVA<0.01)

head(sig[order(-abs(sig$s.dist)),],20) %>%
  kbl(caption="Top pathways found with GMEA/fry+mitch (effect size after 1% FDR filtering)") %>%
  kable_paper("hover", full_width = F)
Top pathways found with GMEA/fry+mitch (effect size after 1% FDR filtering)
set setSize pMANOVA s.bl s.gu p.bl p.gu s.dist SD p.adjustMANOVA
554 Metabolism of RNA 576 4.9e-06 0.1217234 0.0153399 8e-07 0.5335176 0.1226861 0.0752245 0.0061049

Session information

For reproducibility

sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
##  [1] grid      stats4    parallel  stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] tictoc_1.2                                         
##  [2] kableExtra_1.3.4                                   
##  [3] data.table_1.14.8                                  
##  [4] ENmix_1.36.03                                      
##  [5] doParallel_1.0.17                                  
##  [6] qqman_0.1.9                                        
##  [7] RCircos_1.2.2                                      
##  [8] beeswarm_0.4.0                                     
##  [9] forestplot_3.1.3                                   
## [10] abind_1.4-5                                        
## [11] checkmate_2.2.0                                    
## [12] reshape2_1.4.4                                     
## [13] gplots_3.1.3                                       
## [14] eulerr_7.0.0                                       
## [15] GEOquery_2.68.0                                    
## [16] RColorBrewer_1.1-3                                 
## [17] IlluminaHumanMethylation450kmanifest_0.4.0         
## [18] topconfects_1.16.0                                 
## [19] DMRcatedata_2.18.0                                 
## [20] ExperimentHub_2.8.1                                
## [21] AnnotationHub_3.8.0                                
## [22] BiocFileCache_2.8.0                                
## [23] dbplyr_2.3.3                                       
## [24] DMRcate_2.14.1                                     
## [25] limma_3.56.2                                       
## [26] missMethyl_1.34.0                                  
## [27] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 
## [28] R.utils_2.12.2                                     
## [29] R.oo_1.25.0                                        
## [30] R.methodsS3_1.8.2                                  
## [31] plyr_1.8.8                                         
## [32] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [33] minfi_1.46.0                                       
## [34] bumphunter_1.42.0                                  
## [35] locfit_1.5-9.8                                     
## [36] iterators_1.0.14                                   
## [37] foreach_1.5.2                                      
## [38] Biostrings_2.68.1                                  
## [39] XVector_0.40.0                                     
## [40] SummarizedExperiment_1.30.2                        
## [41] Biobase_2.60.0                                     
## [42] MatrixGenerics_1.12.3                              
## [43] matrixStats_1.0.0                                  
## [44] GenomicRanges_1.52.0                               
## [45] GenomeInfoDb_1.36.3                                
## [46] IRanges_2.34.1                                     
## [47] S4Vectors_0.38.1                                   
## [48] BiocGenerics_0.46.0                                
## [49] mitch_1.12.0                                       
## 
## loaded via a namespace (and not attached):
##   [1] DSS_2.48.0                    ProtGenerics_1.32.0          
##   [3] bitops_1.0-7                  httr_1.4.7                   
##   [5] webshot_0.5.5                 dynamicTreeCut_1.63-1        
##   [7] tools_4.3.2                   doRNG_1.8.6                  
##   [9] backports_1.4.1               utf8_1.2.3                   
##  [11] R6_2.5.1                      HDF5Array_1.28.1             
##  [13] lazyeval_0.2.2                Gviz_1.44.1                  
##  [15] rhdf5filters_1.12.1           permute_0.9-7                
##  [17] prettyunits_1.1.1             GGally_2.1.2                 
##  [19] gridExtra_2.3                 base64_2.0.1                 
##  [21] preprocessCore_1.62.1         cli_3.6.1                    
##  [23] sass_0.4.7                    readr_2.1.4                  
##  [25] genefilter_1.82.1             askpass_1.2.0                
##  [27] Rsamtools_2.16.0              systemfonts_1.0.4            
##  [29] foreign_0.8-85                siggenes_1.74.0              
##  [31] illuminaio_0.42.0             svglite_2.1.1                
##  [33] dichromat_2.0-0.1             scrime_1.3.5                 
##  [35] BSgenome_1.68.0               readxl_1.4.3                 
##  [37] impute_1.74.1                 rstudioapi_0.15.0            
##  [39] RSQLite_2.3.1                 generics_0.1.3               
##  [41] BiocIO_1.10.0                 gtools_3.9.4                 
##  [43] dplyr_1.1.3                   Matrix_1.6-1.1               
##  [45] interp_1.1-4                  fansi_1.0.4                  
##  [47] lifecycle_1.0.3               edgeR_3.42.4                 
##  [49] yaml_2.3.7                    rhdf5_2.44.0                 
##  [51] blob_1.2.4                    promises_1.2.1               
##  [53] crayon_1.5.2                  lattice_0.22-5               
##  [55] echarts4r_0.4.5               GenomicFeatures_1.52.2       
##  [57] annotate_1.78.0               KEGGREST_1.40.0              
##  [59] pillar_1.9.0                  knitr_1.44                   
##  [61] beanplot_1.3.1                rjson_0.2.21                 
##  [63] codetools_0.2-19              glue_1.6.2                   
##  [65] vctrs_0.6.3                   png_0.1-8                    
##  [67] cellranger_1.1.0              gtable_0.3.4                 
##  [69] assertthat_0.2.1              cachem_1.0.8                 
##  [71] xfun_0.40                     S4Arrays_1.0.6               
##  [73] mime_0.12                     survival_3.5-7               
##  [75] statmod_1.5.0                 interactiveDisplayBase_1.38.0
##  [77] ellipsis_0.3.2                nlme_3.1-163                 
##  [79] bit64_4.0.5                   bsseq_1.36.0                 
##  [81] progress_1.2.2                filelock_1.0.2               
##  [83] bslib_0.5.1                   nor1mix_1.3-0                
##  [85] KernSmooth_2.23-22            rpart_4.1.21                 
##  [87] colorspace_2.1-0              DBI_1.1.3                    
##  [89] Hmisc_5.1-1                   nnet_7.3-19                  
##  [91] tidyselect_1.2.0              bit_4.0.5                    
##  [93] compiler_4.3.2                curl_5.0.2                   
##  [95] rvest_1.0.3                   htmlTable_2.4.1              
##  [97] xml2_1.3.5                    RPMM_1.25                    
##  [99] DelayedArray_0.26.7           rtracklayer_1.60.1           
## [101] scales_1.2.1                  caTools_1.18.2               
## [103] quadprog_1.5-8                rappdirs_0.3.3               
## [105] stringr_1.5.0                 digest_0.6.33                
## [107] rmarkdown_2.24                htmltools_0.5.6              
## [109] pkgconfig_2.0.3               jpeg_0.1-10                  
## [111] base64enc_0.1-3               sparseMatrixStats_1.12.2     
## [113] highr_0.10                    fastmap_1.1.1                
## [115] ensembldb_2.24.0              rlang_1.1.1                  
## [117] htmlwidgets_1.6.2             shiny_1.7.5                  
## [119] DelayedMatrixStats_1.22.6     jquerylib_0.1.4              
## [121] jsonlite_1.8.7                BiocParallel_1.34.2          
## [123] mclust_6.0.0                  VariantAnnotation_1.46.0     
## [125] RCurl_1.98-1.12               magrittr_2.0.3               
## [127] Formula_1.2-5                 GenomeInfoDbData_1.2.10      
## [129] Rhdf5lib_1.22.1               munsell_0.5.0                
## [131] Rcpp_1.0.11                   stringi_1.7.12               
## [133] zlibbioc_1.46.0               MASS_7.3-60                  
## [135] org.Hs.eg.db_3.17.0           deldir_1.0-9                 
## [137] splines_4.3.2                 multtest_2.56.0              
## [139] hms_1.1.3                     rngtools_1.5.2               
## [141] geneplotter_1.78.0            biomaRt_2.56.1               
## [143] BiocVersion_3.17.1            XML_3.99-0.14                
## [145] evaluate_0.21                 calibrate_1.7.7              
## [147] latticeExtra_0.6-30           biovizBase_1.48.0            
## [149] BiocManager_1.30.22           tzdb_0.4.0                   
## [151] httpuv_1.6.11                 tidyr_1.3.0                  
## [153] openssl_2.1.0                 purrr_1.0.2                  
## [155] reshape_0.8.9                 ggplot2_3.4.3                
## [157] xtable_1.8-4                  restfulr_0.0.15              
## [159] AnnotationFilter_1.24.0       later_1.3.1                  
## [161] viridisLite_0.4.2             tibble_3.2.1                 
## [163] memoise_2.0.1                 AnnotationDbi_1.62.2         
## [165] GenomicAlignments_1.36.0      cluster_2.1.4