Autism Spectrum Disorder Methylation analysis - comparison of Guthrie card and fresh blood
Mark Ziemann
2023-11-08
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.00GMEA 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 elapsedhead(res_bl_df,30) %>%
kbl(caption="Top DM genes in blood identified with GMEA") %>%
kable_paper("hover", full_width = F)| 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 elapsedhead(res_gu_df,30) %>%
kbl(caption="Top DM genes in Guthrie card identified with GMEA") %>%
kable_paper("hover", full_width = F)| 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)| 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)| 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)| 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