Autism Spectrum Disorder Methylation analysis - comparison of Guthrie card and fresh blood
Mark Ziemann
2022-10-19
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/gmea/main/meth_functions.R")
library("data.table")
library("kableExtra")
library("eulerr")
library("tictoc")
})Load data
bl_design <- readRDS(file="bl_design.rds")
bl_mvals <- readRDS(file="bl_mvals.rds")
gu_design <- readRDS(file="gu_design.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.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)## 226.39 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.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)## 206.815 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
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.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.1 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
##
## 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
##
## attached base packages:
## [1] grid stats4 parallel stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] pkgload_1.3.0
## [2] GGally_2.1.2
## [3] ggplot2_3.3.6
## [4] gtools_3.9.3
## [5] tibble_3.1.8
## [6] dplyr_1.0.10
## [7] echarts4r_0.4.4
## [8] tictoc_1.1
## [9] kableExtra_1.3.4
## [10] data.table_1.14.2
## [11] ENmix_1.32.0
## [12] doParallel_1.0.17
## [13] qqman_0.1.8
## [14] RCircos_1.2.2
## [15] beeswarm_0.4.0
## [16] forestplot_3.1.0
## [17] abind_1.4-5
## [18] checkmate_2.1.0
## [19] reshape2_1.4.4
## [20] gplots_3.1.3
## [21] eulerr_6.1.1
## [22] GEOquery_2.64.2
## [23] RColorBrewer_1.1-3
## [24] IlluminaHumanMethylation450kmanifest_0.4.0
## [25] topconfects_1.12.0
## [26] DMRcatedata_2.14.0
## [27] ExperimentHub_2.4.0
## [28] AnnotationHub_3.4.0
## [29] BiocFileCache_2.4.0
## [30] dbplyr_2.2.1
## [31] DMRcate_2.10.0
## [32] limma_3.52.4
## [33] missMethyl_1.30.0
## [34] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [35] R.utils_2.12.0
## [36] R.oo_1.25.0
## [37] R.methodsS3_1.8.2
## [38] plyr_1.8.7
## [39] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [40] minfi_1.42.0
## [41] bumphunter_1.38.0
## [42] locfit_1.5-9.6
## [43] iterators_1.0.14
## [44] foreach_1.5.2
## [45] Biostrings_2.64.1
## [46] XVector_0.36.0
## [47] SummarizedExperiment_1.26.1
## [48] Biobase_2.56.0
## [49] MatrixGenerics_1.8.1
## [50] matrixStats_0.62.0
## [51] GenomicRanges_1.48.0
## [52] GenomeInfoDb_1.32.4
## [53] IRanges_2.30.1
## [54] S4Vectors_0.34.0
## [55] BiocGenerics_0.42.0
## [56] mitch_1.8.0
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 rtracklayer_1.56.1
## [3] tidyr_1.2.1 bit64_4.0.5
## [5] knitr_1.40 DelayedArray_0.22.0
## [7] rpart_4.1.16 KEGGREST_1.36.3
## [9] RCurl_1.98-1.9 AnnotationFilter_1.20.0
## [11] generics_0.1.3 GenomicFeatures_1.48.4
## [13] preprocessCore_1.58.0 RSQLite_2.2.18
## [15] bit_4.0.4 tzdb_0.3.0
## [17] webshot_0.5.4 xml2_1.3.3
## [19] httpuv_1.6.6 assertthat_0.2.1
## [21] xfun_0.33 hms_1.1.2
## [23] jquerylib_0.1.4 evaluate_0.17
## [25] promises_1.2.0.1 fansi_1.0.3
## [27] restfulr_0.0.15 scrime_1.3.5
## [29] progress_1.2.2 caTools_1.18.2
## [31] readxl_1.4.1 DBI_1.1.3
## [33] geneplotter_1.74.0 htmlwidgets_1.5.4
## [35] reshape_0.8.9 purrr_0.3.5
## [37] ellipsis_0.3.2 backports_1.4.1
## [39] permute_0.9-7 calibrate_1.7.7
## [41] annotate_1.74.0 biomaRt_2.52.0
## [43] deldir_1.0-6 sparseMatrixStats_1.8.0
## [45] vctrs_0.4.2 ensembldb_2.20.2
## [47] withr_2.5.0 cachem_1.0.6
## [49] Gviz_1.40.1 BSgenome_1.64.0
## [51] GenomicAlignments_1.32.1 prettyunits_1.1.1
## [53] mclust_5.4.10 svglite_2.1.0
## [55] cluster_2.1.4 RPMM_1.25
## [57] lazyeval_0.2.2 crayon_1.5.2
## [59] genefilter_1.78.0 labeling_0.4.2
## [61] edgeR_3.38.4 pkgconfig_2.0.3
## [63] nlme_3.1-160 ProtGenerics_1.28.0
## [65] nnet_7.3-18 rlang_1.0.6
## [67] lifecycle_1.0.3 filelock_1.0.2
## [69] dichromat_2.0-0.1 cellranger_1.1.0
## [71] rngtools_1.5.2 base64_2.0.1
## [73] Matrix_1.5-1 Rhdf5lib_1.18.2
## [75] base64enc_0.1-3 viridisLite_0.4.1
## [77] png_0.1-7 rjson_0.2.21
## [79] bitops_1.0-7 KernSmooth_2.23-20
## [81] rhdf5filters_1.8.0 blob_1.2.3
## [83] DelayedMatrixStats_1.18.1 doRNG_1.8.2
## [85] stringr_1.4.1 nor1mix_1.3-0
## [87] readr_2.1.3 jpeg_0.1-9
## [89] scales_1.2.1 memoise_2.0.1
## [91] magrittr_2.0.3 zlibbioc_1.42.0
## [93] compiler_4.2.1 BiocIO_1.6.0
## [95] illuminaio_0.38.0 Rsamtools_2.12.0
## [97] cli_3.4.1 DSS_2.44.0
## [99] htmlTable_2.4.1 Formula_1.2-4
## [101] MASS_7.3-58.1 tidyselect_1.2.0
## [103] stringi_1.7.8 highr_0.9
## [105] yaml_2.3.5 askpass_1.1
## [107] latticeExtra_0.6-30 sass_0.4.2
## [109] VariantAnnotation_1.42.1 tools_4.2.1
## [111] rstudioapi_0.14 foreign_0.8-83
## [113] bsseq_1.32.0 gridExtra_2.3
## [115] farver_2.1.1 digest_0.6.29
## [117] BiocManager_1.30.18 shiny_1.7.2
## [119] quadprog_1.5-8 Rcpp_1.0.9
## [121] siggenes_1.70.0 BiocVersion_3.15.2
## [123] later_1.3.0 org.Hs.eg.db_3.15.0
## [125] httr_1.4.4 AnnotationDbi_1.58.0
## [127] biovizBase_1.44.0 colorspace_2.0-3
## [129] rvest_1.0.3 XML_3.99-0.11
## [131] splines_4.2.1 statmod_1.4.37
## [133] multtest_2.52.0 systemfonts_1.0.4
## [135] xtable_1.8-4 jsonlite_1.8.2
## [137] dynamicTreeCut_1.63-1 R6_2.5.1
## [139] Hmisc_4.7-1 pillar_1.8.1
## [141] htmltools_0.5.3 mime_0.12
## [143] glue_1.6.2 fastmap_1.1.0
## [145] BiocParallel_1.30.4 interactiveDisplayBase_1.34.0
## [147] beanplot_1.3.1 codetools_0.2-18
## [149] utf8_1.2.2 lattice_0.20-45
## [151] bslib_0.4.0 curl_4.3.3
## [153] openssl_2.0.3 interp_1.1-3
## [155] survival_3.4-0 rmarkdown_2.17
## [157] munsell_0.5.0 rhdf5_2.40.0
## [159] GenomeInfoDbData_1.2.8 HDF5Array_1.24.2
## [161] impute_1.70.0 gtable_0.3.1