Correlates of ADOS and methylation - enrichment analysis
Mark Ziemann
2023-03-20
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
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("RIdeogram")
library("GenomicRanges")
library("tictoc")
library("globaltest")
library("ebGSEA")
data("dualmap850kEID")
})
source("meth_functions.R")Load data
bl_mvals <- readRDS(file="bl_mvals.rds")
gu_mvals <- readRDS(file="gu_mvals.rds")
bl_design <- readRDS(file="bl_design_ados.rds")
gu_design <- readRDS(file="gu_design_ados.rds")
bl_lim <- read.csv("limma_blood_ADOS.csv")
gu_lim <- read.csv("limma_guthrie_ADOS.csv")Probe sets
For each gene, extract out the probes.
anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","Regulatory_Feature_Group","Islands_Name","Relation_to_Island")])
promoters <- grep("Prom",myann$Regulatory_Feature_Group)
promoters <- myann[promoters,]
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
sets <- split(rep(names(gp2), lengths(gp2)), unlist(gp2))
summary(unlist(lapply(sets,length)))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 9.00 24.00 49.68 55.00 6778.00Gene level enrichment 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
tic()
res_bl <- main(mval=bl_mvals,design=bl_design,sets,cores=detectCores()/2)
res_bl[[1]] <- subset(res_bl[[1]],nprobes!=0) # remove genes with no probes obvs
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)## 523.84 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 | |
|---|---|---|---|---|
| ANKRD29 | 31 | -0.5119264 | 0.0000190 | 0.5185575 |
| NUP43 | 20 | 0.5781906 | 0.0001144 | 0.7930777 |
| CXADRP3 | 2 | -2.2005853 | 0.0004065 | 0.7930777 |
| LOC101926941 | 18 | -0.5469809 | 0.0004103 | 0.7930777 |
| OR10W1 | 6 | -1.7758388 | 0.0007027 | 0.7930777 |
| ERO1LB | 13 | -0.7489410 | 0.0007159 | 0.7930777 |
| ATP5I | 24 | -0.5934763 | 0.0008169 | 0.7930777 |
| C12orf10 | 29 | -0.8328921 | 0.0008862 | 0.7930777 |
| LOC101929284 | 6 | -2.0125155 | 0.0009618 | 0.7930777 |
| GNL2 | 20 | -0.8487822 | 0.0009900 | 0.7930777 |
| C9orf89 | 24 | -0.7213163 | 0.0012176 | 0.7930777 |
| SLC52A2 | 6 | -0.8921214 | 0.0013174 | 0.7930777 |
| TRPC1 | 25 | -0.2836119 | 0.0014168 | 0.7930777 |
| POLR3G | 22 | -0.6011733 | 0.0014716 | 0.7930777 |
| TTC29 | 29 | -1.0539440 | 0.0015060 | 0.7930777 |
| LOC101928663 | 2 | -1.9067665 | 0.0015507 | 0.7930777 |
| CEP135 | 25 | -0.5072602 | 0.0016338 | 0.7930777 |
| LMO4 | 29 | -0.7910709 | 0.0016870 | 0.7930777 |
| LINC01220 | 2 | -1.4641912 | 0.0016962 | 0.7930777 |
| SURF1 | 17 | -0.3737872 | 0.0017982 | 0.7930777 |
| C5orf24 | 22 | -0.5231647 | 0.0018726 | 0.7930777 |
| RELN | 63 | -0.6171020 | 0.0018891 | 0.7930777 |
| ZCCHC2 | 30 | -0.5628906 | 0.0018896 | 0.7930777 |
| C3orf54 | 8 | -0.9636677 | 0.0018978 | 0.7930777 |
| CPLX4 | 13 | -1.4227283 | 0.0019430 | 0.7930777 |
| SNORA28 | 5 | -0.9776395 | 0.0019468 | 0.7930777 |
| NOL3 | 20 | -0.7747351 | 0.0022516 | 0.7930777 |
| C1orf114 | 12 | -0.7054945 | 0.0022650 | 0.7930777 |
| SNF8 | 27 | -0.5223357 | 0.0023644 | 0.7930777 |
| LOC101928823 | 11 | 0.5870690 | 0.0024001 | 0.7930777 |
par(mfrow=c(2,1))
hist(res_bl_top$t,xlab="probe t-stat",main="blood at assessment") ; grid() ; abline(v=0,lty=1,lwd=3)
hist(res_bl_df$median,xlab="gene median t-stat",main="blood at assessment") ; grid() ; abline(v=0,lty=1,lwd=3)
pdf("hist_bl_ados.pdf")
par(mfrow=c(2,1))
hist(res_bl_top$t,xlab="probe t-stat",main="blood at assessment") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
hist(res_bl_df$median,xlab="gene median t-stat",main="blood at assessment") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
dev.off()## png
## 2par(mfrow=c(1,1))Now some visualisations.
volcano_plot(res_bl_df)
probe_bias(res_bl_df)
gmea_boxplot(res_bl,sets)
Gene level enrichment for guthrie cards
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[[1]] <- subset(res_gu[[1]],nprobes!=0) # remove genes with no probes obvs
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)## 514.66 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 | |
|---|---|---|---|---|
| MIR3917 | 5 | -0.8808877 | 0.0000710 | 1 |
| ARL5A | 15 | 1.0478366 | 0.0001918 | 1 |
| ADAMTS1 | 19 | 0.6102581 | 0.0002838 | 1 |
| CHMP1B | 19 | 0.4557641 | 0.0004875 | 1 |
| FOXD1 | 11 | 0.7472067 | 0.0005569 | 1 |
| ZNF322B | 1 | -4.0851785 | 0.0006000 | 1 |
| NAA30 | 21 | -0.3732929 | 0.0006911 | 1 |
| ALPP | 14 | -0.6916978 | 0.0009992 | 1 |
| CAAP1 | 7 | 0.7857603 | 0.0010091 | 1 |
| RPS14 | 21 | 0.4654341 | 0.0010451 | 1 |
| LEPROTL1 | 19 | -0.4020884 | 0.0010934 | 1 |
| ZNF226 | 17 | 0.7738140 | 0.0016195 | 1 |
| MEOX1 | 28 | -0.8202385 | 0.0017083 | 1 |
| SFRS12IP1 | 16 | 0.8951516 | 0.0018488 | 1 |
| MIR196A2 | 9 | 0.6807623 | 0.0021186 | 1 |
| PCP4 | 17 | -0.7132310 | 0.0022340 | 1 |
| SEC14L4 | 21 | -0.8930770 | 0.0025559 | 1 |
| GZMA | 7 | 0.4791645 | 0.0026386 | 1 |
| ARL16 | 20 | 0.4894063 | 0.0027294 | 1 |
| LINC00866 | 8 | 0.9125732 | 0.0029358 | 1 |
| RECQL | 20 | 0.4446618 | 0.0029434 | 1 |
| GNG10 | 10 | -0.8812882 | 0.0036663 | 1 |
| LSM1 | 26 | 0.6178880 | 0.0036880 | 1 |
| CSN1S2B | 4 | -1.8839162 | 0.0037939 | 1 |
| LOC441046 | 11 | 0.5893277 | 0.0038516 | 1 |
| ZNF600 | 5 | 0.5349287 | 0.0038972 | 1 |
| LOC286135 | 4 | -1.5617343 | 0.0039589 | 1 |
| C1QTNF9B-AS1 | 3 | 1.7431989 | 0.0040010 | 1 |
| TMEM223 | 14 | 1.0364105 | 0.0040902 | 1 |
| LOC101928766 | 8 | -1.5651393 | 0.0043081 | 1 |
par(mfrow=c(2,1))
hist(res_gu_top$t,xlab="probe t-stat",main="neonatal Guthrie card") ; grid() ; abline(v=0,lty=1,lwd=3)
hist(res_gu_df$median,xlab="gene median t-stat",main="neonatal Guthrie card") ; grid() ; abline(v=0,lty=1,lwd=3)
pdf("hist_gu_ados.pdf")
par(mfrow=c(2,1))
hist(res_gu_top$t,xlab="probe t-stat",main="neonatal Guthrie card") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
hist(res_gu_df$median,xlab="gene median t-stat",main="neonatal Guthrie card") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
dev.off()## png
## 2par(mfrow=c(1,1))Now some visualisations.
volcano_plot(res_gu_df)
probe_bias(res_gu_df)
gmea_boxplot(res_gu,sets)
Comparison of blood and guthrie card
Now compare blood and guthrie card
res_bl_df$bl <- res_bl_df$median
res_gu_df$gu <- res_gu_df$median
m1 <- merge(res_bl_df,res_gu_df,by=0)
m2 <- m1[,c("Row.names","bl","gu")]
rownames(m2) <- m2$Row.names
m2$Row.names=NULL
plot(m2$bl,m2$gu,xlab="bl",ylab="gu",
col = rgb(red = 0.6, green = 0.6, blue = 0.6, alpha = 0.5) ,
pch=19,cex=0.9)
grid()
abline(v=0,lty=1,lwd=3,col="black") ; abline(h=0,lty=1,lwd=3,col="black")
mtext("aggregated t-statistic for each gene (median)")
#plot(rank(m1$bl),rank(m1$gu),xlab="bl rank",ylab="gu rank",cex=0.6,pch=19) ; grid()
rm2 <- apply(m2,2,rank)
mydiff <- apply(m2,2,function(x) { length(which(x<0)) } )
rm3 <- rm2
rm3[,1] <- rm3[,1] - mydiff[1]
rm3[,2] <- rm3[,2] - mydiff[2]
rnk <- rm3
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 = paste("Rank in contrast", colnames(rnk)[1])
Y_AXIS = paste("Rank in contrast", colnames(rnk)[2])
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 = "Rank-rank plot of all genes", xlab = X_AXIS,
ylab = Y_AXIS)
})
pdf("genedens_ados.pdf")
par(mfrow=c(2,2))
plot(1)
plot(1)
plot(m2$bl,m2$gu,xlab="bl",ylab="gu",
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 (median)")
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 = "Rank-rank plot of all genes", xlab = X_AXIS,
ylab = Y_AXIS)
})
dev.off()## png
## 2Now we need a table of genes with top scores.
m3 <- m2
m3$med <- apply(rnk,1,median)
mg <- merge(res_gu_df,res_bl_df,by=0)
mg$gu = mg$bl = NULL
mg$med <- m3$med
mg2 <- data.frame(mg,rnk)
mg2$med <- apply(rnk,1,median)
mg2 <- mg2[order(mg2$med),]
mg2 <-subset(mg2,nprobes.x>=5)
head(mg2,20) %>%
kbl(caption="Genes with coordinated hypomethylation in neonatal Guthrie card and blood at assessment") %>%
kable_paper("hover", full_width = F)| Row.names | nprobes.x | median.x | PValue.x | FDR.x | nprobes.y | median.y | PValue.y | FDR.y | med | bl | gu | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| MIR1251 | MIR1251 | 6 | -1.4274015 | 0.1774551 | 1 | 6 | -1.792361 | 0.1506930 | 0.7930777 | -19529.5 | -21220.0 | -17839.0 |
| SNORD114-16 | SNORD114-16 | 5 | -1.1449573 | 0.5203532 | 1 | 5 | -1.957956 | 0.0827893 | 0.7930777 | -19468.5 | -21307.0 | -17630.0 |
| LOC101928273 | LOC101928273 | 5 | -1.2127118 | 0.3877187 | 1 | 5 | -1.669061 | 0.0145385 | 0.7930777 | -19415.5 | -21132.0 | -17699.0 |
| LINC00375 | LINC00375 | 5 | -1.3832814 | 0.2943887 | 1 | 6 | -1.512424 | 0.2761556 | 0.7930777 | -19355.0 | -20893.0 | -17817.0 |
| MIR568 | MIR568 | 6 | -0.9928982 | 0.4417367 | 1 | 6 | -1.914089 | 0.0818591 | 0.7930777 | -19332.5 | -21290.0 | -17375.0 |
| BZW1L1 | BZW1L1 | 5 | -1.1110601 | 0.2114301 | 1 | 5 | -1.596645 | 0.1350618 | 0.7930777 | -19312.0 | -21044.0 | -17580.0 |
| XIRP2-AS1 | XIRP2-AS1 | 8 | -1.2066650 | 0.3419524 | 1 | 8 | -1.502519 | 0.0140950 | 0.7930777 | -19283.5 | -20874.0 | -17693.0 |
| LOC400655 | LOC400655 | 9 | -0.9293551 | 0.3457388 | 1 | 9 | -2.085337 | 0.0094342 | 0.7930777 | -19278.0 | -21341.0 | -17215.0 |
| HISPPD1 | HISPPD1 | 5 | -1.1073882 | 0.4421820 | 1 | 5 | -1.551529 | 0.1827707 | 0.7930777 | -19268.0 | -20964.0 | -17572.0 |
| PTH | PTH | 5 | -0.9789510 | 0.4148407 | 1 | 6 | -1.698846 | 0.0563605 | 0.7930777 | -19255.5 | -21159.0 | -17352.0 |
| TAAR9 | TAAR9 | 5 | -1.1557224 | 0.2896427 | 1 | 6 | -1.454520 | 0.1644642 | 0.7930777 | -19211.0 | -20778.0 | -17644.0 |
| OR6T1 | OR6T1 | 7 | -1.0586676 | 0.3326189 | 1 | 8 | -1.510381 | 0.1210877 | 0.7930777 | -19199.0 | -20889.0 | -17509.0 |
| LOC101927356 | LOC101927356 | 5 | -0.9179596 | 0.5002259 | 1 | 6 | -1.785762 | 0.0042609 | 0.7930777 | -19193.5 | -21215.0 | -17172.0 |
| SLC13A1 | SLC13A1 | 8 | -1.0568968 | 0.4250764 | 1 | 9 | -1.477957 | 0.1210331 | 0.7930777 | -19168.0 | -20834.0 | -17502.0 |
| TMPRSS11F | TMPRSS11F | 11 | -0.8784111 | 0.3799398 | 1 | 12 | -1.767057 | 0.0298130 | 0.7930777 | -19133.5 | -21205.0 | -17062.0 |
| DDTL | DDTL | 5 | -1.2451143 | 0.2211730 | 1 | 5 | -1.337351 | 0.1534973 | 0.7930777 | -19107.5 | -20488.0 | -17727.0 |
| CRYGC | CRYGC | 9 | -1.0422063 | 0.1633172 | 1 | 9 | -1.435018 | 0.0982049 | 0.7930777 | -19106.0 | -20736.0 | -17476.0 |
| SNORD114-26 | SNORD114-26 | 5 | -0.8547406 | 0.4835089 | 1 | 5 | -1.817857 | 0.0957504 | 0.7930777 | -19104.5 | -21237.5 | -16971.5 |
| ADGRL4 | ADGRL4 | 5 | -1.2355909 | 0.1427681 | 1 | 5 | -1.336956 | 0.3894606 | 0.8010867 | -19103.0 | -20484.0 | -17722.0 |
| RMST | RMST | 18 | -0.9256270 | 0.2833991 | 1 | 18 | -1.569118 | 0.0336580 | 0.7930777 | -19101.5 | -21000.0 | -17203.0 |
mg2 <- mg2[order(-mg2$med),]
head(mg2,20) %>%
kbl(caption="Genes with coordinated hypermethylation in neonatal Guthrie card and blood at assessment") %>%
kable_paper("hover", full_width = F)| Row.names | nprobes.x | median.x | PValue.x | FDR.x | nprobes.y | median.y | PValue.y | FDR.y | med | bl | gu | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| GPX8 | GPX8 | 9 | 0.7312060 | 0.2336383 | 1 | 9 | 1.4334798 | 0.0458087 | 0.7930777 | 5749.0 | 4227.0 | 7271.0 |
| LOC400027 | LOC400027 | 10 | 0.8814881 | 0.1105989 | 1 | 10 | 0.9710774 | 0.1825455 | 0.7930777 | 5741.5 | 4053.0 | 7430.0 |
| CARD6 | CARD6 | 5 | 1.0981604 | 0.0363444 | 1 | 7 | 0.7900422 | 0.4244594 | 0.8070399 | 5714.5 | 3887.0 | 7542.0 |
| LOC728819 | LOC728819 | 9 | 1.3835518 | 0.0238339 | 1 | 9 | 0.6906842 | 0.2130265 | 0.7930777 | 5712.5 | 3786.0 | 7639.0 |
| HIST1H3E | HIST1H3E | 15 | 0.6501439 | 0.4228645 | 1 | 15 | 1.1093279 | 0.0751159 | 0.7930777 | 5625.0 | 4115.0 | 7135.0 |
| FLJ35776 | FLJ35776 | 9 | 0.7476523 | 0.3584935 | 1 | 9 | 0.8601087 | 0.1781379 | 0.7930777 | 5621.0 | 3953.0 | 7289.0 |
| DNAJB4 | DNAJB4 | 8 | 0.6526773 | 0.0212003 | 1 | 9 | 1.0688812 | 0.2063266 | 0.7930777 | 5618.5 | 4096.0 | 7141.0 |
| ZSCAN16 | ZSCAN16 | 8 | 0.7482512 | 0.1566899 | 1 | 8 | 0.8110517 | 0.1540819 | 0.7930777 | 5601.5 | 3912.0 | 7291.0 |
| CAAP1 | CAAP1 | 7 | 0.7857603 | 0.0010091 | 1 | 8 | 0.7258485 | 0.2840578 | 0.7930777 | 5584.0 | 3825.0 | 7343.0 |
| SMEK2 | SMEK2 | 14 | 0.6589292 | 0.4376319 | 1 | 14 | 0.9050071 | 0.1995825 | 0.7930777 | 5573.5 | 3993.0 | 7154.0 |
| HIST1H2BD | HIST1H2BD | 19 | 0.7566297 | 0.1791553 | 1 | 20 | 0.7336926 | 0.2696345 | 0.7930777 | 5572.0 | 3841.0 | 7303.0 |
| HIST1H2AD | HIST1H2AD | 8 | 0.8016445 | 0.1408507 | 1 | 8 | 0.6835344 | 0.1461026 | 0.7930777 | 5569.5 | 3777.5 | 7361.5 |
| HIST1H2BF | HIST1H2BF | 8 | 0.8016445 | 0.1408507 | 1 | 8 | 0.6835344 | 0.1461026 | 0.7930777 | 5569.5 | 3777.5 | 7361.5 |
| HIST1H2AC | HIST1H2AC | 16 | 0.7649495 | 0.0425652 | 1 | 16 | 0.6878030 | 0.4154244 | 0.8047133 | 5549.0 | 3782.0 | 7316.0 |
| LOC100128822 | LOC100128822 | 9 | 0.6226731 | 0.2257377 | 1 | 9 | 0.9326637 | 0.2320806 | 0.7930777 | 5541.5 | 4016.0 | 7067.0 |
| RBBP5 | RBBP5 | 11 | 0.7479680 | 0.1374463 | 1 | 11 | 0.6800551 | 0.1405112 | 0.7930777 | 5531.5 | 3773.0 | 7290.0 |
| SIRT4 | SIRT4 | 10 | 0.7740743 | 0.1589140 | 1 | 13 | 0.6493541 | 0.3554988 | 0.7937780 | 5528.5 | 3730.0 | 7327.0 |
| ARMC1 | ARMC1 | 15 | 0.6852645 | 0.0916985 | 1 | 16 | 0.6967156 | 0.4803700 | 0.8219941 | 5496.0 | 3792.0 | 7200.0 |
| ZEB2-AS1 | ZEB2-AS1 | 11 | 0.5886547 | 0.2622867 | 1 | 11 | 0.8717299 | 0.2232118 | 0.7930777 | 5472.5 | 3966.0 | 6979.0 |
| PRR15L | PRR15L | 14 | 0.6222995 | 0.1417170 | 1 | 14 | 0.7687025 | 0.4987191 | 0.8275530 | 5469.5 | 3874.0 | 7065.0 |
Gene set level enrichment
genesets <- gmt_import("ReactomePathways.gmt")
cores = 8Mitch enrichment analysis
Here I’m using the median t-statistic for downstream enrichment. Pathways are from REACTOME and analysis using mitch.
res_bl_df$bl <- res_bl_df$median
res_gu_df$gu <- res_gu_df$median
bl <- res_bl_df[,"bl",drop=FALSE]
mitch_bl <- mitch_calc(bl,genesets,priority="effect")
head( mitch_bl$enrichment_result,30) %>%
kbl(caption="BLOOD: Top effect size pathways found with mitch") %>%
kable_paper("hover", full_width = F)
sig <- subset(mitch_bl$enrichment_result,`p.adjustANOVA`<0.05)
head(sig[order(-abs(sig$s.dist)),],30) %>%
kbl(caption="BLOOD: Top effect size pathways found with mitch after 1% FDR filtering") %>%
kable_paper("hover", full_width = F)
gu <- res_gu_df[,"gu",drop=FALSE]
mitch_gu <- mitch_calc(gu,genesets,priority="effect")
head( mitch_gu$enrichment_result,30) %>%
kbl(caption="GUTHRIE: Top effect size pathways found with mitch") %>%
kable_paper("hover", full_width = F)
sig <- subset(mitch_gu$enrichment_result,`p.adjustANOVA`<0.05)
head(sig[order(-abs(sig$s.dist)),],30) %>%
kbl(caption="GUTHRIE: Top effect size pathways found with mitch after 1% FDR filtering") %>%
kable_paper("hover", full_width = F)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(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)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpZzM26I/blgu_mitch.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|. | 3%
|
|.. | 6% (checklibraries)
|
|.... | 9%
|
|..... | 12% (peek)
|
|...... | 15%
|
|....... | 18% (metrics)
|
|........ | 21%
|
|.......... | 24% (scatterplot)
|
|........... | 26%
|
|............ | 29% (contourplot)
|
|............. | 32%
|
|.............. | 35% (input_geneset_metrics1)
|
|................ | 38%
|
|................. | 41% (input_geneset_metrics2)
|
|.................. | 44%
|
|................... | 47% (input_geneset_metrics3)
|
|.................... | 50%
|
|...................... | 53% (echart1d)
|
|....................... | 56% (echart2d)
|
|........................ | 59%
|
|......................... | 62% (heatmap)
|
|........................... | 65%
|
|............................ | 68% (effectsize)
|
|............................. | 71%
|
|.............................. | 74% (results_table)
|
|............................... | 76%
|
|................................. | 79% (results_table_complete)
|
|.................................. | 82%
|
|................................... | 85% (detailed_geneset_reports1d)
|
|.................................... | 88%
|
|..................................... | 91% (detailed_geneset_reports2d)
|
|....................................... | 94%
|
|........................................ | 97% (session_info)
|
|.........................................| 100% ## output file: /home/mdz/projects/asd_meth/mitch.knit.md## /usr/bin/pandoc +RTS -K512m -RTS /home/mdz/projects/asd_meth/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpZzM26I/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /tmp/RtmpZzM26I/rmarkdown-str3a502368b283af.html##
## Output created: /tmp/RtmpZzM26I/mitch_report.html## [1] TRUEmitch_plots(res=res,outfile="blgu_mitch_ados.pdf")## png
## 2sig <- subset(res$enrichment_result,p.adjustMANOVA<0.01)
head(sig[order(-abs(sig$s.dist)),],30) %>%
kbl(caption="Top pathways found with mitch (effect size after 1% FDR filtering)") %>%
kable_paper("hover", full_width = F)| set | setSize | pMANOVA | s.bl | s.gu | p.bl | p.gu | s.dist | SD | p.adjustMANOVA | |
|---|---|---|---|---|---|---|---|---|---|---|
| 212 | Class C/3 (Metabotropic glutamate/pheromone receptors) | 39 | 0.0000000 | -0.5260969 | -0.4451724 | 0.0000000 | 0.0000015 | 0.6891708 | 0.0572222 | 0.0000000 |
| 380 | Eicosanoids | 12 | 0.0000841 | 0.6502295 | -0.1591905 | 0.0000959 | 0.3396372 | 0.6694327 | 0.5723464 | 0.0007698 |
| 516 | Glucuronidation | 25 | 0.0000081 | -0.4846569 | -0.3813647 | 0.0000272 | 0.0009624 | 0.6167101 | 0.0730386 | 0.0000914 |
| 1163 | SARS-CoV-1 modulates host translation machinery | 34 | 0.0000001 | 0.5291174 | 0.3014517 | 0.0000001 | 0.0023452 | 0.6089650 | 0.1609840 | 0.0000022 |
| 1215 | Sensory perception of sweet, bitter, and umami (glutamate) taste | 41 | 0.0000002 | -0.4475161 | -0.3285582 | 0.0000007 | 0.0002713 | 0.5551767 | 0.0841159 | 0.0000030 |
| 903 | Peptide chain elongation | 84 | 0.0000000 | 0.4680166 | 0.2845095 | 0.0000000 | 0.0000065 | 0.5477091 | 0.1297591 | 0.0000000 |
| 395 | Eukaryotic Translation Elongation | 88 | 0.0000000 | 0.4639213 | 0.2776535 | 0.0000000 | 0.0000067 | 0.5406611 | 0.1317112 | 0.0000000 |
| 95 | Aspirin ADME | 44 | 0.0000002 | -0.3723122 | -0.3879403 | 0.0000192 | 0.0000084 | 0.5376933 | 0.0110507 | 0.0000033 |
| 397 | Eukaryotic Translation Termination | 87 | 0.0000000 | 0.4561885 | 0.2622115 | 0.0000000 | 0.0000234 | 0.5261776 | 0.1371625 | 0.0000000 |
| 1208 | Selenocysteine synthesis | 87 | 0.0000000 | 0.4554126 | 0.2617756 | 0.0000000 | 0.0000242 | 0.5252876 | 0.1369221 | 0.0000000 |
| 1496 | Viral mRNA Translation | 84 | 0.0000000 | 0.4449038 | 0.2721301 | 0.0000000 | 0.0000161 | 0.5215306 | 0.1221695 | 0.0000000 |
| 882 | PINK1-PRKN Mediated Mitophagy | 21 | 0.0010561 | 0.3476404 | 0.3816325 | 0.0058111 | 0.0024618 | 0.5162337 | 0.0240361 | 0.0077998 |
| 1172 | SCF(Skp2)-mediated degradation of p27/p21 | 58 | 0.0000000 | 0.3548917 | 0.3678922 | 0.0000029 | 0.0000012 | 0.5111681 | 0.0091928 | 0.0000003 |
| 1145 | Response of EIF2AK4 (GCN2) to amino acid deficiency | 95 | 0.0000000 | 0.4349631 | 0.2537865 | 0.0000000 | 0.0000189 | 0.5035876 | 0.1281112 | 0.0000000 |
| 829 | Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) | 89 | 0.0000000 | 0.4261122 | 0.2632858 | 0.0000000 | 0.0000175 | 0.5008902 | 0.1151357 | 0.0000000 |
| 106 | Autodegradation of Cdh1 by Cdh1:APC/C | 62 | 0.0000000 | 0.3209240 | 0.3683680 | 0.0000123 | 0.0000005 | 0.4885562 | 0.0335480 | 0.0000003 |
| 449 | Formation of a pool of free 40S subunits | 94 | 0.0000000 | 0.4283806 | 0.2339271 | 0.0000000 | 0.0000880 | 0.4880900 | 0.1374994 | 0.0000000 |
| 399 | Expression and translocation of olfactory receptors | 356 | 0.0000000 | -0.3013619 | -0.3772402 | 0.0000000 | 0.0000000 | 0.4828345 | 0.0536541 | 0.0000000 |
| 1216 | Sensory perception of taste | 47 | 0.0000015 | -0.3564503 | -0.3256706 | 0.0000234 | 0.0001116 | 0.4828231 | 0.0217646 | 0.0000220 |
| 224 | Competing endogenous RNAs (ceRNAs) regulate PTEN translation | 16 | 0.0012294 | -0.4560554 | 0.1573322 | 0.0015841 | 0.2758674 | 0.4824313 | 0.4337306 | 0.0086791 |
| 753 | Mitophagy | 27 | 0.0004874 | 0.3124424 | 0.3641463 | 0.0049452 | 0.0010537 | 0.4798154 | 0.0365602 | 0.0038445 |
| 1213 | Senescence-Associated Secretory Phenotype (SASP) | 53 | 0.0000003 | 0.2653677 | 0.3980270 | 0.0008290 | 0.0000005 | 0.4783780 | 0.0938042 | 0.0000047 |
| 13 | APC/C:Cdc20 mediated degradation of Cyclin B | 24 | 0.0011445 | 0.2857018 | 0.3826321 | 0.0153866 | 0.0011729 | 0.4775278 | 0.0685401 | 0.0082466 |
| 853 | Olfactory Signaling Pathway | 363 | 0.0000000 | -0.2929109 | -0.3740650 | 0.0000000 | 0.0000000 | 0.4751015 | 0.0573846 | 0.0000000 |
| 1169 | SARS-CoV-2 modulates host translation machinery | 46 | 0.0000023 | 0.4002583 | 0.2532630 | 0.0000026 | 0.0029527 | 0.4736548 | 0.1039414 | 0.0000312 |
| 492 | GSK3B and BTRC:CUL1-mediated-degradation of NFE2L2 | 50 | 0.0000013 | 0.3210559 | 0.3452131 | 0.0000854 | 0.0000240 | 0.4714329 | 0.0170818 | 0.0000186 |
| 1404 | The role of GTSE1 in G2/M progression after G2 checkpoint | 58 | 0.0000001 | 0.3562880 | 0.3075098 | 0.0000027 | 0.0000507 | 0.4706415 | 0.0344914 | 0.0000028 |
| 493 | GTP hydrolysis and joining of the 60S ribosomal subunit | 104 | 0.0000000 | 0.4097452 | 0.2233373 | 0.0000000 | 0.0000824 | 0.4666591 | 0.1318103 | 0.0000000 |
| 565 | Hh mutants are degraded by ERAD | 54 | 0.0000006 | 0.3489444 | 0.3087141 | 0.0000091 | 0.0000865 | 0.4659041 | 0.0284471 | 0.0000093 |
| 21 | AUF1 (hnRNP D0) binds and destabilizes mRNA | 53 | 0.0000008 | 0.3526228 | 0.3039961 | 0.0000089 | 0.0001283 | 0.4655711 | 0.0343842 | 0.0000117 |
Camera enrichment analysis
Camera is designed to control for correlation bias. Using the median value gives better results than the directional significance or the combined metric. Blood dataset. Need to filter by FDR then sort by effect size.
stat <- res_bl_df$median
names(stat) <- rownames(res_bl_df)
stat[is.na(stat)] <- 0
tic()
cres <- cameraPR(statistic=stat, index=genesets, use.ranks = FALSE, inter.gene.cor=0.01, sort = TRUE)
cres <- subset(cres,NGenes>4)
cres$FDR <- p.adjust(cres$PValue,method="fdr")
mymedians <- sapply(1:nrow(cres),function(i) {
myset <- genesets[[which(names(genesets) == rownames(cres)[i])]]
mystats <- stat[names(stat) %in% myset]
median(mystats)
})
cres$median <- mymedians
toc() # 1.0 sec## 3.604 sec elapsedsig <- subset(cres,FDR<0.05)
head(sig[order(-abs(sig$median)),],20) %>%
kbl(caption = "BLOOD: Top effect size pathways found with CAMERA after 5% FDR filtering") %>%
kable_paper("hover", full_width = F)| NGenes | Direction | PValue | FDR | median | |
|---|---|---|---|---|---|
| Expression and translocation of olfactory receptors | 359 | Down | 0.0000000 | 0.0000557 | -0.8031927 |
| Class C/3 (Metabotropic glutamate/pheromone receptors) | 39 | Down | 0.0000005 | 0.0002979 | -0.8029345 |
| Olfactory Signaling Pathway | 366 | Down | 0.0000001 | 0.0000559 | -0.7970944 |
| Sensory perception of sweet, bitter, and umami (glutamate) taste | 41 | Down | 0.0000123 | 0.0053857 | -0.7224652 |
| Aspirin ADME | 44 | Down | 0.0000705 | 0.0085143 | -0.6442641 |
| Glucuronidation | 25 | Down | 0.0000223 | 0.0053857 | -0.6094985 |
| Sensory perception of taste | 47 | Down | 0.0002054 | 0.0150070 | -0.5983857 |
| Sensory Perception | 570 | Down | 0.0000196 | 0.0053857 | -0.5921152 |
| Metabolism of RNA | 639 | Up | 0.0007046 | 0.0367915 | -0.2416861 |
| Processing of Capped Intron-Containing Pre-mRNA | 232 | Up | 0.0004330 | 0.0246030 | -0.2258106 |
| mRNA Splicing - Major Pathway | 174 | Up | 0.0003687 | 0.0215874 | -0.2159036 |
| mRNA Splicing | 184 | Up | 0.0001359 | 0.0115705 | -0.2019545 |
| Translation | 262 | Up | 0.0003148 | 0.0196209 | -0.1866647 |
| Influenza Infection | 148 | Up | 0.0005360 | 0.0295889 | -0.1749798 |
| mRNA Splicing - Minor Pathway | 51 | Up | 0.0006806 | 0.0365282 | -0.1682665 |
| Regulation of expression of SLITs and ROBOs | 159 | Up | 0.0003659 | 0.0215874 | -0.1597353 |
| rRNA processing | 186 | Up | 0.0002097 | 0.0150070 | -0.1578202 |
| Influenza Viral RNA Transcription and Replication | 129 | Up | 0.0002691 | 0.0173319 | -0.1559051 |
| rRNA processing in the nucleus and cytosol | 180 | Up | 0.0001005 | 0.0092445 | -0.1512278 |
| Selenoamino acid metabolism | 103 | Up | 0.0001377 | 0.0115705 | -0.1483371 |
nrow(subset(cres,FDR<0.05 & Direction=="Up"))## [1] 31nrow(subset(cres,FDR<0.05 & Direction=="Down"))## [1] 8cres_bl <- cresNow drilling down to the top set.
myset <- rownames(head(sig[order(-abs(sig$median)),],1))
myset_members <- genesets[[which(names(genesets) == myset)]]
mystats <- stat[which(names(stat) %in% myset_members)]
myl <- list("allgenes"=stat,myset=mystats)
boxplot(myl,cex=0)
library(vioplot)## Loading required package: sm## Package 'sm', version 2.2-5.7: type help(sm) for summary information## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericvioplot(myl,main=myset)
mygenes <- head(mystats[order(mystats)],5)
mygenes## OR2L2 OR5AK2 OR4C11 OR5T3 OR5R1
## -2.861307 -2.405921 -2.359792 -2.336945 -2.259373myprobes <- unique(unlist(sets[which(names(sets) %in% names(mygenes))]))
myprobedat <- res_bl_top[which(rownames(res_bl_top) %in% myprobes),]
myprobedat %>%
kbl(caption = "Probe data for top 5 genes") %>%
kable_paper("hover", full_width = F)| logFC | AveExpr | t | P.Value | adj.P.Val | B | |
|---|---|---|---|---|---|---|
| cg00938877 | -0.0813302 | 2.8547659 | -3.111247 | 0.0090933 | 0.8764716 | -2.674653 |
| cg08959365 | -0.0535718 | 2.6306172 | -2.861307 | 0.0144397 | 0.8764716 | -3.113687 |
| cg04277698 | -0.0944025 | 2.2013011 | -2.657212 | 0.0210447 | 0.8764716 | -3.468232 |
| cg23984446 | -0.0365309 | 3.5304323 | -2.430747 | 0.0318648 | 0.8764716 | -3.854519 |
| cg14696931 | -0.0743728 | 2.3521299 | -2.381095 | 0.0348756 | 0.8764716 | -3.937880 |
| cg15141401 | -0.0855296 | 1.5597086 | -2.359792 | 0.0362497 | 0.8764716 | -3.973475 |
| cg14284630 | -0.0444939 | 2.0070866 | -2.336945 | 0.0377813 | 0.8764716 | -4.011531 |
| cg17381444 | -0.0404639 | 2.5490545 | -2.316046 | 0.0392364 | 0.8764716 | -4.046229 |
| cg15434232 | -0.0611575 | 0.7671502 | -2.202700 | 0.0481086 | 0.8764716 | -4.232419 |
| cg07804888 | 0.0168582 | 2.7121774 | 1.406162 | 0.1852768 | 0.8764716 | -5.401324 |
| cg20184856 | -0.0877620 | 2.8980060 | -1.227698 | 0.2433189 | 0.8764716 | -5.616157 |
Guthrie dataset.
stat <- res_gu_df$median
names(stat) <- rownames(res_gu_df)
stat[is.na(stat)] <- 0
tic()
cres <- cameraPR(statistic=stat, index=genesets, use.ranks = FALSE, inter.gene.cor=0.01, sort = TRUE)
cres <- subset(cres,NGenes>4)
cres$FDR <- p.adjust(cres$PValue,method="fdr")
mymedians <- sapply(1:nrow(cres),function(i) {
myset <- genesets[[which(names(genesets) == rownames(cres)[i])]]
mystats <- stat[names(stat) %in% myset]
median(mystats)
})
cres$median <- mymedians
toc() # 1.0 sec## 6.877 sec elapsedsig <- subset(cres,FDR<0.05)
head(sig[order(-abs(sig$median)),],20) %>%
kbl(caption = "GUTHRIE: Top effect size pathways found with CAMERA after 5% FDR filtering") %>%
kable_paper("hover", full_width = F)| NGenes | Direction | PValue | FDR | median | |
|---|---|---|---|---|---|
| Class C/3 (Metabotropic glutamate/pheromone receptors) | 39 | Down | 7.67e-05 | 0.0370595 | -0.4694532 |
| Expression and translocation of olfactory receptors | 356 | Down | 0.00e+00 | 0.0000504 | -0.4371060 |
| Olfactory Signaling Pathway | 363 | Down | 1.00e-07 | 0.0000504 | -0.4354951 |
| Sensory Perception | 567 | Down | 5.50e-06 | 0.0035182 | -0.3662568 |
cres_gu <- cresCompare CAMERA data
Compare median set score for both contrasts. In the graph below, the larger points represent larger gene sets. Red fill means significant in blood contrast. Blue ring means significant in guthrie contrast.
m4 <- merge(cres_bl,cres_gu,by=0)
sigx <- subset(m4,FDR.x<0.05)
sigy <- subset(m4,FDR.y<0.05)
plot(m4$median.x,m4$median.y,col="gray",pch=19,cex=log10(m4$NGenes.x),
xlab="blood",ylab="Guthrie",main="median t-stat")
grid()
points(sigx$median.x,sigx$median.y,col="red",pch=19,cex=log10(sigx$NGenes.x))
points(sigy$median.x,sigy$median.y,col="blue",pch=1,cex=log10(sigy$NGenes.x))
abline(v=0,h=0,lty=2,lwd=3)
par(mfrow=c(1,1))
pdf("camera_pw_ados.pdf")
plot(m4$median.x,m4$median.y,col="gray",pch=19,cex=log10(m4$NGenes.x),
xlab="blood",ylab="Guthrie",main="median t-stat")
grid()
points(sigx$median.x,sigx$median.y,col="red",pch=19,cex=log10(sigx$NGenes.x))
points(sigy$median.x,sigy$median.y,col="blue",pch=1,cex=log10(sigy$NGenes.x))
abline(v=0,h=0,lty=2,lwd=3)
dev.off()## png
## 2subset(m4,FDR.x<0.05 & FDR.y<0.05 ) %>%
kbl(caption = "Pathways significant in blood (x) and guthrie cards (y) at 5%FDR") %>%
kable_paper("hover", full_width = F)| Row.names | NGenes.x | Direction.x | PValue.x | FDR.x | median.x | NGenes.y | Direction.y | PValue.y | FDR.y | median.y | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 262 | Class C/3 (Metabotropic glutamate/pheromone receptors) | 39 | Down | 5.00e-07 | 0.0002979 | -0.8029345 | 39 | Down | 7.67e-05 | 0.0370595 | -0.4694532 |
| 522 | Expression and translocation of olfactory receptors | 359 | Down | 0.00e+00 | 0.0000557 | -0.8031927 | 356 | Down | 0.00e+00 | 0.0000504 | -0.4371060 |
| 1096 | Olfactory Signaling Pathway | 366 | Down | 1.00e-07 | 0.0000559 | -0.7970944 | 363 | Down | 1.00e-07 | 0.0000504 | -0.4354951 |
| 1512 | Sensory Perception | 570 | Down | 1.96e-05 | 0.0053857 | -0.5921152 | 567 | Down | 5.50e-06 | 0.0035182 | -0.3662568 |
make scatterplots of the top sets
XMAX=max(res$ranked_profile[,1])
XMIN=min(res$ranked_profile[,1])
YMAX=max(res$ranked_profile[,2])
YMIN=min(res$ranked_profile[,2])
plot(res$detailed_sets[[1]], pch=19,cex=2,col="lightgray",
xlim=c(XMIN,XMAX),ylim=c(YMIN,YMAX),
xlab="blood",ylab="Guthrie")
text(res$detailed_sets[[1]],labels=rownames(res$detailed_sets[[1]]))
mtext(names(res$detailed_sets)[[1]])
grid() ; abline(v=0,h=0,lty=2,lwd=3)
pdf("c3_metabotropic_ados.pdf")
plot(res$detailed_sets[[1]], pch=19,cex=2,col="lightgray",
xlim=c(XMIN,XMAX),ylim=c(YMIN,YMAX),
xlab="blood",ylab="Guthrie")
text(res$detailed_sets[[1]],labels=rownames(res$detailed_sets[[1]]))
mtext(names(res$detailed_sets)[[1]])
grid() ; abline(v=0,h=0,lty=2,lwd=3)
dev.off()## png
## 2Make heatmap of top set.
topset<- rownames(res$detailed_sets[[1]])
probes <- unlist(sets[which(names(sets) %in% topset)])
topm <- bl_mvals[which(rownames(bl_mvals) %in% probes),]
colfunc <- colorRampPalette(c("yellow", "red"))
col_pal <- colfunc(length(min(bl_design[,ncol(bl_design)]):max(bl_design[,ncol(bl_design)])))
mycols <- col_pal[bl_design[,ncol(bl_design)]]
colfunc2 <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(topm,scale="row",trace="none",ColSideColors=mycols,col=colfunc2(25),
main="C/3 receptor methylation in blood")
topset<- rownames(res$detailed_sets[[1]])
probes <- unlist(sets[which(names(sets) %in% topset)])
topm <- gu_mvals[which(rownames(gu_mvals) %in% probes),]
colfunc <- colorRampPalette(c("yellow", "red"))
col_pal <- colfunc(length(min(gu_design[,ncol(gu_design)]):max(gu_design[,ncol(gu_design)])))
mycols <- col_pal[gu_design[,ncol(gu_design)]]
colfunc2 <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(topm,scale="row",trace="none",ColSideColors=mycols,col=colfunc2(25),
main="C/3 receptor methylation in Guthrie cards")
Distinguishing ASD from birth
Here I will use MDS plot to cluster samples.
First I will use all probes, then I will focus on a subset of candidate probes.
plot(cmdscale(dist(t(bl_mvals))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "n",main="Blood")
mtext("Numbers indicate ADOS score")
text(cmdscale(dist(t(bl_mvals))), labels=bl_design[,ncol(bl_design)], )
plot(cmdscale(dist(t(gu_mvals))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "n",main="Guthrie card")
mtext("Numbers indicate ADOS score")
text(cmdscale(dist(t(gu_mvals))), labels=gu_design[,ncol(gu_design)], )
There is no apparent clustering by ADOS score, indicating a subtle effect.
Now I will examine the probes associated with the pathway of interest.
myprobes## [1] "cg20184856" "cg00938877" "cg08959365" "cg10435092" "cg15141401"
## [6] "cg23984446" "cg22123267" "cg14696931" "cg15434232" "cg04277698"
## [11] "cg07804888" "cg17381444" "cg14284630" "cg17626801"myprobedat <- bl_mvals[which(rownames(bl_mvals) %in% myprobes),]
plot(cmdscale(dist(t(myprobedat))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "n",main="Blood")
text(cmdscale(dist(t(myprobedat))), labels=bl_design[,ncol(bl_design)], )
myprobedat <- gu_mvals[which(rownames(gu_mvals) %in% myprobes),]
plot(cmdscale(dist(t(myprobedat))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "n",main="Guthrie card")
text(cmdscale(dist(t(myprobedat))), labels=gu_design[,ncol(gu_design)], )
Again, there is no apparent clustering by ADOS score. This indicates that the top probes cannot be used diagnostically.
myprobes2 <- head(bl_lim$Name,50)
myprobedat2 <- bl_mvals[which(rownames(bl_mvals) %in% myprobes2),]
plot(cmdscale(dist(t(myprobedat2))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "n",main="Blood")
text(cmdscale(dist(t(myprobedat2))), labels=bl_design[,ncol(bl_design)], )
myprobes2 <- head(gu_lim$Name,50)
myprobedat2 <- gu_mvals[which(rownames(gu_mvals) %in% myprobes2),]
plot(cmdscale(dist(t(myprobedat2))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "n", main="Guthrie")
text(cmdscale(dist(t(myprobedat2))), labels=gu_design[,ncol(gu_design)], )
MDS analysis of top probes identified without correction for covariates.
bl_design2 <- bl_design[,c(1,10)]
fit <- lmFit(bl_mvals, bl_design2)
fit <- eBayes(fit)
summary(decideTests(fit))## (Intercept) ADOS
## Down 248755 0
## NotSig 39844 802647
## Up 514048 0top <- topTable(fit,coef="ADOS",num=Inf, sort.by = "P")
head(top)## logFC AveExpr t P.Value adj.P.Val B
## cg13784557 0.04779472 -4.739442 5.456349 2.595342e-05 0.9999998 2.647442
## cg03716591 0.08873078 -3.321034 5.282177 3.833417e-05 0.9999998 2.267504
## cg11078128 -0.05723514 2.941155 -5.154247 5.114886e-05 0.9999998 1.986608
## cg06091288 0.44586721 -1.482468 5.126156 5.450383e-05 0.9999998 1.924737
## cg00320749 0.07116909 -3.962859 5.116515 5.570628e-05 0.9999998 1.903487
## cg00572586 -0.07365197 -3.390343 -5.061720 6.307166e-05 0.9999998 1.782568myprobes3 <- rownames(top)[1:20]
myprobedat3 <- bl_mvals[which(rownames(bl_mvals) %in% myprobes3),]
plot(cmdscale(dist(t(myprobedat3))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "n",main="Blood")
text(cmdscale(dist(t(myprobedat3))), labels=bl_design[,ncol(bl_design)], )
gu_design2 <- gu_design[,c(1,12)]
fit <- lmFit(gu_mvals, gu_design2)
fit <- eBayes(fit)
summary(decideTests(fit))## (Intercept) ADOS
## Down 251749 0
## NotSig 57584 790658
## Up 481325 0top <- topTable(fit,coef="ADOS",num=Inf, sort.by = "P")
head(top)## logFC AveExpr t P.Value adj.P.Val B
## cg24046505 0.05799757 -1.9843389 5.232968 2.418307e-05 0.9999981 2.662615
## cg04098297 0.06458003 2.5127605 5.024096 4.082524e-05 0.9999981 2.158205
## cg00779672 0.06285769 2.1754652 4.996222 4.378776e-05 0.9999981 2.090748
## cg09303399 -0.08027977 -0.3905299 -4.884562 5.799307e-05 0.9999981 1.820269
## cg15310518 -0.04903224 -0.3222181 -4.883002 5.822148e-05 0.9999981 1.816486
## cg01692639 -0.05839269 -3.8720322 -4.783392 7.483579e-05 0.9999981 1.574938myprobes3 <- rownames(top)[1:20]
myprobedat3 <- gu_mvals[which(rownames(gu_mvals) %in% myprobes3),]
plot(cmdscale(dist(t(myprobedat3))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "n",main="Guthrie")
text(cmdscale(dist(t(myprobedat3))), labels=gu_design[,ncol(gu_design)], )
Again it was unsuccessful.
Try ebGSEA
ebGSEA uses globaltest to identify DMGs from infinium data. Some points to note:
only 1 sided test - does not distinguish between up and down methylated genes.
cannot use complex experiment designs.
Therefore I modify how the ebGSEA code uses globaltest.
item mapEIDto850k.lv : A list mapping Entrez Gene ID to 850k probes
#bl_gt <- doGT(bl_design,bl_mvals,array="850k",ncores=16)
data.m <- bl_mvals
subsets <- mclapply(mapEIDto850k.lv,intersect,rownames(data.m),mc.cores = 16)
nrep.v <- unlist(lapply(subsets,length))
selG.idx <- which(nrep.v>0)
gt.o <- gt(response=bl_design[,"ADOS"], alternative=t(data.m), model="linear",
directional = TRUE, standardize = FALSE, permutations = 0,
subsets=subsets[selG.idx],trace=F)
gt_res <- result(gt.o)
gt_res <- gt_res[order(gt_res$`p-value`),]Gene set level enrichment - ADOS
genesets <- gmt_import("ReactomePathways.gmt")
cores = 8Aggregate mval to genes and examine pathway enrichment in blood and guthries.
tic()
bl_gsmea <- gsmea(mval=bl_mvals,design=bl_design,probesets=sets,genesets=genesets,cores=8)
toc()## 198.519 sec elapsedbl_gsmea$res <- subset(bl_gsmea$res,nprobes>=5)
bl_gsmea$res$FDR <- p.adjust(bl_gsmea$res$PValue,method="fdr")
head(bl_gsmea$res,20)## nprobes
## RAS processing 19
## Signaling by MAPK mutants 6
## LDL remodeling 6
## Packaging Of Telomere Ends 7
## Removal of the Flap Intermediate from the C-strand 17
## POLB-Dependent Long Patch Base Excision Repair 7
## Integration of provirus 9
## Early Phase of HIV Life Cycle 14
## Protein repair 6
## 2-LTR circle formation 7
## Phase 1 - inactivation of fast Na+ channels 6
## Pyroptosis 26
## Processive synthesis on the C-strand of the telomere 19
## Synthesis of Ketone Bodies 8
## tRNA modification in the mitochondrion 7
## APOBEC3G mediated resistance to HIV-1 infection 5
## Synthesis of diphthamide-EEF2 8
## PCNA-Dependent Long Patch Base Excision Repair 21
## Resolution of AP sites via the multiple-nucleotide patch replacement pathway 24
## Activation of RAS in B cells 5
## median
## RAS processing -0.6367022
## Signaling by MAPK mutants -0.8910984
## LDL remodeling -1.4384307
## Packaging Of Telomere Ends -0.9230999
## Removal of the Flap Intermediate from the C-strand -0.9620088
## POLB-Dependent Long Patch Base Excision Repair -0.9952393
## Integration of provirus -0.7583022
## Early Phase of HIV Life Cycle -0.6760405
## Protein repair -0.9225705
## 2-LTR circle formation -0.7583022
## Phase 1 - inactivation of fast Na+ channels -0.5750352
## Pyroptosis -0.7617448
## Processive synthesis on the C-strand of the telomere -0.9952393
## Synthesis of Ketone Bodies -1.2918676
## tRNA modification in the mitochondrion -0.2242023
## APOBEC3G mediated resistance to HIV-1 infection -1.1021444
## Synthesis of diphthamide-EEF2 -0.8318361
## PCNA-Dependent Long Patch Base Excision Repair -0.9538704
## Resolution of AP sites via the multiple-nucleotide patch replacement pathway -0.9745548
## Activation of RAS in B cells -0.6367022
## PValue
## RAS processing 0.003037974
## Signaling by MAPK mutants 0.005847454
## LDL remodeling 0.006521048
## Packaging Of Telomere Ends 0.006747719
## Removal of the Flap Intermediate from the C-strand 0.010727635
## POLB-Dependent Long Patch Base Excision Repair 0.011627318
## Integration of provirus 0.012098183
## Early Phase of HIV Life Cycle 0.012385914
## Protein repair 0.015444495
## 2-LTR circle formation 0.015522860
## Phase 1 - inactivation of fast Na+ channels 0.015581495
## Pyroptosis 0.016724018
## Processive synthesis on the C-strand of the telomere 0.017142568
## Synthesis of Ketone Bodies 0.017952402
## tRNA modification in the mitochondrion 0.018124031
## APOBEC3G mediated resistance to HIV-1 infection 0.019387763
## Synthesis of diphthamide-EEF2 0.019415129
## PCNA-Dependent Long Patch Base Excision Repair 0.019466886
## Resolution of AP sites via the multiple-nucleotide patch replacement pathway 0.022145109
## Activation of RAS in B cells 0.023422280
## FDR
## RAS processing 0.5514271
## Signaling by MAPK mutants 0.5514271
## LDL remodeling 0.5514271
## Packaging Of Telomere Ends 0.5514271
## Removal of the Flap Intermediate from the C-strand 0.5514271
## POLB-Dependent Long Patch Base Excision Repair 0.5514271
## Integration of provirus 0.5514271
## Early Phase of HIV Life Cycle 0.5514271
## Protein repair 0.5514271
## 2-LTR circle formation 0.5514271
## Phase 1 - inactivation of fast Na+ channels 0.5514271
## Pyroptosis 0.5514271
## Processive synthesis on the C-strand of the telomere 0.5514271
## Synthesis of Ketone Bodies 0.5514271
## tRNA modification in the mitochondrion 0.5514271
## APOBEC3G mediated resistance to HIV-1 infection 0.5514271
## Synthesis of diphthamide-EEF2 0.5514271
## PCNA-Dependent Long Patch Base Excision Repair 0.5514271
## Resolution of AP sites via the multiple-nucleotide patch replacement pathway 0.5514271
## Activation of RAS in B cells 0.5514271tic()
gu_gsmea <- gsmea(mval=gu_mvals,design=gu_design,probesets=sets,genesets=genesets,cores=8)
toc()## 197.278 sec elapsedgu_gsmea$res <- subset(gu_gsmea$res,nprobes>=5)
gu_gsmea$res$FDR <- p.adjust(gu_gsmea$res$PValue,method="fdr")
head(gu_gsmea$res,20)## nprobes
## APC-Cdc20 mediated degradation of Nek2A 26
## Inactivation of APC/C via direct inhibition of the APC/C complex 21
## Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components 21
## KSRP (KHSRP) binds and destabilizes mRNA 17
## Signaling by NTRK3 (TRKC) 17
## Phosphorylation of Emi1 6
## APC/C:Cdc20 mediated degradation of Cyclin B 24
## POLB-Dependent Long Patch Base Excision Repair 7
## Phosphorylation of the APC/C 20
## RUNX3 regulates WNT signaling 8
## Transport and synthesis of PAPS 6
## Activated NTRK3 signals through PI3K 6
## Binding of TCF/LEF:CTNNB1 to target gene promoters 8
## Diseases associated with surfactant metabolism 9
## Negative regulation of TCF-dependent signaling by DVL-interacting proteins 5
## Formation of apoptosome 10
## Regulation of the apoptosome activity 10
## Serine biosynthesis 9
## Sodium-coupled phosphate cotransporters 5
## CLEC7A/inflammasome pathway 6
## median
## APC-Cdc20 mediated degradation of Nek2A 0.3593637
## Inactivation of APC/C via direct inhibition of the APC/C complex 0.3512109
## Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components 0.3512109
## KSRP (KHSRP) binds and destabilizes mRNA -0.7596917
## Signaling by NTRK3 (TRKC) -0.6409496
## Phosphorylation of Emi1 0.5763514
## APC/C:Cdc20 mediated degradation of Cyclin B 0.3355649
## POLB-Dependent Long Patch Base Excision Repair 0.3560613
## Phosphorylation of the APC/C 0.3355649
## RUNX3 regulates WNT signaling -1.0689451
## Transport and synthesis of PAPS -0.9923302
## Activated NTRK3 signals through PI3K -1.1584291
## Binding of TCF/LEF:CTNNB1 to target gene promoters -1.2414620
## Diseases associated with surfactant metabolism -0.8414110
## Negative regulation of TCF-dependent signaling by DVL-interacting proteins -0.7514719
## Formation of apoptosome -0.6161805
## Regulation of the apoptosome activity -0.6161805
## Serine biosynthesis 0.7882701
## Sodium-coupled phosphate cotransporters -1.0204257
## CLEC7A/inflammasome pathway 0.7691143
## PValue
## APC-Cdc20 mediated degradation of Nek2A 0.01210049
## Inactivation of APC/C via direct inhibition of the APC/C complex 0.01809397
## Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components 0.01809397
## KSRP (KHSRP) binds and destabilizes mRNA 0.02866691
## Signaling by NTRK3 (TRKC) 0.02939644
## Phosphorylation of Emi1 0.03293071
## APC/C:Cdc20 mediated degradation of Cyclin B 0.03341884
## POLB-Dependent Long Patch Base Excision Repair 0.03829817
## Phosphorylation of the APC/C 0.04069662
## RUNX3 regulates WNT signaling 0.04218950
## Transport and synthesis of PAPS 0.04847664
## Activated NTRK3 signals through PI3K 0.05098899
## Binding of TCF/LEF:CTNNB1 to target gene promoters 0.05245429
## Diseases associated with surfactant metabolism 0.05449754
## Negative regulation of TCF-dependent signaling by DVL-interacting proteins 0.06702115
## Formation of apoptosome 0.06938447
## Regulation of the apoptosome activity 0.06938447
## Serine biosynthesis 0.07735325
## Sodium-coupled phosphate cotransporters 0.07848826
## CLEC7A/inflammasome pathway 0.08895789
## FDR
## APC-Cdc20 mediated degradation of Nek2A 0.9999532
## Inactivation of APC/C via direct inhibition of the APC/C complex 0.9999532
## Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components 0.9999532
## KSRP (KHSRP) binds and destabilizes mRNA 0.9999532
## Signaling by NTRK3 (TRKC) 0.9999532
## Phosphorylation of Emi1 0.9999532
## APC/C:Cdc20 mediated degradation of Cyclin B 0.9999532
## POLB-Dependent Long Patch Base Excision Repair 0.9999532
## Phosphorylation of the APC/C 0.9999532
## RUNX3 regulates WNT signaling 0.9999532
## Transport and synthesis of PAPS 0.9999532
## Activated NTRK3 signals through PI3K 0.9999532
## Binding of TCF/LEF:CTNNB1 to target gene promoters 0.9999532
## Diseases associated with surfactant metabolism 0.9999532
## Negative regulation of TCF-dependent signaling by DVL-interacting proteins 0.9999532
## Formation of apoptosome 0.9999532
## Regulation of the apoptosome activity 0.9999532
## Serine biosynthesis 0.9999532
## Sodium-coupled phosphate cotransporters 0.9999532
## CLEC7A/inflammasome pathway 0.9999532Gene set level enrichment - diagnosis
tic()
bl_gsmea <- gsmea(mval=bl_mvals,design=bl_design,probesets=sets,genesets=genesets,cores=8)
toc()## 199.292 sec elapsedbl_gsmea$res <- subset(bl_gsmea$res,nprobes>=5)
bl_gsmea$res$FDR <- p.adjust(bl_gsmea$res$PValue,method="fdr")
head(bl_gsmea$res,20)## nprobes
## RAS processing 19
## Signaling by MAPK mutants 6
## LDL remodeling 6
## Packaging Of Telomere Ends 7
## Removal of the Flap Intermediate from the C-strand 17
## POLB-Dependent Long Patch Base Excision Repair 7
## Integration of provirus 9
## Early Phase of HIV Life Cycle 14
## Protein repair 6
## 2-LTR circle formation 7
## Phase 1 - inactivation of fast Na+ channels 6
## Pyroptosis 26
## Processive synthesis on the C-strand of the telomere 19
## Synthesis of Ketone Bodies 8
## tRNA modification in the mitochondrion 7
## APOBEC3G mediated resistance to HIV-1 infection 5
## Synthesis of diphthamide-EEF2 8
## PCNA-Dependent Long Patch Base Excision Repair 21
## Resolution of AP sites via the multiple-nucleotide patch replacement pathway 24
## Activation of RAS in B cells 5
## median
## RAS processing -0.6367022
## Signaling by MAPK mutants -0.8910984
## LDL remodeling -1.4384307
## Packaging Of Telomere Ends -0.9230999
## Removal of the Flap Intermediate from the C-strand -0.9620088
## POLB-Dependent Long Patch Base Excision Repair -0.9952393
## Integration of provirus -0.7583022
## Early Phase of HIV Life Cycle -0.6760405
## Protein repair -0.9225705
## 2-LTR circle formation -0.7583022
## Phase 1 - inactivation of fast Na+ channels -0.5750352
## Pyroptosis -0.7617448
## Processive synthesis on the C-strand of the telomere -0.9952393
## Synthesis of Ketone Bodies -1.2918676
## tRNA modification in the mitochondrion -0.2242023
## APOBEC3G mediated resistance to HIV-1 infection -1.1021444
## Synthesis of diphthamide-EEF2 -0.8318361
## PCNA-Dependent Long Patch Base Excision Repair -0.9538704
## Resolution of AP sites via the multiple-nucleotide patch replacement pathway -0.9745548
## Activation of RAS in B cells -0.6367022
## PValue
## RAS processing 0.003037974
## Signaling by MAPK mutants 0.005847454
## LDL remodeling 0.006521048
## Packaging Of Telomere Ends 0.006747719
## Removal of the Flap Intermediate from the C-strand 0.010727635
## POLB-Dependent Long Patch Base Excision Repair 0.011627318
## Integration of provirus 0.012098183
## Early Phase of HIV Life Cycle 0.012385914
## Protein repair 0.015444495
## 2-LTR circle formation 0.015522860
## Phase 1 - inactivation of fast Na+ channels 0.015581495
## Pyroptosis 0.016724018
## Processive synthesis on the C-strand of the telomere 0.017142568
## Synthesis of Ketone Bodies 0.017952402
## tRNA modification in the mitochondrion 0.018124031
## APOBEC3G mediated resistance to HIV-1 infection 0.019387763
## Synthesis of diphthamide-EEF2 0.019415129
## PCNA-Dependent Long Patch Base Excision Repair 0.019466886
## Resolution of AP sites via the multiple-nucleotide patch replacement pathway 0.022145109
## Activation of RAS in B cells 0.023422280
## FDR
## RAS processing 0.5514271
## Signaling by MAPK mutants 0.5514271
## LDL remodeling 0.5514271
## Packaging Of Telomere Ends 0.5514271
## Removal of the Flap Intermediate from the C-strand 0.5514271
## POLB-Dependent Long Patch Base Excision Repair 0.5514271
## Integration of provirus 0.5514271
## Early Phase of HIV Life Cycle 0.5514271
## Protein repair 0.5514271
## 2-LTR circle formation 0.5514271
## Phase 1 - inactivation of fast Na+ channels 0.5514271
## Pyroptosis 0.5514271
## Processive synthesis on the C-strand of the telomere 0.5514271
## Synthesis of Ketone Bodies 0.5514271
## tRNA modification in the mitochondrion 0.5514271
## APOBEC3G mediated resistance to HIV-1 infection 0.5514271
## Synthesis of diphthamide-EEF2 0.5514271
## PCNA-Dependent Long Patch Base Excision Repair 0.5514271
## Resolution of AP sites via the multiple-nucleotide patch replacement pathway 0.5514271
## Activation of RAS in B cells 0.5514271tic()
gu_gsmea <- gsmea(mval=gu_mvals,design=gu_design,probesets=sets,genesets=genesets,cores=8)
toc()## 197.187 sec elapsedgu_gsmea$res <- subset(gu_gsmea$res,nprobes>=5)
gu_gsmea$res$FDR <- p.adjust(gu_gsmea$res$PValue,method="fdr")
head(gu_gsmea$res,20)## nprobes
## APC-Cdc20 mediated degradation of Nek2A 26
## Inactivation of APC/C via direct inhibition of the APC/C complex 21
## Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components 21
## KSRP (KHSRP) binds and destabilizes mRNA 17
## Signaling by NTRK3 (TRKC) 17
## Phosphorylation of Emi1 6
## APC/C:Cdc20 mediated degradation of Cyclin B 24
## POLB-Dependent Long Patch Base Excision Repair 7
## Phosphorylation of the APC/C 20
## RUNX3 regulates WNT signaling 8
## Transport and synthesis of PAPS 6
## Activated NTRK3 signals through PI3K 6
## Binding of TCF/LEF:CTNNB1 to target gene promoters 8
## Diseases associated with surfactant metabolism 9
## Negative regulation of TCF-dependent signaling by DVL-interacting proteins 5
## Formation of apoptosome 10
## Regulation of the apoptosome activity 10
## Serine biosynthesis 9
## Sodium-coupled phosphate cotransporters 5
## CLEC7A/inflammasome pathway 6
## median
## APC-Cdc20 mediated degradation of Nek2A 0.3593637
## Inactivation of APC/C via direct inhibition of the APC/C complex 0.3512109
## Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components 0.3512109
## KSRP (KHSRP) binds and destabilizes mRNA -0.7596917
## Signaling by NTRK3 (TRKC) -0.6409496
## Phosphorylation of Emi1 0.5763514
## APC/C:Cdc20 mediated degradation of Cyclin B 0.3355649
## POLB-Dependent Long Patch Base Excision Repair 0.3560613
## Phosphorylation of the APC/C 0.3355649
## RUNX3 regulates WNT signaling -1.0689451
## Transport and synthesis of PAPS -0.9923302
## Activated NTRK3 signals through PI3K -1.1584291
## Binding of TCF/LEF:CTNNB1 to target gene promoters -1.2414620
## Diseases associated with surfactant metabolism -0.8414110
## Negative regulation of TCF-dependent signaling by DVL-interacting proteins -0.7514719
## Formation of apoptosome -0.6161805
## Regulation of the apoptosome activity -0.6161805
## Serine biosynthesis 0.7882701
## Sodium-coupled phosphate cotransporters -1.0204257
## CLEC7A/inflammasome pathway 0.7691143
## PValue
## APC-Cdc20 mediated degradation of Nek2A 0.01210049
## Inactivation of APC/C via direct inhibition of the APC/C complex 0.01809397
## Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components 0.01809397
## KSRP (KHSRP) binds and destabilizes mRNA 0.02866691
## Signaling by NTRK3 (TRKC) 0.02939644
## Phosphorylation of Emi1 0.03293071
## APC/C:Cdc20 mediated degradation of Cyclin B 0.03341884
## POLB-Dependent Long Patch Base Excision Repair 0.03829817
## Phosphorylation of the APC/C 0.04069662
## RUNX3 regulates WNT signaling 0.04218950
## Transport and synthesis of PAPS 0.04847664
## Activated NTRK3 signals through PI3K 0.05098899
## Binding of TCF/LEF:CTNNB1 to target gene promoters 0.05245429
## Diseases associated with surfactant metabolism 0.05449754
## Negative regulation of TCF-dependent signaling by DVL-interacting proteins 0.06702115
## Formation of apoptosome 0.06938447
## Regulation of the apoptosome activity 0.06938447
## Serine biosynthesis 0.07735325
## Sodium-coupled phosphate cotransporters 0.07848826
## CLEC7A/inflammasome pathway 0.08895789
## FDR
## APC-Cdc20 mediated degradation of Nek2A 0.9999532
## Inactivation of APC/C via direct inhibition of the APC/C complex 0.9999532
## Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components 0.9999532
## KSRP (KHSRP) binds and destabilizes mRNA 0.9999532
## Signaling by NTRK3 (TRKC) 0.9999532
## Phosphorylation of Emi1 0.9999532
## APC/C:Cdc20 mediated degradation of Cyclin B 0.9999532
## POLB-Dependent Long Patch Base Excision Repair 0.9999532
## Phosphorylation of the APC/C 0.9999532
## RUNX3 regulates WNT signaling 0.9999532
## Transport and synthesis of PAPS 0.9999532
## Activated NTRK3 signals through PI3K 0.9999532
## Binding of TCF/LEF:CTNNB1 to target gene promoters 0.9999532
## Diseases associated with surfactant metabolism 0.9999532
## Negative regulation of TCF-dependent signaling by DVL-interacting proteins 0.9999532
## Formation of apoptosome 0.9999532
## Regulation of the apoptosome activity 0.9999532
## Serine biosynthesis 0.9999532
## Sodium-coupled phosphate cotransporters 0.9999532
## CLEC7A/inflammasome pathway 0.9999532Session information
For reproducibility
sessionInfo()## R version 4.2.2 Patched (2022-11-10 r83330)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 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
##
## attached base packages:
## [1] grid stats4 parallel stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] vioplot_0.4.0
## [2] zoo_1.8-11
## [3] sm_2.2-5.7.1
## [4] pkgload_1.3.2
## [5] GGally_2.1.2
## [6] ggplot2_3.4.1
## [7] gtools_3.9.4
## [8] tibble_3.2.0
## [9] dplyr_1.1.0
## [10] echarts4r_0.4.4
## [11] ebGSEA_0.1.0
## [12] globaltest_5.52.0
## [13] survival_3.5-5
## [14] tictoc_1.1
## [15] RIdeogram_0.2.2
## [16] kableExtra_1.3.4
## [17] data.table_1.14.8
## [18] ENmix_1.34.0
## [19] doParallel_1.0.17
## [20] qqman_0.1.8
## [21] RCircos_1.2.2
## [22] beeswarm_0.4.0
## [23] forestplot_3.1.1
## [24] abind_1.4-5
## [25] checkmate_2.1.0
## [26] reshape2_1.4.4
## [27] gplots_3.1.3
## [28] eulerr_7.0.0
## [29] GEOquery_2.66.0
## [30] RColorBrewer_1.1-3
## [31] IlluminaHumanMethylation450kmanifest_0.4.0
## [32] topconfects_1.14.0
## [33] DMRcatedata_2.16.0
## [34] ExperimentHub_2.6.0
## [35] AnnotationHub_3.6.0
## [36] BiocFileCache_2.6.0
## [37] dbplyr_2.3.1
## [38] DMRcate_2.12.0
## [39] limma_3.54.0
## [40] missMethyl_1.32.0
## [41] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [42] R.utils_2.12.2
## [43] R.oo_1.25.0
## [44] R.methodsS3_1.8.2
## [45] plyr_1.8.8
## [46] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [47] minfi_1.44.0
## [48] bumphunter_1.40.0
## [49] locfit_1.5-9.7
## [50] iterators_1.0.14
## [51] foreach_1.5.2
## [52] Biostrings_2.66.0
## [53] XVector_0.38.0
## [54] SummarizedExperiment_1.28.0
## [55] Biobase_2.58.0
## [56] MatrixGenerics_1.10.0
## [57] matrixStats_0.63.0
## [58] GenomicRanges_1.50.2
## [59] GenomeInfoDb_1.34.6
## [60] IRanges_2.32.0
## [61] S4Vectors_0.36.1
## [62] BiocGenerics_0.44.0
## [63] mitch_1.10.0
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 rtracklayer_1.58.0
## [3] tidyr_1.3.0 bit64_4.0.5
## [5] knitr_1.42 DelayedArray_0.24.0
## [7] rpart_4.1.19 KEGGREST_1.38.0
## [9] RCurl_1.98-1.10 AnnotationFilter_1.22.0
## [11] generics_0.1.3 GenomicFeatures_1.50.3
## [13] preprocessCore_1.60.2 RSQLite_2.3.0
## [15] bit_4.0.5 tzdb_0.3.0
## [17] webshot_0.5.4 xml2_1.3.3
## [19] httpuv_1.6.9 assertthat_0.2.1
## [21] xfun_0.37 hms_1.1.2
## [23] jquerylib_0.1.4 evaluate_0.20
## [25] promises_1.2.0.1 fansi_1.0.4
## [27] restfulr_0.0.15 scrime_1.3.5
## [29] progress_1.2.2 caTools_1.18.2
## [31] readxl_1.4.2 DBI_1.1.3
## [33] geneplotter_1.76.0 htmlwidgets_1.6.2
## [35] reshape_0.8.9 purrr_1.0.1
## [37] ellipsis_0.3.2 backports_1.4.1
## [39] permute_0.9-7 calibrate_1.7.7
## [41] grImport2_0.2-0 annotate_1.76.0
## [43] biomaRt_2.54.0 deldir_1.0-6
## [45] sparseMatrixStats_1.10.0 vctrs_0.6.0
## [47] ensembldb_2.22.0 withr_2.5.0
## [49] cachem_1.0.7 Gviz_1.42.0
## [51] BSgenome_1.66.2 GenomicAlignments_1.34.0
## [53] prettyunits_1.1.1 mclust_6.0.0
## [55] svglite_2.1.1 cluster_2.1.4
## [57] RPMM_1.25 lazyeval_0.2.2
## [59] crayon_1.5.2 genefilter_1.80.3
## [61] labeling_0.4.2 edgeR_3.40.2
## [63] pkgconfig_2.0.3 nlme_3.1-162
## [65] ProtGenerics_1.30.0 nnet_7.3-18
## [67] rlang_1.1.0 lifecycle_1.0.3
## [69] filelock_1.0.2 dichromat_2.0-0.1
## [71] rsvg_2.4.0 cellranger_1.1.0
## [73] rngtools_1.5.2 base64_2.0.1
## [75] Matrix_1.5-3 Rhdf5lib_1.20.0
## [77] base64enc_0.1-3 viridisLite_0.4.1
## [79] png_0.1-8 rjson_0.2.21
## [81] bitops_1.0-7 KernSmooth_2.23-20
## [83] rhdf5filters_1.10.0 blob_1.2.4
## [85] DelayedMatrixStats_1.20.0 doRNG_1.8.6
## [87] stringr_1.5.0 nor1mix_1.3-0
## [89] readr_2.1.4 jpeg_0.1-10
## [91] scales_1.2.1 memoise_2.0.1
## [93] magrittr_2.0.3 zlibbioc_1.44.0
## [95] compiler_4.2.2 BiocIO_1.8.0
## [97] illuminaio_0.40.0 Rsamtools_2.14.0
## [99] cli_3.6.0 DSS_2.46.0
## [101] htmlTable_2.4.1 Formula_1.2-5
## [103] MASS_7.3-58.3 tidyselect_1.2.0
## [105] stringi_1.7.12 highr_0.10
## [107] yaml_2.3.7 askpass_1.1
## [109] latticeExtra_0.6-30 sass_0.4.5
## [111] VariantAnnotation_1.44.0 tools_4.2.2
## [113] rstudioapi_0.14 foreign_0.8-84
## [115] bsseq_1.34.0 gridExtra_2.3
## [117] farver_2.1.1 digest_0.6.31
## [119] BiocManager_1.30.20 shiny_1.7.4
## [121] quadprog_1.5-8 Rcpp_1.0.10
## [123] siggenes_1.72.0 BiocVersion_3.16.0
## [125] later_1.3.0 org.Hs.eg.db_3.16.0
## [127] httr_1.4.5 AnnotationDbi_1.60.0
## [129] biovizBase_1.46.0 colorspace_2.1-0
## [131] rvest_1.0.3 XML_3.99-0.13
## [133] splines_4.2.2 statmod_1.5.0
## [135] kpmt_0.1.0 multtest_2.54.0
## [137] systemfonts_1.0.4 xtable_1.8-4
## [139] jsonlite_1.8.4 dynamicTreeCut_1.63-1
## [141] R6_2.5.1 Hmisc_5.0-1
## [143] pillar_1.8.1 htmltools_0.5.4
## [145] mime_0.12 glue_1.6.2
## [147] fastmap_1.1.1 BiocParallel_1.32.5
## [149] interactiveDisplayBase_1.36.0 beanplot_1.3.1
## [151] codetools_0.2-19 utf8_1.2.3
## [153] lattice_0.20-45 bslib_0.4.2
## [155] curl_5.0.0 openssl_2.0.6
## [157] interp_1.1-3 rmarkdown_2.20
## [159] munsell_0.5.0 rhdf5_2.42.0
## [161] GenomeInfoDbData_1.2.9 HDF5Array_1.26.0
## [163] impute_1.72.3 gtable_0.3.2