ASD Gene Methylation Enrichment Analysis
Mark Ziemann
2022-10-12
Introduction
Previously, Mandhri and I analysed the B-PROOF 450K data, trying to understand whether vitamin supplementation caused changes in gene methylation. We used Limma and some basic analyses, which showed no specific probes with FDR<0.05, nor any DMRs.
In this analysis we will use the principle of Gene Set Enrichment Analysis, applying it to many probes belonging to genes. If the probes are trending in concert, then we can make some judgement about the enrichment of those probes. The statistical test used is the 1 sided t-test, which can be applied to competitive or self contained test types.
library("parallel")
library("dplyr")
library("kableExtra")
library("eulerr")
library("mitch")
library("IlluminaHumanMethylationEPICanno.ilm10b2.hg19")
library("tictoc")Get array annotation data
ann <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b2.hg19)
myann <- data.frame(ann[,c("UCSC_RefGene_Name","Regulatory_Feature_Group")])Load ASD data
Previously Namitha performed this analysis - there weren’t many significant genes so I thought it would be a good idea to use GMEA.
Read in and inspect the data
dm <- read.csv("ASD_blood_top_dmps_onADOS_withintw_Nov27_limma.csv")
head(dm,50) %>%
kbl(caption = "Top significant genes with limma") %>%
kable_paper("hover", full_width = F)| Name | chr | pos | strand | Islands_Name | Relation_to_Island | UCSC_RefGene_Name | UCSC_RefGene_Group | logFC | AveExpr | t | P.Value | adj.P.Val | B |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg07655025 | chr19 | 12835001 |
|
chr19:12833104-12833574 | S_Shore | TNPO2 | TSS200 | -0.0854949 | -4.3361600 | -6.706650 | 0.0000228 | 0.8764716 | 3.0481274 |
| cg18253799 | chr17 | 78069191 |
|
OpenSea | CCDC40 | Body | 0.1005063 | 4.6851023 | 5.976848 | 0.0000670 | 0.8764716 | 2.0368010 | |
| cg25032595 | chr13 | 96204978 |
|
chr13:96204691-96205496 | Island | CLDN10;CLDN10;CLDN10;CLDN10 | 5’UTR;1stExon;Body;Body | -0.0705806 | -4.5040270 | -5.801969 | 0.0000877 | 0.8764716 | 1.7818845 |
| cg16956665 | chr11 | 125479362 |
|
OpenSea | STT3A | Body | 0.1228970 | 0.6660224 | 5.770674 | 0.0000921 | 0.8764716 | 1.7357493 | |
| cg06705930 | chr11 | 31823191 |
|
chr11:31820060-31821416 | S_Shore | PAX6;PAX6;PAX6 | Body;Body;Body | -0.0611646 | 0.6547421 | -5.605520 | 0.0001193 | 0.8764716 | 1.4896913 |
| cg18845598 | chr7 | 143746594 |
|
OpenSea | OR2A5 | TSS1500 | -0.0875782 | 0.5607065 | -5.586211 | 0.0001230 | 0.8764716 | 1.4606389 | |
| cg13069346 | chr13 | 66437506 |
|
OpenSea | -0.0778752 | 3.3712970 | -5.545615 | 0.0001311 | 0.8764716 | 1.3993644 | |||
| cg12322146 | chr3 | 29197502 |
|
OpenSea | -0.0824150 | 3.0448357 | -5.451486 | 0.0001524 | 0.8764716 | 1.2562812 | |||
| cg11568374 | chr2 | 64015338 |
|
OpenSea | -0.0685263 | 2.5345142 | -5.340095 | 0.0001823 | 0.8764716 | 1.0851481 | |||
| cg22492966 | chr10 | 65028929 |
|
OpenSea | JMJD1C;JMJD1C;JMJD1C | Body;1stExon;5’UTR | 0.0655305 | -4.3710723 | 5.329706 | 0.0001853 | 0.8764716 | 1.0690876 | |
| cg09119495 | chr5 | 90252465 |
|
OpenSea | ADGRV1;ADGRV1 | Body;Body | 0.0716463 | -1.4703628 | 5.329690 | 0.0001854 | 0.8764716 | 1.0690644 | |
| cg02812571 | chr16 | 25377575 |
|
OpenSea | 0.0854464 | -3.3133349 | 5.318120 | 0.0001889 | 0.8764716 | 1.0511579 | |||
| cg04156077 | chr3 | 149421196 |
|
OpenSea | WWTR1 | TSS200 | -0.1231433 | 1.9400856 | -5.291017 | 0.0001973 | 0.8764716 | 1.0091305 | |
| cg24911519 | chr7 | 123295209 |
|
OpenSea | LMOD2 | TSS1500 | -0.0670507 | 2.4406162 | -5.246528 | 0.0002121 | 0.8764716 | 0.9398935 | |
| cg20573508 | chr4 | 75578663 |
|
OpenSea | -0.0853904 | 3.0762924 | -5.210624 | 0.0002249 | 0.8764716 | 0.8837929 | |||
| cg01350895 | chr11 | 19719938 |
|
OpenSea | NAV2 | Body | 0.0797398 | 2.1580027 | 5.204202 | 0.0002273 | 0.8764716 | 0.8737372 | |
| cg00368055 | chr4 | 11429689 |
|
chr4:11430235-11431256 | N_Shore | HS3ST1 | 5’UTR | -0.0777281 | -4.0918634 | -5.165993 | 0.0002420 | 0.8764716 | 0.8137773 |
| cg13253692 | chr11 | 92245175 |
|
OpenSea | FAT3 | Body | -0.0705305 | 3.0953779 | -5.131886 | 0.0002559 | 0.8764716 | 0.7600645 | |
| cg02746635 | chr5 | 129723351 |
|
OpenSea | -0.0661462 | 3.1495539 | -5.123567 | 0.0002595 | 0.8764716 | 0.7469363 | |||
| cg26566874 | chr7 | 22602547 |
|
OpenSea | LOC100506178 | TSS1500 | -0.0734783 | -3.5861013 | -5.027003 | 0.0003044 | 0.8764716 | 0.5937751 | |
| cg24827992 | chr6 | 33260021 |
|
chr6:33256539-33256973 | S_Shelf | RGL2;RGL2 | Body;Body | -0.0685574 | 3.8052028 | -5.016494 | 0.0003097 | 0.8764716 | 0.5770218 |
| cg12497004 | chr21 | 34287140 |
|
OpenSea | 0.0848047 | -4.6598131 | 4.993825 | 0.0003216 | 0.8764716 | 0.5408255 | |||
| cg15980797 | chr6 | 137813648 |
|
chr6:137814355-137815202 | N_Shore | OLIG3;OLIG3 | 1stExon;3’UTR | 0.0844454 | -3.5461377 | 4.992543 | 0.0003223 | 0.8764716 | 0.5387767 |
| cg19204328 | chr18 | 66379577 |
|
chr18:66381934-66382575 | N_Shelf | TMX3 | Body | 0.0531719 | 3.7266528 | 4.943507 | 0.0003498 | 0.8764716 | 0.4602054 |
| cg13494636 | chr15 | 89918122 |
|
chr15:89920793-89922768 | N_Shelf | 0.1004856 | -2.1071034 | 4.933684 | 0.0003556 | 0.8764716 | 0.4444221 | ||
| cg27474897 | chr2 | 199480893 |
|
OpenSea | -0.0818969 | 0.9703255 | -4.918467 | 0.0003648 | 0.8764716 | 0.4199445 | |||
| cg01343313 | chr6 | 15736169 |
|
OpenSea | 0.0738354 | -3.6617380 | 4.891154 | 0.0003819 | 0.8764716 | 0.3759221 | |||
| cg16942213 | chr20 | 49104270 |
|
OpenSea | 0.0564384 | 3.0326105 | 4.884606 | 0.0003861 | 0.8764716 | 0.3653527 | |||
| cg08189124 | chr11 | 1980957 |
|
chr11:1983124-1983351 | N_Shelf | -0.0581177 | 2.8185106 | -4.883602 | 0.0003867 | 0.8764716 | 0.3637311 | ||
| cg10435618 | chr10 | 838112 |
|
chr10:839608-840210 | N_Shore | -0.0682292 | 3.1880854 | -4.878911 | 0.0003898 | 0.8764716 | 0.3561547 | ||
| cg09518009 | chr14 | 27243163 |
|
OpenSea | -0.0690079 | 2.4212842 | -4.873288 | 0.0003935 | 0.8764716 | 0.3470674 | |||
| cg16361775 | chr4 | 66567883 |
|
OpenSea | -0.0667119 | 0.9583160 | -4.854893 | 0.0004059 | 0.8764716 | 0.3173101 | |||
| cg10354730 | chr16 | 357249 |
|
chr16:354243-354446 | S_Shelf | AXIN1;AXIN1 | Body;Body | -0.0653508 | 3.4895542 | -4.846812 | 0.0004114 | 0.8764716 | 0.3042207 |
| cg07158816 | chr10 | 24564294 |
|
OpenSea | MIR603;KIAA1217;KIAA1217;KIAA1217 | TSS1500;Body;Body;Body | -0.0660411 | 2.9786330 | -4.845037 | 0.0004127 | 0.8764716 | 0.3013448 | |
| cg05145608 | chr7 | 4757991 |
|
chr7:4757802-4758037 | Island | FOXK1 | Body | -0.0616320 | 3.7940401 | -4.833596 | 0.0004207 | 0.8764716 | 0.2827952 |
| cg07988820 | chr12 | 82153109 |
|
chr12:82152320-82152674 | S_Shore | PPFIA2 | TSS200 | -0.0855187 | -4.5012730 | -4.808826 | 0.0004387 | 0.8764716 | 0.2425682 |
| cg25425801 | chr4 | 87514871 |
|
chr4:87515183-87515735 | N_Shore | PTPN13;PTPN13;PTPN13;PTPN13 | TSS1500;TSS1500;TSS1500;TSS1500 | 0.0678365 | -3.4414351 | 4.798957 | 0.0004460 | 0.8764716 | 0.2265171 |
| cg23348270 | chr2 | 176969612 |
|
chr2:176969217-176969895 | Island | 0.0629276 | -4.6889767 | 4.798547 | 0.0004463 | 0.8764716 | 0.2258488 | ||
| cg01577079 | chr3 | 160389186 |
|
OpenSea | -0.0654179 | 4.1147963 | -4.792741 | 0.0004508 | 0.8764716 | 0.2163991 | |||
| cg25629947 | chr10 | 118609001 |
|
chr10:118608599-118609377 | Island | ENO4 | TSS200 | -0.0635783 | -4.7721423 | -4.768611 | 0.0004696 | 0.8764716 | 0.1770708 |
| cg06484220 | chr9 | 89900174 |
|
OpenSea | -0.0569198 | 1.9686843 | -4.763269 | 0.0004738 | 0.8764716 | 0.1683530 | |||
| cg13866747 | chr13 | 91965808 |
|
OpenSea | -0.0842283 | 1.6722604 | -4.757431 | 0.0004785 | 0.8764716 | 0.1588197 | |||
| cg08873424 | chr6 | 28945492 |
|
OpenSea | -0.0833208 | -3.4190101 | -4.756008 | 0.0004797 | 0.8764716 | 0.1564960 | |||
| cg05890541 | chr4 | 100384195 |
|
OpenSea | 0.0734458 | -4.4514453 | 4.749529 | 0.0004850 | 0.8764716 | 0.1459097 | |||
| cg15207358 | chr15 | 44184056 |
|
OpenSea | FRMD5;FRMD5;FRMD5;FRMD5 | Body;Body;Body;Body | -0.0597480 | 3.9712959 | -4.739788 | 0.0004931 | 0.8764716 | 0.1299828 | |
| cg22467932 | chr15 | 49912175 |
|
chr15:49912985-49913482 | N_Shore | FAM227B;DTWD1;DTWD1 | 5’UTR;TSS1500;TSS1500 | -0.0569199 | 3.9417413 | -4.732563 | 0.0004992 | 0.8764716 | 0.1181624 |
| cg01965598 | chr5 | 53607153 |
|
chr5:53605932-53606798 | S_Shore | ARL15 | TSS1500 | 0.0551092 | -4.6660925 | 4.728086 | 0.0005030 | 0.8764716 | 0.1108327 |
| cg19531247 | chr6 | 11056480 |
|
OpenSea | ELOVL2-AS1;ELOVL2-AS1;ELOVL2-AS1 | Body;Body;Body | -0.1128130 | 3.0913449 | -4.698135 | 0.0005293 | 0.8764716 | 0.0617280 | |
| cg09809932 | chr1 | 6515597 |
|
chr1:6514108-6516325 | Island | ESPN | Body | -0.0727228 | -5.4865756 | -4.697547 | 0.0005298 | 0.8764716 | 0.0607630 |
| cg04835122 | chr2 | 74056056 |
|
chr2:74055941-74056320 | Island | STAMBP;STAMBP;STAMBP | TSS200;TSS200;TSS200 | 0.0720488 | -4.7231146 | 4.687265 | 0.0005392 | 0.8764716 | 0.0438751 |
rownames(dm) <- dm$Name
dm <- dm[,c("UCSC_RefGene_Name","t")]
# histogram of t values
hist(dm$t,breaks=seq(from=-7,to=6,by=1))
Whole gene analysis
# set cores to used for parallel execution
CORES= detectCores()/2
calc_sc <- function(dm) {
gn <- unique(unlist(strsplit( dm$UCSC_RefGene_Name ,";")))
gnl <- strsplit( dm$UCSC_RefGene_Name ,";")
gnl <- mclapply(gnl,unique,mc.cores=CORES)
dm$UCSC_RefGene_Name <- gnl
l <- mclapply(1:nrow(dm), function(i) {
a <- dm[i,]
len <- length(a[[1]][[1]])
tvals <- as.numeric(rep(a[2],len))
genes <- a[[1]][[1]]
data.frame(genes,tvals)
},mc.cores=CORES)
df <- do.call(rbind,l)
gme_res <- mclapply( 1:length(gn), function(i) {
g <- gn[i]
tstats <- df[which(df$genes==g),"tvals"]
myn <- length(tstats)
mymean <- mean(tstats)
mymedian <- median(tstats)
if ( length(tstats) < 2 ) {
pval=1
} else {
wtselfcont <- t.test(tstats)
pval=wtselfcont$p.value
}
res <- c("gene"=g,"nprobes"=myn,"mean"=mymean,"median"=mymedian,
"p-value(sc)"=pval)
} , mc.cores=CORES )
gme_res_df <- do.call(rbind, gme_res)
rownames(gme_res_df) <- gme_res_df[,1]
gme_res_df <- gme_res_df[,-1]
tmp <- apply(gme_res_df,2,as.numeric)
rownames(tmp) <- rownames(gme_res_df)
gme_res_df <- as.data.frame(tmp)
gme_res_df$sig <- -log10(gme_res_df[,4])
gme_res_df <- gme_res_df[order(-gme_res_df$sig),]
gme_res_df$`fdr(sc)` <- p.adjust(gme_res_df$`p-value(sc)`)
out <- list("df"=df,"gme_res_df"=gme_res_df)
return(out)
}
tic()
gme_res_wholegene <- calc_sc(dm)
time2 <- toc() #4core: 239.0 247.0## 111.895 sec elapsedtime2## $tic
## elapsed
## 5854.499
##
## $toc
## elapsed
## 5966.394
##
## $msg
## logical(0)
##
## $callback_msg
## [1] "111.895 sec elapsed"df <- gme_res_wholegene[[1]]
res <- gme_res_wholegene[[2]]
write.table(res,file="gmea_wholegene.tsv")Table of top results.
head(res,50) %>%
kbl(caption = "Top significant genes with GMEA") %>%
kable_paper("hover", full_width = F)| nprobes | mean | median | p-value(sc) | sig | fdr(sc) | |
|---|---|---|---|---|---|---|
| MAD1L1 | 776 | -0.4256926 | -0.4364163 | 0 | 34.23025 | 0 |
| HDAC4 | 504 | -0.4618768 | -0.5754865 | 0 | 27.54266 | 0 |
| ATP11A | 394 | -0.5100916 | -0.5798780 | 0 | 26.47413 | 0 |
| GNAS | 212 | -0.5959431 | -0.6734883 | 0 | 24.38223 | 0 |
| TRAPPC9 | 416 | -0.4741987 | -0.5443188 | 0 | 23.04434 | 0 |
| PDE4D | 282 | -0.6894219 | -0.8931016 | 0 | 22.50192 | 0 |
| TBCD | 403 | -0.4443596 | -0.4683355 | 0 | 22.44445 | 0 |
| CMIP | 226 | -0.5893563 | -0.6298971 | 0 | 22.42396 | 0 |
| FNDC3B | 148 | -0.7673632 | -0.9453880 | 0 | 20.45488 | 0 |
| NCOR2 | 296 | -0.4978154 | -0.5559144 | 0 | 20.00876 | 0 |
| PTPRN2 | 1399 | -0.2586132 | -0.3012425 | 0 | 20.00749 | 0 |
| FILIP1L | 112 | -0.8824234 | -1.0614623 | 0 | 19.44360 | 0 |
| MGMT | 203 | -0.5857988 | -0.6627790 | 0 | 19.41758 | 0 |
| CUX1 | 316 | -0.5078786 | -0.5618714 | 0 | 19.26698 | 0 |
| LRBA | 126 | -0.9035878 | -1.1232366 | 0 | 18.91055 | 0 |
| SMYD3 | 215 | -0.6163169 | -0.7822458 | 0 | 18.67922 | 0 |
| MCPH1 | 123 | -0.9542774 | -0.9878345 | 0 | 18.37229 | 0 |
| ARHGAP15 | 123 | -0.8761692 | -1.0956881 | 0 | 18.02107 | 0 |
| C7orf50 | 350 | -0.4623870 | -0.5656618 | 0 | 17.82372 | 0 |
| MIR1273H | 88 | -0.9863959 | -1.0719828 | 0 | 17.18627 | 0 |
| TCF12 | 103 | -0.9138222 | -1.0737424 | 0 | 16.97517 | 0 |
| ANK3 | 159 | -0.7989265 | -1.0181212 | 0 | 16.32550 | 0 |
| C6orf10 | 90 | -0.9727345 | -1.0712590 | 0 | 16.26369 | 0 |
| C15orf41 | 67 | -0.9433907 | -1.0536925 | 0 | 16.20077 | 0 |
| SND1 | 178 | -0.6434235 | -0.7809250 | 0 | 16.19118 | 0 |
| AGAP1 | 415 | -0.4260793 | -0.4793557 | 0 | 16.18358 | 0 |
| DIP2C | 578 | -0.3304536 | -0.4264368 | 0 | 15.61955 | 0 |
| SLC12A7 | 176 | -0.5535554 | -0.5598834 | 0 | 15.56060 | 0 |
| RPTOR | 540 | -0.3430784 | -0.4420126 | 0 | 15.25391 | 0 |
| TRIO | 244 | -0.5383156 | -0.6808292 | 0 | 15.19742 | 0 |
| EXOC4 | 135 | -0.7767376 | -1.0130211 | 0 | 14.98450 | 0 |
| KCNQ1 | 376 | -0.3809910 | -0.4578382 | 0 | 14.92425 | 0 |
| PRDM16 | 655 | -0.2962127 | -0.3785919 | 0 | 14.79446 | 0 |
| ARHGAP26 | 191 | -0.6176872 | -0.7371087 | 0 | 14.62190 | 0 |
| CMSS1 | 76 | -0.9864276 | -1.1832835 | 0 | 14.43242 | 0 |
| JARID2 | 188 | -0.5602058 | -0.7118861 | 0 | 13.93999 | 0 |
| INPP5A | 418 | -0.3737891 | -0.4026773 | 0 | 13.73301 | 0 |
| CCNY | 111 | -0.6935753 | -0.7901824 | 0 | 13.63168 | 0 |
| BANP | 168 | -0.5405059 | -0.5550284 | 0 | 13.49174 | 0 |
| WDR7 | 66 | -1.0531256 | -1.2199720 | 0 | 13.33680 | 0 |
| NSMCE2 | 118 | -0.7633090 | -0.9206554 | 0 | 13.31774 | 0 |
| FOXK1 | 193 | -0.5129186 | -0.5258607 | 0 | 13.12964 | 0 |
| SYNE2 | 141 | -0.6840765 | -0.7923164 | 0 | 13.06841 | 0 |
| TFDP1 | 96 | -0.5653363 | -0.5984338 | 0 | 13.02035 | 0 |
| SDK1 | 403 | -0.4015345 | -0.4242803 | 0 | 12.96274 | 0 |
| EIF4G3 | 85 | -0.8015040 | -0.9482999 | 0 | 12.84626 | 0 |
| XYLT1 | 144 | -0.6067358 | -0.8054171 | 0 | 12.79535 | 0 |
| SVIL | 140 | -0.6968659 | -0.8878289 | 0 | 12.75763 | 0 |
| IQSEC1 | 165 | -0.5032054 | -0.5493782 | 0 | 12.73063 | 0 |
| TMCC1 | 69 | -0.8443072 | -1.0097807 | 0 | 12.71731 | 0 |
Volcano plots
# volcano selfcont
sig <- subset(res,`fdr(sc)` < 0.05)
plot(res$median , -log10(res$`p-value(sc)`) ,
xlab="effect size (mean t-stat)", ylab="-log10(p-value)",
pch=19, cex=0.5, col="gray",main="self contained test")
grid()
points(sig$median , -log10(sig$`p-value(sc)`) ,
pch=19, cex=0.5, col="red")
# nprobes plot
plot(res$nprobes,res$sig,log="x",ylim=c(0,50),pch=19,cex=0.6)
points(sig$nprobes,sig$sig,col="red",pch=19,cex=0.62)
MIN = min(sig$nprobes)
LEFT = nrow(subset(res,nprobes<MIN))
RIGHT = nrow(subset(res,nprobes>MIN))
SIG = nrow(sig)
TOT = nrow(res)
HEADER <- paste(TOT, "genes in total.", SIG, "with FDR<0.05.",
RIGHT, "well covered and", LEFT, "poorly covered")
mtext(HEADER)
abline(v=MIN,lty=2)
Boxplots
Boxplots smallest pvalue.
par(mfrow=c(1,2))
n=50
# self contained
gs <- head(rownames(res),50)
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="smallest p-val(selfcont)",cex.axis=0.6,
xlab="t-statistic")
grid()
n=50
# effect size median
sig <- subset(res,`fdr(sc)` < 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))Mitch
dmscore <- data.frame( res$median * res$sig)
rownames(dmscore) <- rownames(res)
colnames(dmscore) <- "metric"
hist(dmscore$metric)
if ( ! file.exists("ReactomePathways.gmt") ) {
download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip",
destfile="ReactomePathways.gmt.zip")
unzip("ReactomePathways.gmt.zip")
}
file.info("ReactomePathways.gmt")## size isdir mode mtime ctime
## ReactomePathways.gmt 903549 FALSE 664 2022-09-26 16:10:50 2022-09-26 16:10:50
## atime uid gid uname grname
## ReactomePathways.gmt 2022-10-12 14:35:21 1000 1000 mdz mdzgenesets <- gmt_import("ReactomePathways.gmt")
length(genesets)## [1] 2584mres <- mitch_calc(x=dmscore, genesets=genesets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(mres$enrichment_result,20) %>%
kbl(caption = "Top enriched gene sets with GMEA-Mitch") %>%
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 47 | Activation of RAC1 | 11 | 0.0000294 | -0.7273929 | 0.0003722 |
| 380 | Eicosanoids | 12 | 0.0001456 | 0.6331610 | 0.0013890 |
| 1105 | Regulation of TP53 Activity through Association with Co-factors | 14 | 0.0001453 | -0.5862779 | 0.0013890 |
| 87 | Apoptotic cleavage of cell adhesion proteins | 11 | 0.0009672 | -0.5745138 | 0.0068905 |
| 1163 | SARS-CoV-1 modulates host translation machinery | 34 | 0.0000000 | 0.5479364 | 0.0000011 |
| 695 | MET activates PTK2 signaling | 18 | 0.0001205 | -0.5233640 | 0.0012100 |
| 1490 | VEGFR2 mediated vascular permeability | 26 | 0.0000042 | -0.5213855 | 0.0000906 |
| 171 | Carnitine metabolism | 11 | 0.0031804 | -0.5135615 | 0.0177458 |
| 457 | Free fatty acids regulate insulin secretion | 10 | 0.0081269 | -0.4832990 | 0.0377304 |
| 668 | Laminin interactions | 23 | 0.0000752 | -0.4767933 | 0.0008371 |
| 895 | PTK6 Regulates RHO GTPases, RAS GTPase and MAP kinases | 14 | 0.0022587 | -0.4713258 | 0.0136403 |
| 436 | Fatty acids | 15 | 0.0023283 | 0.4539937 | 0.0139864 |
| 903 | Peptide chain elongation | 84 | 0.0000000 | 0.4531341 | 0.0000000 |
| 1098 | Regulation of RUNX1 Expression and Activity | 25 | 0.0001133 | -0.4459065 | 0.0011445 |
| 397 | Eukaryotic Translation Termination | 87 | 0.0000000 | 0.4445115 | 0.0000000 |
| 1397 | Termination of O-glycan biosynthesis | 25 | 0.0001239 | -0.4433569 | 0.0012362 |
| 395 | Eukaryotic Translation Elongation | 88 | 0.0000000 | 0.4414337 | 0.0000000 |
| 697 | MET activates RAS signaling | 11 | 0.0112755 | -0.4411960 | 0.0489659 |
| 1017 | RHOJ GTPase cycle | 51 | 0.0000001 | -0.4389186 | 0.0000020 |
| 1019 | RHOU GTPase cycle | 37 | 0.0000044 | -0.4362632 | 0.0000934 |
mitch_report(mres,outfile="mitch_gmea_wg.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpNdEwHt/mitch_gmea_wg.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|.. | 3%
## inline R code fragments
##
##
|
|.... | 6%
## label: checklibraries (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...... | 9%
## ordinary text without R code
##
##
|
|........ | 12%
## label: peek (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.......... | 15%
## ordinary text without R code
##
##
|
|............ | 18%
## label: metrics (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|................ | 24%
## label: scatterplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ message : logi FALSE
## $ warning : logi FALSE##
|
|................... | 26%
## ordinary text without R code
##
##
|
|..................... | 29%
## label: contourplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ warning : logi FALSE
## $ message : logi FALSE## Contour plot does not apply to unidimensional analysis.##
|
|....................... | 32%
## ordinary text without R code
##
##
|
|......................... | 35%
## label: input_geneset_metrics1 (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|........................... | 38%
## ordinary text without R code
##
##
|
|............................. | 41%
## label: input_geneset_metrics2 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"##
|
|............................... | 44%
## ordinary text without R code
##
##
|
|................................. | 47%
## label: input_geneset_metrics3 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ message : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7##
|
|................................... | 50%
## ordinary text without R code
##
##
|
|..................................... | 53%
## label: echart1d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|....................................... | 56%
## label: echart2d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|......................................... | 59%
## ordinary text without R code
##
##
|
|........................................... | 62%
## label: heatmap (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 10
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|............................................. | 65%
## ordinary text without R code
##
##
|
|............................................... | 68%
## label: effectsize (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|................................................. | 71%
## ordinary text without R code
##
##
|
|................................................... | 74%
## label: results_table (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|...................................................... | 76%
## ordinary text without R code
##
##
|
|........................................................ | 79%
## label: results_table_complete (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|.......................................................... | 82%
## ordinary text without R code
##
##
|
|............................................................ | 85%
## label: detailed_geneset_reports1d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE##
|
|.............................................................. | 88%
## ordinary text without R code
##
##
|
|................................................................ | 91%
## label: detailed_geneset_reports2d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 5
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE
##
##
|
|.................................................................. | 94%
## ordinary text without R code
##
##
|
|.................................................................... | 97%
## label: session_info (with options)
## List of 3
## $ include: logi TRUE
## $ echo : logi TRUE
## $ results: chr "markup"
##
##
|
|......................................................................| 100%
## ordinary text without R code## 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/RtmpNdEwHt/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/RtmpNdEwHt/rmarkdown-str1ab483ec39dca.html##
## Output created: /tmp/RtmpNdEwHt/mitch_report.html## [1] TRUE1-sided T gene set test
Consider using 1-sided t-test if mitch is biased.
m1 <- dmscore
gse_res <- mclapply( 1:length(genesets), function(i) {
gs <- genesets[i]
name <- names(gs)
n_members <- length(which(rownames(m1) %in% gs[[1]]))
if ( n_members > 4 ) {
tstats <- m1[which(rownames(m1) %in% gs[[1]]),]
myn <- length(tstats)
mymean <- mean(tstats)
mymedian <- median(tstats)
wt <- t.test(tstats)
res <- c(name,myn,mymean,mymedian,wt$p.value)
}
} , mc.cores = CORES)
gse_res_df <- do.call(rbind, gse_res)
rownames(gse_res_df) <- gse_res_df[,1]
gse_res_df <- gse_res_df[,-1]
colnames(gse_res_df) <- c("n_genes","t_mean","t_median","p.value")
tmp <- apply(gse_res_df,2,as.numeric)
rownames(tmp) <- rownames(gse_res_df)
gse_res_df <- tmp
gse_res_df <- as.data.frame(gse_res_df)
gse_res_df <- gse_res_df[order(gse_res_df$p.value),]
gse_res_df$logp <- -log10(gse_res_df$p.value )
gse_res_df$fdr <- p.adjust(gse_res_df$p.value,method="fdr")
head(gse_res_df)## n_genes t_mean t_median p.value
## Signal Transduction 2409 -1.1654700 -0.4533366 2.870333e-167
## Metabolism 1963 -0.9290406 -0.3419812 2.985521e-129
## Immune System 1870 -0.9220815 -0.3634758 2.886018e-125
## Metabolism of proteins 1767 -0.8825232 -0.3284561 3.051954e-105
## Disease 1567 -1.0153933 -0.3618446 5.479151e-103
## Gene expression (Transcription) 1420 -0.8983687 -0.3354673 1.824764e-88
## logp fdr
## Signal Transduction 166.54207 5.545484e-164
## Metabolism 128.52498 2.884013e-126
## Immune System 124.53970 1.858595e-122
## Metabolism of proteins 104.51542 1.474094e-102
## Disease 102.26129 2.117144e-100
## Gene expression (Transcription) 87.73879 5.875740e-86head(gse_res_df[order(-abs(gse_res_df$t_median)),],10)## n_genes t_mean
## MET interacts with TNS proteins 5 -2.390064
## MET activates PTPN11 5 -2.598720
## Activation of RAC1 11 -2.815608
## Role of ABL in ROBO-SLIT signaling 8 -3.578171
## Loss of MECP2 binding ability to the NCoR/SMRT complex 5 -3.540679
## Maturation of protein 3a_9683673 9 -2.173795
## Maturation of protein 3a_9694719 9 -2.173795
## CREB phosphorylation 6 -1.998658
## RUNX1 regulates estrogen receptor mediated transcription 6 -2.379295
## ALK mutants bind TKIs 12 -2.184915
## t_median p.value
## MET interacts with TNS proteins -2.959258 0.044173058
## MET activates PTPN11 -2.959258 0.055318710
## Activation of RAC1 -2.876952 0.001299711
## Role of ABL in ROBO-SLIT signaling -2.743888 0.019862207
## Loss of MECP2 binding ability to the NCoR/SMRT complex -2.733640 0.149391741
## Maturation of protein 3a_9683673 -2.717529 0.001137865
## Maturation of protein 3a_9694719 -2.717529 0.001137865
## CREB phosphorylation -2.579213 0.022122160
## RUNX1 regulates estrogen receptor mediated transcription -2.254571 0.009024103
## ALK mutants bind TKIs -2.217244 0.004922192
## logp fdr
## MET interacts with TNS proteins 1.3548425 0.057937779
## MET activates PTPN11 1.2571280 0.069264905
## Activation of RAC1 2.8861532 0.003402496
## Role of ABL in ROBO-SLIT signaling 1.7019725 0.031147552
## Loss of MECP2 binding ability to the NCoR/SMRT complex 0.8256734 0.159695606
## Maturation of protein 3a_9683673 2.9439094 0.003015229
## Maturation of protein 3a_9694719 2.9439094 0.003015229
## CREB phosphorylation 1.6551725 0.033866889
## RUNX1 regulates estrogen receptor mediated transcription 2.0445960 0.016572782
## ALK mutants bind TKIs 2.3078414 0.010041895head(gse_res_df,30) %>%
kbl(caption = "Top significant genesets") %>%
kable_paper("hover", full_width = F)| n_genes | t_mean | t_median | p.value | logp | fdr | |
|---|---|---|---|---|---|---|
| Signal Transduction | 2409 | -1.1654700 | -0.4533366 | 0 | 166.54207 | 0 |
| Metabolism | 1963 | -0.9290406 | -0.3419812 | 0 | 128.52498 | 0 |
| Immune System | 1870 | -0.9220815 | -0.3634758 | 0 | 124.53970 | 0 |
| Metabolism of proteins | 1767 | -0.8825232 | -0.3284561 | 0 | 104.51542 | 0 |
| Disease | 1567 | -1.0153933 | -0.3618446 | 0 | 102.26129 | 0 |
| Gene expression (Transcription) | 1420 | -0.8983687 | -0.3354673 | 0 | 87.73879 | 0 |
| Post-translational protein modification | 1289 | -0.9780871 | -0.3695283 | 0 | 86.09808 | 0 |
| RNA Polymerase II Transcription | 1284 | -0.8820708 | -0.3230669 | 0 | 77.33624 | 0 |
| Generic Transcription Pathway | 1163 | -0.8989980 | -0.3229938 | 0 | 70.11017 | 0 |
| Innate Immune System | 999 | -0.9531670 | -0.3610626 | 0 | 65.14931 | 0 |
| Developmental Biology | 1027 | -1.0429386 | -0.3563097 | 0 | 61.13019 | 0 |
| Infectious disease | 830 | -0.9049617 | -0.2819514 | 0 | 51.38012 | 0 |
| Metabolism of lipids | 696 | -1.1134649 | -0.4803375 | 0 | 51.34453 | 0 |
| Transport of small molecules | 689 | -1.0558477 | -0.4259990 | 0 | 50.62418 | 0 |
| Sensory Perception | 570 | -1.0155699 | -0.4639070 | 0 | 49.55422 | 0 |
| Adaptive Immune System | 709 | -0.9474317 | -0.3625209 | 0 | 49.20157 | 0 |
| Signaling by GPCR | 670 | -1.1308502 | -0.4309666 | 0 | 46.77641 | 0 |
| Signaling by Rho GTPases, Miro GTPases and RHOBTB3 | 603 | -1.3656777 | -0.5793046 | 0 | 46.02353 | 0 |
| Vesicle-mediated transport | 623 | -1.1925647 | -0.4805157 | 0 | 45.73560 | 0 |
| Signaling by Rho GTPases | 589 | -1.3812741 | -0.5848349 | 0 | 45.21005 | 0 |
| GPCR downstream signalling | 603 | -1.1867695 | -0.4491977 | 0 | 43.14734 | 0 |
| Membrane Trafficking | 585 | -1.1842461 | -0.4769073 | 0 | 42.59354 | 0 |
| Hemostasis | 555 | -1.2662532 | -0.4868853 | 0 | 41.98358 | 0 |
| Cellular responses to stimuli | 696 | -0.8450616 | -0.2091962 | 0 | 38.76568 | 0 |
| Cytokine Signaling in Immune system | 670 | -0.8361616 | -0.2421952 | 0 | 38.50294 | 0 |
| Cellular responses to stress | 682 | -0.8562656 | -0.2091962 | 0 | 38.30069 | 0 |
| Signaling by Receptor Tyrosine Kinases | 493 | -1.2327809 | -0.5046052 | 0 | 37.18123 | 0 |
| RHO GTPase cycle | 420 | -1.5410841 | -0.6616470 | 0 | 36.08182 | 0 |
| Neutrophil degranulation | 458 | -1.0214350 | -0.4515763 | 0 | 32.50754 | 0 |
| Nervous system development | 527 | -1.0607187 | -0.2914606 | 0 | 32.39422 | 0 |
head(gse_res_df[order(-abs(gse_res_df$t_median)),],30) %>%
kbl(caption = "Top effect size genesets") %>%
kable_paper("hover", full_width = F)| n_genes | t_mean | t_median | p.value | logp | fdr | |
|---|---|---|---|---|---|---|
| MET interacts with TNS proteins | 5 | -2.390064 | -2.959258 | 0.0441731 | 1.3548425 | 0.0579378 |
| MET activates PTPN11 | 5 | -2.598720 | -2.959258 | 0.0553187 | 1.2571280 | 0.0692649 |
| Activation of RAC1 | 11 | -2.815608 | -2.876952 | 0.0012997 | 2.8861532 | 0.0034025 |
| Role of ABL in ROBO-SLIT signaling | 8 | -3.578171 | -2.743888 | 0.0198622 | 1.7019725 | 0.0311476 |
| Loss of MECP2 binding ability to the NCoR/SMRT complex | 5 | -3.540679 | -2.733640 | 0.1493917 | 0.8256734 | 0.1596956 |
| Maturation of protein 3a_9683673 | 9 | -2.173795 | -2.717529 | 0.0011379 | 2.9439094 | 0.0030152 |
| Maturation of protein 3a_9694719 | 9 | -2.173795 | -2.717529 | 0.0011379 | 2.9439094 | 0.0030152 |
| CREB phosphorylation | 6 | -1.998658 | -2.579213 | 0.0221222 | 1.6551725 | 0.0338669 |
| RUNX1 regulates estrogen receptor mediated transcription | 6 | -2.379295 | -2.254571 | 0.0090241 | 2.0445960 | 0.0165728 |
| ALK mutants bind TKIs | 12 | -2.184915 | -2.217244 | 0.0049222 | 2.3078414 | 0.0100419 |
| LRR FLII-interacting protein 1 (LRRFIP1) activates type I IFN production | 5 | -2.348674 | -2.184857 | 0.0963299 | 1.0162388 | 0.1094761 |
| Post-transcriptional silencing by small RNAs | 7 | -1.994432 | -2.118831 | 0.0194235 | 1.7116722 | 0.0306142 |
| Loss of Function of TGFBR1 in Cancer | 7 | -1.811190 | -2.103727 | 0.0092307 | 2.0347673 | 0.0169152 |
| VEGFR2 mediated vascular permeability | 26 | -2.339323 | -2.069563 | 0.0000200 | 4.6988419 | 0.0000976 |
| RHO GTPases activate IQGAPs | 11 | -2.103286 | -2.052799 | 0.0067380 | 2.1714685 | 0.0131228 |
| DNA methylation | 6 | -2.056379 | -2.007306 | 0.0356601 | 1.4478172 | 0.0489661 |
| Regulation of TP53 Activity through Association with Co-factors | 14 | -2.739041 | -1.986231 | 0.0017039 | 2.7685659 | 0.0042421 |
| RUNX2 regulates genes involved in cell migration | 8 | -2.928541 | -1.977132 | 0.0169978 | 1.7696081 | 0.0275732 |
| RUNX1 regulates transcription of genes involved in WNT signaling | 5 | -2.483370 | -1.918885 | 0.0157033 | 1.8040091 | 0.0258863 |
| RUNX1 regulates expression of components of tight junctions | 5 | -1.821143 | -1.918885 | 0.0230733 | 1.6368905 | 0.0348737 |
| AKT phosphorylates targets in the nucleus | 9 | -2.200881 | -1.872909 | 0.0264791 | 1.5770960 | 0.0386680 |
| SMAD2/3 Phosphorylation Motif Mutants in Cancer | 6 | -1.762434 | -1.794554 | 0.0262222 | 1.5813311 | 0.0383797 |
| TGFBR1 KD Mutants in Cancer | 6 | -1.762434 | -1.794554 | 0.0262222 | 1.5813311 | 0.0383797 |
| RUNX2 regulates chondrocyte maturation | 5 | -4.191764 | -1.765786 | 0.2263833 | 0.6451557 | 0.2331410 |
| Defective ST3GAL3 causes MCT12 and EIEE15 | 8 | -2.302302 | -1.761819 | 0.0421586 | 1.3751139 | 0.0558645 |
| Signaling by TGF-beta Receptor Complex in Cancer | 8 | -1.648131 | -1.746331 | 0.0077290 | 2.1118752 | 0.0146829 |
| MET activates PI3K/AKT signaling | 6 | -2.203227 | -1.742876 | 0.0555608 | 1.2552316 | 0.0694779 |
| Inactivation of CDC42 and RAC1 | 8 | -2.347666 | -1.681234 | 0.0448249 | 1.3484803 | 0.0585543 |
| MET activates PTK2 signaling | 18 | -2.304248 | -1.630240 | 0.0002453 | 3.6103670 | 0.0008600 |
| RORA activates gene expression | 17 | -2.283110 | -1.610898 | 0.0029097 | 2.5361556 | 0.0065672 |
sig <- subset(gse_res_df,fdr<0.05)
nsig <- nrow(sig)
# volcano
plot(gse_res_df$t_median , gse_res_df$logp ,
xlab="effect size (mean t-stat)", ylab="-log10(p-value)",
pch=19, cex=0.5, col="gray")
points(sig$t_median , sig$logp ,
pch=19, cex=0.5, col="red")
abline(h=5,lty=2)
# boxplot most significant
n=50
sig2 <- head(gse_res_df,n)
sig2 <- sig2[order(-sig2$t_median),]
gs <- rownames(sig2)[1:n]
tstats <- lapply(1:n, function(i) {
g <- rownames(sig2)[i]
genes <- genesets[[which(names(genesets) == g)]]
m1[which(rownames(m1) %in% genes),"metric"]
} )
names(tstats) <- gs
# boxplot
par(mar=c(5,20,3,1))
boxplot(tstats,horizontal=TRUE,las=1,
main="most significant gene sets",cex.axis=0.6,
xlab="mean t-statistic")
grid()
# boxplot biggest effect size
n=50
sig2 <- tail(sig[order(abs(sig$t_median)),],n)
gs <- rownames(sig2)[1:n]
tstats <- lapply(1:n, function(i) {
g <- rownames(sig2)[i]
genes <- genesets[[which(names(genesets) == g)]]
m1[which(rownames(m1) %in% genes),"metric"]
} )
names(tstats) <- gs
# boxplot
boxplot(tstats,horizontal=TRUE,las=1,
main="largest effect size",cex.axis=0.6,
xlab="mean t-statistic")
grid()
par(mfrow=c(1,1))Promoter analysis only
dm <- read.csv("ASD_blood_top_dmps_onADOS_withintw_Nov27_limma.csv")
rownames(dm) <- dm$Name
dm <- dm[,9:ncol(dm)]
dm <- merge(dm,myann,by=0)
rownames(dm) <- dm[,1]
dm[,1] = NULL
dm <- dm[grep("Promo",dm$Regulatory_Feature_Group),]
dim(dm)## [1] 110628 8dm <- dm[order(dm$P.Value),]
head(dm,50) %>%
kbl(caption = "Top significant genes with limma") %>%
kable_paper("hover", full_width = F)| logFC | AveExpr | t | P.Value | adj.P.Val | B | UCSC_RefGene_Name | Regulatory_Feature_Group | |
|---|---|---|---|---|---|---|---|---|
| cg07655025 | -0.0854949 | -4.336160 | -6.706650 | 0.0000228 | 0.8764716 | 3.0481274 | TNPO2 | Promoter_Associated |
| cg22492966 | 0.0655305 | -4.371072 | 5.329706 | 0.0001853 | 0.8764716 | 1.0690876 | JMJD1C;JMJD1C;JMJD1C | Promoter_Associated |
| cg04835122 | 0.0720488 | -4.723115 | 4.687265 | 0.0005392 | 0.8764716 | 0.0438751 | STAMBP;STAMBP;STAMBP | Promoter_Associated |
| cg01972157 | 0.0675974 | -4.682192 | 4.612851 | 0.0006124 | 0.8764716 | -0.0787834 | GSPT1;GSPT1;GSPT1;GSPT1 | Promoter_Associated |
| cg01287672 | 0.0675047 | -4.451346 | 4.608839 | 0.0006167 | 0.8764716 | -0.0854170 | UTP6 | Promoter_Associated |
| cg26034408 | 0.0879174 | -5.323799 | 4.585468 | 0.0006419 | 0.8764716 | -0.1241105 | MTFR1;MTFR1 | Promoter_Associated |
| cg17699947 | 0.0767633 | -5.150053 | 4.552240 | 0.0006797 | 0.8764716 | -0.1792516 | KBTBD8 | Promoter_Associated |
| cg20885078 | 0.0573758 | -5.187835 | 4.547303 | 0.0006855 | 0.8764716 | -0.1874573 | MIR548N | Promoter_Associated |
| cg02713138 | 0.0610751 | -3.543557 | 4.545234 | 0.0006880 | 0.8764716 | -0.1908963 | FAM55C;FAM55C;NFKBIZ | Promoter_Associated |
| cg23146534 | 0.0616873 | -3.744097 | 4.429817 | 0.0008402 | 0.8764716 | -0.3836713 | C7orf38 | Promoter_Associated |
| cg01494348 | 0.0634881 | -3.309309 | 4.427392 | 0.0008438 | 0.8764716 | -0.3877411 | NAPRT1 | Promoter_Associated |
| cg09367367 | -0.0655440 | -4.792395 | -4.407008 | 0.0008743 | 0.8764716 | -0.4219733 | RPL6;RPL6 | Promoter_Associated_Cell_type_specific |
| cg03128972 | -0.0633462 | -2.729192 | -4.404630 | 0.0008779 | 0.8764716 | -0.4259702 | EPN2;EPN2;EPN2 | Promoter_Associated_Cell_type_specific |
| cg16101257 | 0.0556665 | -6.090360 | 4.371233 | 0.0009306 | 0.8764716 | -0.4821788 | ZNF517 | Promoter_Associated |
| cg25173405 | 0.0582739 | -3.184171 | 4.362890 | 0.0009443 | 0.8764716 | -0.4962423 | C17orf57 | Promoter_Associated |
| cg02834886 | 0.0491714 | -4.603489 | 4.361484 | 0.0009466 | 0.8764716 | -0.4986137 | SERINC1;PKIB;PKIB;PKIB | Promoter_Associated |
| cg01474603 | -0.0718588 | -4.032645 | -4.360659 | 0.0009479 | 0.8764716 | -0.5000042 | UHRF1BP1L;UHRF1BP1L | Promoter_Associated |
| cg11806890 | -0.0476313 | -4.815429 | -4.254866 | 0.0011413 | 0.8764716 | -0.6791000 | MRPL4;MRPL4;MRPL4 | Promoter_Associated |
| cg18395777 | -0.0819403 | -5.323954 | -4.245810 | 0.0011596 | 0.8764716 | -0.6944920 | ACADM;ACADM | Promoter_Associated |
| cg02429684 | 0.0562522 | -4.698572 | 4.223465 | 0.0012062 | 0.8764716 | -0.7325127 | NAPA | Promoter_Associated |
| cg04304927 | 0.0474021 | -4.420114 | 4.219361 | 0.0012150 | 0.8764716 | -0.7395025 | SLC39A8;SLC39A8;SLC39A8;SLC39A8 | Promoter_Associated |
| cg13018448 | -0.0595578 | -3.046104 | -4.205468 | 0.0012452 | 0.8764716 | -0.7631772 | CUTA;CUTA;CUTA;CUTA;CUTA | Promoter_Associated |
| cg24019599 | 0.0686707 | -2.179191 | 4.198263 | 0.0012611 | 0.8764716 | -0.7754633 | Promoter_Associated | |
| cg09911982 | -0.0638613 | -5.265260 | -4.176913 | 0.0013097 | 0.8764716 | -0.8119062 | Promoter_Associated | |
| cg23594990 | 0.0640998 | -5.440064 | 4.168794 | 0.0013286 | 0.8764716 | -0.8257776 | KIAA1731 | Promoter_Associated |
| cg20290983 | 0.1044379 | -4.494748 | 4.131288 | 0.0014200 | 0.8764716 | -0.8899512 | MRPS18A | Promoter_Associated |
| cg16318349 | -0.0613010 | 2.531760 | -4.115089 | 0.0014614 | 0.8764716 | -0.9177145 | PBXIP1 | Promoter_Associated |
| cg02436796 | -0.0576997 | -2.867308 | -4.103308 | 0.0014924 | 0.8764716 | -0.9379246 | ST7;ST7;ST7OT4;ST7OT1 | Promoter_Associated |
| cg13969672 | 0.0547328 | -6.357855 | 4.095700 | 0.0015127 | 0.8764716 | -0.9509829 | SNRPC;SNRPC;SNRPC | Promoter_Associated |
| cg08505684 | -0.0558090 | -4.449155 | -4.072379 | 0.0015768 | 0.8764716 | -0.9910498 | SAR1B;SAR1B | Promoter_Associated |
| cg06604625 | -0.0547269 | -3.700016 | -4.063440 | 0.0016021 | 0.8764716 | -1.0064225 | CCDC46;CCDC46 | Promoter_Associated |
| cg04662451 | 0.0467976 | -5.045994 | 4.061632 | 0.0016073 | 0.8764716 | -1.0095338 | RAD51;RAD51;RAD51;RAD51 | Promoter_Associated |
| cg10127679 | -0.0479377 | -5.023762 | -4.058561 | 0.0016161 | 0.8764716 | -1.0148161 | UBXN1 | Promoter_Associated |
| cg03355757 | -0.0558783 | -4.437940 | -4.038967 | 0.0016736 | 0.8764716 | -1.0485517 | ZMYM6NB | Promoter_Associated |
| cg15424910 | -0.1416113 | -2.768785 | -4.029890 | 0.0017010 | 0.8764716 | -1.0641911 | THOC5;THOC5;THOC5;THOC5 | Promoter_Associated |
| cg15992940 | 0.0578560 | -4.190814 | 4.027011 | 0.0017098 | 0.8764716 | -1.0691547 | SNHG4;MATR3 | Promoter_Associated |
| cg11526113 | -0.1076416 | -5.114586 | -4.001903 | 0.0017883 | 0.8764716 | -1.1124671 | Promoter_Associated | |
| cg09554789 | -0.0817259 | -4.312655 | -3.997035 | 0.0018039 | 0.8764716 | -1.1208724 | PNPLA6;PNPLA6;PNPLA6;PNPLA6 | Promoter_Associated |
| cg04168050 | 0.0750196 | -4.075329 | 3.977709 | 0.0018674 | 0.8764716 | -1.1542621 | BPNT1 | Promoter_Associated |
| cg00236696 | 0.0713155 | -4.689500 | 3.955640 | 0.0019428 | 0.8764716 | -1.1924329 | HCFC2 | Promoter_Associated |
| cg14978691 | 0.0560474 | -4.056448 | 3.952466 | 0.0019539 | 0.8764716 | -1.1979258 | ARPC1B | Promoter_Associated |
| cg13363438 | 0.0565006 | -4.221095 | 3.952033 | 0.0019554 | 0.8764716 | -1.1986760 | CLDN10;CLDN10 | Promoter_Associated |
| cg01538982 | 0.0626260 | -4.996327 | 3.929109 | 0.0020376 | 0.8764716 | -1.2383810 | CCDC51;CCDC51;CCDC72 | Promoter_Associated |
| cg04406808 | -0.0771244 | -4.829050 | -3.920231 | 0.0020704 | 0.8764716 | -1.2537714 | HSD17B12;HSD17B12 | Promoter_Associated |
| cg15448144 | 0.0829955 | -4.645548 | 3.909055 | 0.0021124 | 0.8764716 | -1.2731541 | OTUD6B | Promoter_Associated |
| cg01811895 | 0.0495059 | -4.747199 | 3.902728 | 0.0021366 | 0.8764716 | -1.2841313 | CBLB | Promoter_Associated_Cell_type_specific |
| cg00320749 | 0.0943612 | -3.962859 | 3.884037 | 0.0022097 | 0.8764716 | -1.3165840 | SNORD43;RPL3;RPL3 | Promoter_Associated |
| cg24147849 | -0.0887884 | -4.040269 | -3.883108 | 0.0022134 | 0.8764716 | -1.3181972 | NIT2 | Promoter_Associated |
| cg22259932 | 0.0667546 | -5.804034 | 3.881047 | 0.0022216 | 0.8764716 | -1.3217771 | UCRC;UCRC;ZMAT5;ZMAT5;ZMAT5;ZMAT5 | Promoter_Associated |
| cg26301297 | 0.0610574 | -3.965274 | 3.873738 | 0.0022511 | 0.8764716 | -1.3344761 | NSMCE4A;NSMCE4A | Promoter_Associated |
dm <- dm[,c("UCSC_RefGene_Name","t")]
# histogram of t values
hist(dm$t,breaks=seq(from=-7,to=6,by=1))
# set cores to used for parallel execution
tic()
gme_res_promoter <- calc_sc(dm)
time2 <- toc() #4core 30.3 ## 20.256 sec elapsedtime2## $tic
## elapsed
## 6049.461
##
## $toc
## elapsed
## 6069.717
##
## $msg
## logical(0)
##
## $callback_msg
## [1] "20.256 sec elapsed"df <- gme_res_promoter[[1]]
res <- gme_res_promoter[[2]]
write.table(res ,file="gmea_promo.tsv")Table of top results.
head(res,50) %>%
kbl(caption = "Top significant genes with GMEA") %>%
kable_paper("hover", full_width = F)| nprobes | mean | median | p-value(sc) | sig | fdr(sc) | |
|---|---|---|---|---|---|---|
| KDM2B | 22 | -0.7242432 | -0.6992326 | 0.0000039 | 5.403606 | 0.0473620 |
| MGMT | 18 | -0.8956944 | -0.7985600 | 0.0000052 | 5.288002 | 0.0618014 |
| AXIN1 | 12 | -0.9561582 | -1.0506627 | 0.0000150 | 4.823293 | 0.1801654 |
| ATP10D | 10 | 0.9510437 | 1.0408392 | 0.0000174 | 4.759895 | 0.2084647 |
| KCNQ1OT1 | 12 | -0.5537933 | -0.6313157 | 0.0000217 | 4.662628 | 0.2607735 |
| KCNQ1 | 12 | -0.5537933 | -0.6313157 | 0.0000217 | 4.662628 | 0.2607735 |
| SLC16A3 | 27 | -0.6203555 | -0.7571377 | 0.0000250 | 4.602548 | 0.2994136 |
| SNURF | 10 | -0.7622984 | -0.8388893 | 0.0000256 | 4.592486 | 0.3064061 |
| SNRPN | 10 | -0.7622984 | -0.8388893 | 0.0000256 | 4.592486 | 0.3064061 |
| CTSZ | 11 | -1.1381244 | -1.2008368 | 0.0000435 | 4.361868 | 0.5210053 |
| UFC1 | 12 | 0.9633917 | 1.0438509 | 0.0000561 | 4.251310 | 0.6719928 |
| FBRSL1 | 24 | -0.6404954 | -0.7479343 | 0.0000640 | 4.193761 | 0.7671437 |
| GAS8 | 11 | -0.7776392 | -0.8257174 | 0.0000672 | 4.172664 | 0.8052621 |
| RAB43 | 19 | 0.7528307 | 0.8482663 | 0.0000757 | 4.120914 | 0.9070916 |
| NAP1L5 | 13 | -0.9015721 | -0.9706174 | 0.0000797 | 4.098494 | 0.9550704 |
| NUP205 | 15 | 0.9053235 | 0.7222074 | 0.0000854 | 4.068605 | 1.0000000 |
| FAM50B | 29 | -0.5053885 | -0.5059175 | 0.0001047 | 3.979919 | 1.0000000 |
| FANCC | 7 | 0.7024691 | 0.7947616 | 0.0001182 | 3.927504 | 1.0000000 |
| CWC22 | 14 | 1.2784392 | 1.3768257 | 0.0001316 | 3.880880 | 1.0000000 |
| JDP2 | 15 | -0.8062562 | -0.8637757 | 0.0001359 | 3.866657 | 1.0000000 |
| FILIP1L | 10 | -0.8726546 | -0.7357412 | 0.0001527 | 3.816168 | 1.0000000 |
| LOC100130522 | 6 | -0.9049565 | -0.8877822 | 0.0001744 | 3.758403 | 1.0000000 |
| PLEKHM3 | 12 | -0.6744385 | -0.7368880 | 0.0001833 | 3.736901 | 1.0000000 |
| SH3PXD2A | 17 | -0.6487840 | -0.6934406 | 0.0001866 | 3.729036 | 1.0000000 |
| ARMC9 | 13 | 0.7003096 | 0.6342171 | 0.0001873 | 3.727445 | 1.0000000 |
| LOC728323 | 4 | -0.6829558 | -0.7029425 | 0.0001891 | 3.723349 | 1.0000000 |
| DALRD3 | 25 | -0.6204661 | -0.6269258 | 0.0001991 | 3.700973 | 1.0000000 |
| KRTCAP3 | 7 | -0.9465435 | -0.8742049 | 0.0002097 | 3.678474 | 1.0000000 |
| RWDD3 | 16 | -0.6943266 | -0.6251181 | 0.0002231 | 3.651469 | 1.0000000 |
| SPI1 | 10 | -1.0012682 | -0.8879562 | 0.0003107 | 3.507708 | 1.0000000 |
| ZEB2 | 11 | 1.0543943 | 1.0933811 | 0.0003175 | 3.498239 | 1.0000000 |
| VTI1B | 13 | -0.6502728 | -0.6256825 | 0.0003302 | 3.481211 | 1.0000000 |
| ZNF227 | 14 | 0.9679130 | 1.1816602 | 0.0003447 | 3.462505 | 1.0000000 |
| C3orf38 | 18 | 0.5139933 | 0.4731383 | 0.0003946 | 3.403871 | 1.0000000 |
| LOC401431 | 22 | -0.8395424 | -0.8354688 | 0.0003968 | 3.401401 | 1.0000000 |
| NDUFA9 | 15 | 0.7133179 | 0.7457509 | 0.0004066 | 3.390781 | 1.0000000 |
| GNB1 | 5 | -1.3460734 | -1.4217668 | 0.0004081 | 3.389201 | 1.0000000 |
| LOC101928069 | 4 | -0.7171417 | -0.7260065 | 0.0004180 | 3.378816 | 1.0000000 |
| TRIB3 | 6 | 0.4092316 | 0.3987227 | 0.0004380 | 3.358493 | 1.0000000 |
| PSIMCT-1 | 8 | -0.6515312 | -0.6752221 | 0.0004677 | 3.330044 | 1.0000000 |
| LCK | 13 | 0.7611715 | 0.8018580 | 0.0004862 | 3.313228 | 1.0000000 |
| KIFC2 | 16 | -0.7720185 | -0.8862975 | 0.0005026 | 3.298794 | 1.0000000 |
| IRF5 | 18 | -0.7433843 | -0.8865655 | 0.0005496 | 3.259966 | 1.0000000 |
| ATP6V0E2 | 22 | -0.8217652 | -0.8354688 | 0.0005534 | 3.256944 | 1.0000000 |
| ANGEL2 | 13 | 0.5984283 | 0.6229189 | 0.0005644 | 3.248414 | 1.0000000 |
| DHX16 | 39 | 0.3925620 | 0.3827674 | 0.0005930 | 3.226956 | 1.0000000 |
| P2RX1 | 8 | -0.9525725 | -1.0223215 | 0.0006342 | 3.197760 | 1.0000000 |
| ARL6IP4 | 17 | -0.8798686 | -0.8238255 | 0.0006682 | 3.175121 | 1.0000000 |
| SORD | 13 | -0.4463849 | -0.5443812 | 0.0006816 | 3.166495 | 1.0000000 |
| ZNF235 | 10 | 0.7900609 | 0.7104605 | 0.0006939 | 3.158718 | 1.0000000 |
Volcano plots
# volcano selfcont
sig <- subset(res,`fdr(sc)` < 0.05)
plot(res$median , -log10(res$`p-value(sc)`) ,
xlab="effect size (mean t-stat)", ylab="-log10(p-value)",
pch=19, cex=0.5, col="gray",main="self contained test")
grid()
points(sig$median , -log10(sig$`p-value(sc)`) ,
pch=19, cex=0.5, col="red")
Boxplots
Boxplots smallest pvalue.
par(mfrow=c(1,2))
n=50
# self contained
gs <- head(rownames(res),50)
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="smallest p-val(selfcont)",cex.axis=0.6,
xlab="t-statistic")
grid()
n=50
# effect size median
sig <- subset(res,`fdr(sc)` < 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))Mitch
dmscore <- data.frame( res$median * res$sig)
rownames(dmscore) <- rownames(res)
colnames(dmscore) <- "metric"
genesets <- gmt_import("ReactomePathways.gmt")
mres <- mitch_calc(x=dmscore, genesets=genesets,priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(mres$enrichment_result,20) %>%
kbl(caption = "Top enriched gene sets with GMEA-Mitch (promoter only)") %>%
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 1058 | Synthesis of PIPs at the early endosome membrane | 12 | 0.0028154 | 0.4980251 | 0.0429693 |
| 906 | SARS-CoV-1 targets host intracellular signalling and regulatory pathways | 15 | 0.0020173 | -0.4604791 | 0.0361536 |
| 1073 | TNF receptor superfamily (TNFSF) members mediating non-canonical NF-kB pathway | 10 | 0.0118452 | 0.4596196 | 0.0975918 |
| 1072 | TICAM1-dependent activation of IRF3/IRF7 | 12 | 0.0089312 | 0.4359354 | 0.0845349 |
| 333 | Formation of tubulin folding intermediates by CCT/TriC | 14 | 0.0049309 | 0.4340082 | 0.0567981 |
| 549 | Metabolism of nitric oxide: NOS3 activation and regulation | 12 | 0.0093294 | -0.4334460 | 0.0850091 |
| 912 | SARS-CoV-2 targets host intracellular signalling and regulatory pathways | 11 | 0.0223410 | -0.3978230 | 0.1498812 |
| 208 | Defective HDR through Homologous Recombination (HRR) due to BRCA1 loss-of-function | 22 | 0.0012476 | 0.3975887 | 0.0267254 |
| 209 | Defective HDR through Homologous Recombination (HRR) due to PALB2 loss of function | 22 | 0.0012476 | 0.3975887 | 0.0267254 |
| 210 | Defective HDR through Homologous Recombination Repair (HRR) due to PALB2 loss of BRCA1 binding function | 22 | 0.0012476 | 0.3975887 | 0.0267254 |
| 211 | Defective HDR through Homologous Recombination Repair (HRR) due to PALB2 loss of BRCA2/RAD51/RAD51C binding function | 22 | 0.0012476 | 0.3975887 | 0.0267254 |
| 437 | Impaired BRCA2 binding to PALB2 | 22 | 0.0012476 | 0.3975887 | 0.0267254 |
| 182 | Cytochrome c-mediated apoptotic response | 11 | 0.0254156 | 0.3892062 | 0.1623267 |
| 689 | PTK6 Regulates RHO GTPases, RAS GTPase and MAP kinases | 10 | 0.0368899 | -0.3811530 | 0.2066175 |
| 730 | Prefoldin mediated transfer of substrate to CCT/TriC | 17 | 0.0097817 | 0.3619619 | 0.0871787 |
| 585 | Molecules associated with elastic fibres | 11 | 0.0380061 | -0.3613077 | 0.2111241 |
| 604 | NOTCH3 Activation and Transmission of Signal to the Nucleus | 15 | 0.0156135 | 0.3606099 | 0.1199002 |
| 1 | A tetrasaccharide linker sequence is required for GAG synthesis | 10 | 0.0487791 | -0.3598865 | 0.2457379 |
| 286 | Elevation of cytosolic Ca2+ levels | 10 | 0.0506498 | -0.3569498 | 0.2505371 |
| 1043 | Stabilization of p53 | 49 | 0.0000170 | 0.3553432 | 0.0069153 |
mitch_report(mres,outfile="gmea_mitch_promo.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpNdEwHt/gmea_mitch_promo.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|.. | 3%
## inline R code fragments
##
##
|
|.... | 6%
## label: checklibraries (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...... | 9%
## ordinary text without R code
##
##
|
|........ | 12%
## label: peek (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.......... | 15%
## ordinary text without R code
##
##
|
|............ | 18%
## label: metrics (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|................ | 24%
## label: scatterplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ message : logi FALSE
## $ warning : logi FALSE##
|
|................... | 26%
## ordinary text without R code
##
##
|
|..................... | 29%
## label: contourplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ warning : logi FALSE
## $ message : logi FALSE## Contour plot does not apply to unidimensional analysis.##
|
|....................... | 32%
## ordinary text without R code
##
##
|
|......................... | 35%
## label: input_geneset_metrics1 (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|........................... | 38%
## ordinary text without R code
##
##
|
|............................. | 41%
## label: input_geneset_metrics2 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"##
|
|............................... | 44%
## ordinary text without R code
##
##
|
|................................. | 47%
## label: input_geneset_metrics3 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ message : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7##
|
|................................... | 50%
## ordinary text without R code
##
##
|
|..................................... | 53%
## label: echart1d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|....................................... | 56%
## label: echart2d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|......................................... | 59%
## ordinary text without R code
##
##
|
|........................................... | 62%
## label: heatmap (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 10
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|............................................. | 65%
## ordinary text without R code
##
##
|
|............................................... | 68%
## label: effectsize (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|................................................. | 71%
## ordinary text without R code
##
##
|
|................................................... | 74%
## label: results_table (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|...................................................... | 76%
## ordinary text without R code
##
##
|
|........................................................ | 79%
## label: results_table_complete (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|.......................................................... | 82%
## ordinary text without R code
##
##
|
|............................................................ | 85%
## label: detailed_geneset_reports1d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE##
|
|.............................................................. | 88%
## ordinary text without R code
##
##
|
|................................................................ | 91%
## label: detailed_geneset_reports2d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 5
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE
##
##
|
|.................................................................. | 94%
## ordinary text without R code
##
##
|
|.................................................................... | 97%
## label: session_info (with options)
## List of 3
## $ include: logi TRUE
## $ echo : logi TRUE
## $ results: chr "markup"
##
##
|
|......................................................................| 100%
## ordinary text without R code## 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/RtmpNdEwHt/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/RtmpNdEwHt/rmarkdown-str1ab482446ae1c.html##
## Output created: /tmp/RtmpNdEwHt/mitch_report.html## [1] TRUE1-sided t gene set test
Consider using 1-sided t-test if mitch is biased.
m1 <- res[,"median",drop=FALSE]
gse_res <- mclapply( 1:length(genesets), function(i) {
gs <- genesets[i]
name <- names(gs)
n_members <- length(which(rownames(m1) %in% gs[[1]]))
if ( n_members > 4 ) {
tstats <- m1[which(rownames(m1) %in% gs[[1]]),]
myn <- length(tstats)
mymean <- mean(tstats)
mymedian <- median(tstats)
wt <- t.test(tstats)
res <- c(name,myn,mymean,mymedian,wt$p.value)
}
} , mc.cores = CORES)
gse_res_df <- do.call(rbind, gse_res)
rownames(gse_res_df) <- gse_res_df[,1]
gse_res_df <- gse_res_df[,-1]
colnames(gse_res_df) <- c("n_genes","t_mean","t_median","p.value")
tmp <- apply(gse_res_df,2,as.numeric)
rownames(tmp) <- rownames(gse_res_df)
gse_res_df <- tmp
gse_res_df <- as.data.frame(gse_res_df)
gse_res_df <- gse_res_df[order(gse_res_df$p.value),]
gse_res_df$logp <- -log10(gse_res_df$p.value )
gse_res_df$fdr <- p.adjust(gse_res_df$p.value,method="fdr")
head(gse_res_df)## n_genes
## Metabolism of RNA 555
## DNA Repair 250
## SLC-mediated transmembrane transport 105
## Neuronal System 139
## SARS-CoV-1 targets host intracellular signalling and regulatory pathways 15
## Induction of Cell-Cell Fusion 5
## t_mean
## Metabolism of RNA 0.07534063
## DNA Repair 0.08829874
## SLC-mediated transmembrane transport -0.14804880
## Neuronal System -0.14329149
## SARS-CoV-1 targets host intracellular signalling and regulatory pathways -0.29711481
## Induction of Cell-Cell Fusion -0.35150239
## t_median
## Metabolism of RNA 0.06137902
## DNA Repair 0.06415508
## SLC-mediated transmembrane transport -0.12101471
## Neuronal System -0.10801697
## SARS-CoV-1 targets host intracellular signalling and regulatory pathways -0.32442357
## Induction of Cell-Cell Fusion -0.35282317
## p.value
## Metabolism of RNA 1.827231e-05
## DNA Repair 6.657578e-04
## SLC-mediated transmembrane transport 7.674414e-04
## Neuronal System 8.222341e-04
## SARS-CoV-1 targets host intracellular signalling and regulatory pathways 8.659384e-04
## Induction of Cell-Cell Fusion 1.055422e-03
## logp
## Metabolism of RNA 4.738207
## DNA Repair 3.176684
## SLC-mediated transmembrane transport 3.114955
## Neuronal System 3.085005
## SARS-CoV-1 targets host intracellular signalling and regulatory pathways 3.062513
## Induction of Cell-Cell Fusion 2.976574
## fdr
## Metabolism of RNA 0.02965596
## DNA Repair 0.25524290
## SLC-mediated transmembrane transport 0.25524290
## Neuronal System 0.25524290
## SARS-CoV-1 targets host intracellular signalling and regulatory pathways 0.25524290
## Induction of Cell-Cell Fusion 0.25524290Session 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] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] pkgload_1.3.0
## [2] GGally_2.1.2
## [3] ggplot2_3.3.6
## [4] reshape2_1.4.4
## [5] beeswarm_0.4.0
## [6] gplots_3.1.3
## [7] gtools_3.9.3
## [8] tibble_3.1.8
## [9] echarts4r_0.4.4
## [10] tictoc_1.1
## [11] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0
## [12] minfi_1.42.0
## [13] bumphunter_1.38.0
## [14] locfit_1.5-9.6
## [15] iterators_1.0.14
## [16] foreach_1.5.2
## [17] Biostrings_2.64.1
## [18] XVector_0.36.0
## [19] SummarizedExperiment_1.26.1
## [20] Biobase_2.56.0
## [21] MatrixGenerics_1.8.1
## [22] matrixStats_0.62.0
## [23] GenomicRanges_1.48.0
## [24] GenomeInfoDb_1.32.4
## [25] IRanges_2.30.1
## [26] S4Vectors_0.34.0
## [27] BiocGenerics_0.42.0
## [28] mitch_1.8.0
## [29] eulerr_6.1.1
## [30] kableExtra_1.3.4
## [31] dplyr_1.0.10
##
## loaded via a namespace (and not attached):
## [1] BiocFileCache_2.4.0 systemfonts_1.0.4
## [3] plyr_1.8.7 splines_4.2.1
## [5] BiocParallel_1.30.3 digest_0.6.29
## [7] htmltools_0.5.3 fansi_1.0.3
## [9] magrittr_2.0.3 memoise_2.0.1
## [11] tzdb_0.3.0 limma_3.52.4
## [13] readr_2.1.3 annotate_1.74.0
## [15] svglite_2.1.0 askpass_1.1
## [17] siggenes_1.70.0 prettyunits_1.1.1
## [19] colorspace_2.0-3 blob_1.2.3
## [21] rvest_1.0.3 rappdirs_0.3.3
## [23] xfun_0.33 crayon_1.5.2
## [25] RCurl_1.98-1.9 jsonlite_1.8.2
## [27] genefilter_1.78.0 GEOquery_2.64.2
## [29] survival_3.4-0 glue_1.6.2
## [31] gtable_0.3.1 zlibbioc_1.42.0
## [33] webshot_0.5.4 DelayedArray_0.22.0
## [35] Rhdf5lib_1.18.2 HDF5Array_1.24.2
## [37] scales_1.2.1 DBI_1.1.3
## [39] rngtools_1.5.2 Rcpp_1.0.9
## [41] viridisLite_0.4.1 xtable_1.8-4
## [43] progress_1.2.2 bit_4.0.4
## [45] mclust_5.4.10 preprocessCore_1.58.0
## [47] htmlwidgets_1.5.4 httr_1.4.4
## [49] RColorBrewer_1.1-3 ellipsis_0.3.2
## [51] pkgconfig_2.0.3 reshape_0.8.9
## [53] XML_3.99-0.11 sass_0.4.2
## [55] dbplyr_2.2.1 utf8_1.2.2
## [57] tidyselect_1.2.0 rlang_1.0.6
## [59] later_1.3.0 AnnotationDbi_1.58.0
## [61] munsell_0.5.0 tools_4.2.1
## [63] cachem_1.0.6 cli_3.4.1
## [65] generics_0.1.3 RSQLite_2.2.18
## [67] evaluate_0.17 stringr_1.4.1
## [69] fastmap_1.1.0 yaml_2.3.5
## [71] knitr_1.40 bit64_4.0.5
## [73] beanplot_1.3.1 scrime_1.3.5
## [75] caTools_1.18.2 purrr_0.3.5
## [77] KEGGREST_1.36.3 nlme_3.1-160
## [79] doRNG_1.8.2 sparseMatrixStats_1.8.0
## [81] mime_0.12 nor1mix_1.3-0
## [83] xml2_1.3.3 biomaRt_2.52.0
## [85] compiler_4.2.1 rstudioapi_0.14
## [87] filelock_1.0.2 curl_4.3.3
## [89] png_0.1-7 bslib_0.4.0
## [91] stringi_1.7.8 highr_0.9
## [93] GenomicFeatures_1.48.4 lattice_0.20-45
## [95] Matrix_1.5-1 multtest_2.52.0
## [97] vctrs_0.4.2 pillar_1.8.1
## [99] lifecycle_1.0.3 rhdf5filters_1.8.0
## [101] jquerylib_0.1.4 data.table_1.14.2
## [103] bitops_1.0-7 httpuv_1.6.6
## [105] rtracklayer_1.56.1 R6_2.5.1
## [107] BiocIO_1.6.0 promises_1.2.0.1
## [109] KernSmooth_2.23-20 gridExtra_2.3
## [111] codetools_0.2-18 MASS_7.3-58.1
## [113] assertthat_0.2.1 rhdf5_2.40.0
## [115] openssl_2.0.3 rjson_0.2.21
## [117] withr_2.5.0 GenomicAlignments_1.32.1
## [119] Rsamtools_2.12.0 GenomeInfoDbData_1.2.8
## [121] hms_1.1.2 quadprog_1.5-8
## [123] grid_4.2.1 tidyr_1.2.1
## [125] base64_2.0.1 rmarkdown_2.17
## [127] DelayedMatrixStats_1.18.1 illuminaio_0.38.0
## [129] shiny_1.7.2 restfulr_0.0.15