Pathway analysis powered with R/bioconductor
Mark Ziemann & Anusuiya Bora
2026-05-04
Source: https://github.com/markziemann/pathway-workshop
Introduction
#knitr::opts_chunk$set(dev = 'svg') # set output device to svg
suppressPackageStartupMessages({
library("getDEE2")
library("DESeq2")
library("kableExtra")
library("eulerr")
library("gplots")
library("mitch")
library("fgsea")
library("RhpcBLASctl")
library("parallel")
library("beeswarm")
library("vioplot")
})
RhpcBLASctl::blas_set_num_threads(1)Here we demonstrate how to do over-representation analysis (ORA) and functional class scoring (FCS).
Load DE data
For information on the dataset we are using, see the script called “dataprep.Rmd”.
de <- readRDS("data/de.Rds")
head(de) %>%
kbl(caption="To DEGs for Aza treatment in AML3 cells") %>%
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | SRR1171523 | SRR1171524 | SRR1171525 | SRR1171526 | SRR1171527 | SRR1171528 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000165949 IFI27 | 1960.1970 | -3.384492 | 0.0938869 | -36.04861 | 0 | 0 | 4689 | 3583 | 2758 | 309 | 384 | 334 |
| ENSG00000090382 LYZ | 7596.0299 | -1.650342 | 0.0561143 | -29.41036 | 0 | 0 | 14212 | 10237 | 10789 | 3476 | 3764 | 3738 |
| ENSG00000115461 IGFBP5 | 531.2217 | -5.071157 | 0.1795239 | -28.24781 | 0 | 0 | 1320 | 823 | 1025 | 23 | 26 | 43 |
| ENSG00000157601 MX1 | 827.1511 | -2.877795 | 0.1047823 | -27.46450 | 0 | 0 | 1732 | 1501 | 1206 | 184 | 223 | 186 |
| ENSG00000111331 OAS3 | 2127.2010 | -2.661214 | 0.0972124 | -27.37525 | 0 | 0 | 4204 | 3977 | 2972 | 562 | 614 | 560 |
| ENSG00000070915 SLC12A3 | 424.5509 | -3.374852 | 0.1298671 | -25.98697 | 0 | 0 | 1012 | 721 | 653 | 63 | 85 | 76 |
Gene sets were obtained from Reactome.
gs <- gmtPathways("ref/ReactomePathways_2025-09-02.gmt")
head(gs,3)## $`2-LTR circle formation`
## [1] "BANF1" "HMGA1" "LIG4" "PSIP1" "XRCC4" "XRCC5" "XRCC6"
## [8] "gag" "gag-pol" "rev" "vif" "vpr" "vpu"
##
## $`3-Methylcrotonyl-CoA carboxylase deficiency`
## [1] "MCCC1" "MCCC2"
##
## $`3-hydroxyisobutyryl-CoA hydrolase deficiency`
## [1] "HIBCH"summary(unlist(lapply(gs,length)))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 5.00 15.00 48.61 43.00 2617.00ORA: Filter genes
Here we select up and down-regulated genes as well as the background. The background was previously selected based on the presence of an average of 10 reads per sample.
up <- rownames(subset(de, padj < 0.01 & log2FoldChange > 0))
up <- unique(sapply(strsplit(up," "),"[[",2))
length(up)## [1] 1066dn <- rownames(subset(de, padj < 0.01 & log2FoldChange < 0))
dn <- unique(sapply(strsplit(dn," "),"[[",2))
length(dn)## [1] 1432# background selection was done earlier
#mxf <- mx[which(rowMeans(mx)>=10),]
dim(de)## [1] 13168 12bg <- rownames(de)
bg <- unique(sapply(strsplit(bg," "),"[[",2))
length(bg)## [1] 13164Now we write the gene lists to text files that we can use for ORA with web tools like iDEP.
writeLines(up,"data/up.txt")
writeLines(dn,"data/down.txt")
writeLines(bg,"data/background.txt")ORA: run the test with reactome gene sets
The goal is to understand whether certain gene pathways are enriched at a higher level in the DE list as compared to the background list of detected genes.
We are running separate tests for up and down DEGs.
First, up.
ora_up <- fora(pathways=gs, genes=up, universe=bg, minSize = 5, maxSize = Inf)
nrow(ora_up)## [1] 1838message("number of up-regulated gene sets with FDR<0.01")## number of up-regulated gene sets with FDR<0.01nrow(subset(ora_up,padj<0.01))## [1] 145head(ora_up,10) %>%
kbl(caption="Top upregulated pathways with ORA (significance)") %>%
kable_paper("hover", full_width = F)| pathway | pval | padj | foldEnrichment | overlap | size | overlapGenes |
|---|---|---|---|---|---|---|
| Cell Cycle | 0 | 0 | 2.964605 | 139 | 579 | AHCTF1, …. |
| Cell Cycle, Mitotic | 0 | 0 | 3.113629 | 118 | 468 | AHCTF1, …. |
| Metabolism of RNA | 0 | 0 | 2.207637 | 128 | 716 | ADAT1, A…. |
| M Phase | 0 | 0 | 2.916781 | 77 | 326 | AHCTF1, …. |
| Cell Cycle Checkpoints | 0 | 0 | 3.188463 | 63 | 244 | AHCTF1, …. |
| Mitotic Metaphase and Anaphase | 0 | 0 | 3.205982 | 54 | 208 | AHCTF1, …. |
| Mitotic Anaphase | 0 | 0 | 3.161813 | 53 | 207 | AHCTF1, …. |
| Mitotic Prometaphase | 0 | 0 | 3.218614 | 49 | 188 | AHCTF1, …. |
| DNA Repair | 0 | 0 | 2.744215 | 60 | 270 | ACTL6A, …. |
| Mitotic G1 phase and G1/S transition | 0 | 0 | 3.450447 | 38 | 136 | AKT3, CC…. |
ora_up <- ora_up[order(-ora_up$foldEnrichment),]
head(ora_up,10) %>%
kbl(caption="Top upregulated pathways with ORA (fold enrichment)") %>%
kable_paper("hover", full_width = F)| pathway | pval | padj | foldEnrichment | overlap | size | overlapGenes |
|---|---|---|---|---|---|---|
| Condensation of Prometaphase Chromosomes | 0.0000002 | 0.0000125 | 8.981068 | 8 | 11 | CCNB1, C…. |
| Phosphorylation of Emi1 | 0.0005616 | 0.0075672 | 8.232645 | 4 | 6 | CCNB1, C…. |
| Activation of NOXA and translocation to mitochondria | 0.0046755 | 0.0434519 | 7.409381 | 3 | 5 | PMAIP1, …. |
| Aspartate and asparagine metabolism | 0.0046755 | 0.0434519 | 7.409381 | 3 | 5 | GOT2, NA…. |
| TFAP2A acts as a transcriptional repressor during retinoic acid induced cell differentiation | 0.0046755 | 0.0434519 | 7.409381 | 3 | 5 | HSPD1, M…. |
| G2/M DNA replication checkpoint | 0.0012263 | 0.0145416 | 7.056553 | 4 | 7 | CCNA2, C…. |
| G1/S-Specific Transcription | 0.0000000 | 0.0000001 | 6.387397 | 15 | 29 | CCNE1, C…. |
| Condensation of Prophase Chromosomes | 0.0000462 | 0.0008499 | 6.174484 | 7 | 14 | CCNB1, C…. |
| DNA replication initiation | 0.0022956 | 0.0245312 | 6.174484 | 4 | 8 | POLA1, P…. |
| CREB3 factors activate genes | 0.0087893 | 0.0730986 | 6.174484 | 3 | 6 | CREB3L2,…. |
Next, down.
ora_dn <- fora(pathways=gs, genes=dn, universe=bg, minSize = 5, maxSize = Inf)
nrow(ora_dn)## [1] 1838message("number of down-regulated gene sets with FDR<0.01")## number of down-regulated gene sets with FDR<0.01nrow(subset(ora_dn,padj<0.01))## [1] 68head(ora_dn,10) %>%
kbl(caption="Top downregulated pathways with ORA (significance)") %>%
kable_paper("hover", full_width = F)| pathway | pval | padj | foldEnrichment | overlap | size | overlapGenes |
|---|---|---|---|---|---|---|
| Immune System | 0 | 0 | 2.103040 | 326 | 1425 | ACAA1, A…. |
| Cytokine Signaling in Immune system | 0 | 0 | 2.517446 | 155 | 566 | ADAR, AI…. |
| Neutrophil degranulation | 0 | 0 | 2.852918 | 117 | 377 | ACAA1, A…. |
| Innate Immune System | 0 | 0 | 2.242131 | 180 | 738 | ACAA1, A…. |
| Interferon alpha/beta signaling | 0 | 0 | 6.128492 | 40 | 60 | ADAR, BS…. |
| Interferon Signaling | 0 | 0 | 3.206769 | 75 | 215 | ADAR, B2…. |
| Interferon gamma signaling | 0 | 0 | 5.213791 | 38 | 67 | B2M, GBP…. |
| Extracellular matrix organization | 0 | 0 | 3.104041 | 52 | 154 | ACTA2, A…. |
| Signal Transduction | 0 | 0 | 1.485193 | 269 | 1665 | ACTA2, A…. |
| Signaling by Interleukins | 0 | 0 | 2.214104 | 79 | 328 | AIP, ALO…. |
ora_dn <- ora_dn[order(-ora_dn$foldEnrichment),]
head(ora_dn,10) %>%
kbl(caption="Top dnregulated pathways with ORA (fold enrichment)") %>%
kable_paper("hover", full_width = F)| pathway | pval | padj | foldEnrichment | overlap | size | overlapGenes |
|---|---|---|---|---|---|---|
| Interleukin-20 family signaling | 0.0000000 | 0.0000015 | 7.660615 | 10 | 12 | IL10RB, …. |
| Interleukin-9 signaling | 0.0000826 | 0.0033746 | 7.660615 | 5 | 6 | JAK3, ST…. |
| Prostanoid ligand receptors | 0.0006370 | 0.0143829 | 7.354190 | 4 | 5 | PTGER2, …. |
| Endosomal/Vacuolar pathway | 0.0000003 | 0.0000327 | 6.894553 | 9 | 12 | B2M, CTS…. |
| Interleukin-21 signaling | 0.0000378 | 0.0018298 | 6.894553 | 6 | 8 | JAK3, ST…. |
| Interleukin-2 signaling | 0.0002631 | 0.0084842 | 6.566241 | 5 | 7 | JAK3, PT…. |
| Interferon alpha/beta signaling | 0.0000000 | 0.0000000 | 6.128492 | 40 | 60 | ADAR, BS…. |
| STAT5 Activation | 0.0017457 | 0.0320862 | 6.128492 | 4 | 6 | FLT3, GA…. |
| Specification of primordial germ cells | 0.0017457 | 0.0320862 | 6.128492 | 4 | 6 | CXCR4, E…. |
| Signaling by Leptin | 0.0001030 | 0.0040935 | 6.128492 | 6 | 9 | IRS1, JA…. |
Barplot
npw=15
ora_up_vals <- head(ora_up$foldEnrichment,15)
names(ora_up_vals) <- head(ora_up$pathway,15)
ora_dn_vals <- head(ora_dn$foldEnrichment,15)
names(ora_dn_vals) <- head(ora_dn$pathway,15)
myvals <- c(ora_dn_vals,rev(ora_up_vals))
names(myvals) <- ifelse(nchar(names(myvals)) > 53, paste0(strtrim(names(myvals), 50), '...'), names(myvals))
par(mar=c(5.1,25.1,4.1,2.1))
mycols <- c( rep("blue",length(ora_dn_vals)) , rep("red",length(ora_up_vals)) )
barplot(myvals, col=mycols,
horiz=TRUE,las=2,xlab="Fold enrichment",main="Effect of Aza treatment")
mtext("Red=up, Blue=down")
par(mar=c(5.1,4.1,4.1,2.1))Functional class scoring
We will rank the genes based on the “stat” value from DESeq2, which is the Wald test statistic.
First step is to remove redundant gene symbols (there are 4).
de2 <- de
de2$gene <- sapply(strsplit(rownames(de)," "),"[[",2)
which(duplicated(de2$gene))## [1] 6651 9733 9985 10372de2ag <- aggregate(de2, stat ~ gene, FUN=sum)
destat <- de2ag$stat
names(destat) <- de2ag$gene
head(destat)## A1BG-AS1 A2M AAAS AACS AAED1 AAGAB
## -0.7815770 -2.9964032 0.1881224 1.8320025 -0.0862575 -0.2837587fcs_res <- fgsea(pathways = gs,
stats = destat,
minSize = 5,
maxSize = 50000,
nPermSimple = 10000,
nproc = 1)## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%fcs_res <- fcs_res[order(fcs_res$padj),]Now let’s take a look at the top results based on statistical significance.
head( subset(fcs_res,ES>0) , 10) %>%
kbl(caption="Top upregulated pathways with FCS (significance)") %>%
kable_paper("hover", full_width = F)| pathway | pval | padj | log2err | ES | NES | size | leadingEdge |
|---|---|---|---|---|---|---|---|
| Metabolism of RNA | 0 | 0 | 1.684336 | 0.4717290 | 2.491212 | 716 | GEMIN5, …. |
| Cell Cycle | 0 | 0 | 1.640742 | 0.4937691 | 2.594585 | 579 | CCNA2, C…. |
| Cell Cycle, Mitotic | 0 | 0 | 1.488540 | 0.4972344 | 2.554244 | 468 | CCNA2, C…. |
| M Phase | 0 | 0 | 1.287104 | 0.5100265 | 2.534465 | 326 | CENPF, B…. |
| Cell Cycle Checkpoints | 0 | 0 | 1.212545 | 0.5424191 | 2.617614 | 244 | CCNA2, C…. |
| Mitotic Metaphase and Anaphase | 0 | 0 | 1.169070 | 0.5490495 | 2.607559 | 208 | CENPF, B…. |
| Processing of Capped Intron-Containing Pre-mRNA | 0 | 0 | 1.160180 | 0.5106675 | 2.500565 | 272 | DDX46, N…. |
| Mitotic Anaphase | 0 | 0 | 1.133090 | 0.5456177 | 2.588219 | 207 | CENPF, B…. |
| Mitotic Prometaphase | 0 | 0 | 1.114664 | 0.5591510 | 2.628564 | 188 | CENPF, B…. |
| DNA Repair | 0 | 0 | 1.037696 | 0.4731121 | 2.318012 | 270 | CCNA2, P…. |
head( subset(fcs_res,ES<0) , 10) %>%
kbl(caption="Top downregulated pathways with FCS (significance)") %>%
kable_paper("hover", full_width = F)| pathway | pval | padj | log2err | ES | NES | size | leadingEdge |
|---|---|---|---|---|---|---|---|
| Immune System | 0 | 0 | 1.6781779 | -0.4854837 | -2.034965 | 1425 | IFI27, L…. |
| Cytokine Signaling in Immune system | 0 | 0 | 1.3267161 | -0.5420388 | -2.224470 | 566 | IFI27, M…. |
| Innate Immune System | 0 | 0 | 1.2790340 | -0.4971244 | -2.060399 | 738 | LYZ, HLA…. |
| Interferon alpha/beta signaling | 0 | 0 | 1.1690700 | -0.8501596 | -2.790631 | 60 | IFI27, M…. |
| Interferon Signaling | 0 | 0 | 1.1053366 | -0.6221241 | -2.400821 | 215 | IFI27, M…. |
| Neutrophil degranulation | 0 | 0 | 1.0864405 | -0.5429526 | -2.188181 | 377 | LYZ, HLA…. |
| Interferon gamma signaling | 0 | 0 | 1.0175448 | -0.7821332 | -2.611943 | 67 | OAS3, HL…. |
| Extracellular matrix organization | 0 | 0 | 0.9865463 | -0.6342692 | -2.370969 | 154 | COL1A2, …. |
| Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell | 0 | 0 | 0.8513391 | -0.7453242 | -2.446511 | 60 | HLA-B, H…. |
| Modulation of host responses by IFN-stimulated genes | 0 | 0 | 0.8266573 | -0.9166258 | -2.379377 | 18 | IFI27, I…. |
Now let’s take a look at the results based on ES after removal of nonsignificant sets.
sig <- subset(fcs_res,padj<0.01)
fcs_up <- tail(sig[order(sig$ES),],15)
fcs_up %>%
kbl(caption="Top upregulated pathways with FCS (ES prioritisation)") %>%
kable_paper("hover", full_width = F)| pathway | pval | padj | log2err | ES | NES | size | leadingEdge |
|---|---|---|---|---|---|---|---|
| Regulation of Glucokinase by Glucokinase Regulatory Protein | 0.0000000 | 0.0000001 | 0.7749390 | 0.7718391 | 2.550402 | 29 | NDC1, NU…. |
| Export of Viral Ribonucleoproteins from Nucleus | 0.0000000 | 0.0000000 | 0.8012156 | 0.7760569 | 2.607806 | 31 | NDC1, NU…. |
| NEP/NS2 Interacts with the Cellular Export Machinery | 0.0000000 | 0.0000000 | 0.8012156 | 0.7760569 | 2.607806 | 31 | NDC1, NU…. |
| Impaired BRCA2 binding to PALB2 | 0.0000000 | 0.0000011 | 0.7195128 | 0.7767850 | 2.444409 | 24 | XRCC2, B…. |
| Transport of Ribonucleoproteins into the Host Nucleus | 0.0000000 | 0.0000000 | 0.8140358 | 0.7776911 | 2.613297 | 31 | KPNB1, N…. |
| mRNA decay by 3’ to 5’ exoribonuclease | 0.0002585 | 0.0023290 | 0.4984931 | 0.7787064 | 2.061422 | 13 | DIS3, HB…. |
| SLBP independent Processing of Histone Pre-mRNAs | 0.0006325 | 0.0050547 | 0.4772708 | 0.7801965 | 1.918543 | 10 | NCBP1, L…. |
| SLBP Dependent Processing of Replication-Dependent Histone Pre-mRNAs | 0.0003973 | 0.0034122 | 0.4984931 | 0.7809511 | 1.982672 | 11 | NCBP1, L…. |
| Rev-mediated nuclear export of HIV RNA | 0.0000000 | 0.0000000 | 0.8513391 | 0.7821972 | 2.692677 | 34 | NDC1, NU…. |
| Nuclear import of Rev protein | 0.0000000 | 0.0000000 | 0.8266573 | 0.7824653 | 2.664039 | 33 | KPNB1, N…. |
| Interactions of Rev with host cellular proteins | 0.0000000 | 0.0000000 | 0.8986712 | 0.7869546 | 2.735182 | 36 | KPNB1, N…. |
| Postmitotic nuclear pore complex (NPC) reformation | 0.0000000 | 0.0000001 | 0.7749390 | 0.7968889 | 2.586058 | 27 | KPNB1, N…. |
| Condensation of Prometaphase Chromosomes | 0.0000217 | 0.0002676 | 0.5756103 | 0.8415269 | 2.136462 | 11 | SMC4, NC…. |
| Phosphorylation of Emi1 | 0.0005327 | 0.0043711 | 0.4772708 | 0.8925998 | 1.861919 | 6 | FBXO5, C…. |
| TFAP2A acts as a transcriptional repressor during retinoic acid induced cell differentiation | 0.0000531 | 0.0005913 | 0.5573322 | 0.9442486 | 1.854579 | 5 | HSPD1, M…. |
fcs_dn <- head(sig[order(sig$ES),],15)
fcs_dn %>%
kbl(caption="Top downregulated pathways with FCS (ES prioritisation)") %>%
kable_paper("hover", full_width = F)| pathway | pval | padj | log2err | ES | NES | size | leadingEdge |
|---|---|---|---|---|---|---|---|
| Modulation of host responses by IFN-stimulated genes | 0.0000000 | 0.0000000 | 0.8266573 | -0.9166258 | -2.379377 | 18 | IFI27, I…. |
| Endosomal/Vacuolar pathway | 0.0000032 | 0.0000502 | 0.6272567 | -0.8981315 | -2.111452 | 12 | HLA-B, H…. |
| Constitutive Signaling by NOTCH1 t(7;9)(NOTCH1:M1580_K2555) Translocation Mutant | 0.0008144 | 0.0062894 | 0.4772708 | -0.8967871 | -1.838417 | 7 | NOTCH1, …. |
| Signaling by NOTCH1 t(7;9)(NOTCH1:M1580_K2555) Translocation Mutant | 0.0008144 | 0.0062894 | 0.4772708 | -0.8967871 | -1.838417 | 7 | NOTCH1, …. |
| OAS antiviral response | 0.0010439 | 0.0076746 | 0.4550599 | -0.8621346 | -1.834117 | 8 | OAS3, OA…. |
| Interferon alpha/beta signaling | 0.0000000 | 0.0000000 | 1.1690700 | -0.8501596 | -2.790631 | 60 | IFI27, M…. |
| Interleukin-20 family signaling | 0.0005883 | 0.0047636 | 0.4772708 | -0.8141392 | -1.913991 | 12 | STAT1, S…. |
| Maturation of nucleoprotein_9683610 | 0.0013013 | 0.0092225 | 0.4550599 | -0.7952852 | -1.824819 | 11 | PARP9, P…. |
| MyD88 deficiency (TLR2/4) | 0.0008983 | 0.0068511 | 0.4772708 | -0.7914249 | -1.899945 | 13 | MYD88, S…. |
| Interferon gamma signaling | 0.0000000 | 0.0000000 | 1.0175448 | -0.7821332 | -2.611943 | 67 | OAS3, HL…. |
| Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell | 0.0000000 | 0.0000000 | 0.8513391 | -0.7453242 | -2.446511 | 60 | HLA-B, H…. |
| Collagen degradation | 0.0000392 | 0.0004530 | 0.5573322 | -0.7328445 | -2.071703 | 27 | COL1A2, …. |
| Amyloid fiber formation | 0.0000125 | 0.0001663 | 0.5933255 | -0.7306513 | -2.122351 | 31 | LYZ, TGF…. |
| Antigen Presentation: Folding, assembly and peptide loading of class I MHC | 0.0000120 | 0.0001652 | 0.5933255 | -0.7306172 | -2.113575 | 30 | HLA-B, H…. |
| Smooth Muscle Contraction | 0.0001122 | 0.0011145 | 0.5384341 | -0.7191355 | -2.018117 | 26 | TPM2, TP…. |
Now make a barplot.
fcs_vals <- c(fcs_dn$ES, fcs_up$ES)
names(fcs_vals) <- c(fcs_dn$pathway, fcs_up$pathway)
names(fcs_vals) <- ifelse(nchar(names(fcs_vals)) > 53, paste0(strtrim(names(fcs_vals), 50), '...'), names(fcs_vals))
par(mar=c(5.1,24.1,4.1,2.1))
barplot(fcs_vals, horiz=TRUE,las=2,xlab="ES",main="Effect of Aza treatment")
par(mar=c(5.1,4.1,4.1,2.1))Compare ORA and FCS approaches
ora_sig_up <- subset(ora_up,padj<0.01)$pathway
ora_sig_dn <- subset(ora_dn,padj<0.01)$pathway
fcs_sig_up <- subset(fcs_res,ES>0 & padj<0.01)$pathway
fcs_sig_dn <- subset(fcs_res,ES<0 & padj<0.01)$pathway
v1 <- list("ORA up"=ora_sig_up, "ORA down"=ora_sig_dn,
"FCS up"=fcs_sig_up, "FCS down"=fcs_sig_dn)
plot(euler(v1),quantities = list(cex = 2), labels = list(cex = 2))
Set analysis.
message("Common upregulated")## Common upregulatedintersect(ora_sig_up, fcs_sig_up)## [1] "Condensation of Prometaphase Chromosomes"
## [2] "Phosphorylation of Emi1"
## [3] "G1/S-Specific Transcription"
## [4] "Condensation of Prophase Chromosomes"
## [5] "Interactions of Rev with host cellular proteins"
## [6] "Nuclear import of Rev protein"
## [7] "G0 and Early G1"
## [8] "Postmitotic nuclear pore complex (NPC) reformation"
## [9] "Rev-mediated nuclear export of HIV RNA"
## [10] "Polo-like kinase mediated events"
## [11] "Nuclear Pore Complex (NPC) Disassembly"
## [12] "Transcription of E2F targets under negative control by DREAM complex"
## [13] "Export of Viral Ribonucleoproteins from Nucleus"
## [14] "NEP/NS2 Interacts with the Cellular Export Machinery"
## [15] "Transport of Ribonucleoproteins into the Host Nucleus"
## [16] "Impaired BRCA2 binding to PALB2"
## [17] "Defective TPR may confer susceptibility towards thyroid papillary carcinoma (TPC)"
## [18] "Regulation of Glucokinase by Glucokinase Regulatory Protein"
## [19] "NS1 Mediated Effects on Host Pathways"
## [20] "Defective HDR through Homologous Recombination Repair (HRR) due to PALB2 loss of BRCA1 binding function"
## [21] "Defective HDR through Homologous Recombination Repair (HRR) due to PALB2 loss of BRCA2/RAD51/RAD51C binding function"
## [22] "Defective homologous recombination repair (HRR) due to BRCA1 loss of function"
## [23] "Defective homologous recombination repair (HRR) due to PALB2 loss of function"
## [24] "Vpr-mediated nuclear import of PICs"
## [25] "Transport of the SLBP independent Mature mRNA"
## [26] "SUMOylation of DNA replication proteins"
## [27] "Transport of the SLBP Dependant Mature mRNA"
## [28] "Resolution of D-loop Structures through Synthesis-Dependent Strand Annealing (SDSA)"
## [29] "Amino acid transport across the plasma membrane"
## [30] "Interactions of Vpr with host cellular proteins"
## [31] "SUMOylation of SUMOylation proteins"
## [32] "Defective homologous recombination repair (HRR) due to BRCA2 loss of function"
## [33] "Diseases of DNA Double-Strand Break Repair"
## [34] "Activation of the pre-replicative complex"
## [35] "Transport of Mature mRNA Derived from an Intronless Transcript"
## [36] "Diseases of DNA repair"
## [37] "Recognition of DNA damage by PCNA-containing replication complex"
## [38] "Homologous DNA Pairing and Strand Exchange"
## [39] "Transport of Mature mRNAs Derived from Intronless Transcripts"
## [40] "Nuclear Envelope Breakdown"
## [41] "Cyclin A/B1/B2 associated events during G2/M transition"
## [42] "SUMOylation of ubiquitinylation proteins"
## [43] "Viral Messenger RNA Synthesis"
## [44] "tRNA processing in the nucleus"
## [45] "Metabolism of non-coding RNA"
## [46] "snRNP Assembly"
## [47] "Presynaptic phase of homologous DNA pairing and strand exchange"
## [48] "Activation of ATR in response to replication stress"
## [49] "Amplification of signal from unattached kinetochores via a MAD2 inhibitory signal"
## [50] "Amplification of signal from the kinetochores"
## [51] "SUMOylation of RNA binding proteins"
## [52] "Transcriptional regulation by small RNAs"
## [53] "Resolution of D-loop Structures through Holliday Junction Intermediates"
## [54] "Resolution of Sister Chromatid Cohesion"
## [55] "Impaired BRCA2 binding to RAD51"
## [56] "Resolution of D-Loop Structures"
## [57] "Transport of Mature mRNA derived from an Intron-Containing Transcript"
## [58] "G2/M DNA damage checkpoint"
## [59] "Gene Silencing by RNA"
## [60] "HDR through Homologous Recombination (HRR)"
## [61] "Nuclear Envelope (NE) Reassembly"
## [62] "SUMOylation of chromatin organization proteins"
## [63] "G1/S Transition"
## [64] "DNA strand elongation"
## [65] "Mitotic G1 phase and G1/S transition"
## [66] "EML4 and NUDC in mitotic spindle formation"
## [67] "Mitotic Spindle Checkpoint"
## [68] "Dual Incision in GG-NER"
## [69] "TP53 Regulates Transcription of DNA Repair Genes"
## [70] "Transport of Mature Transcript to Cytoplasm"
## [71] "G2/M Checkpoints"
## [72] "Mitotic Prophase"
## [73] "Gap-filling DNA repair synthesis and ligation in TC-NER"
## [74] "rRNA modification in the nucleus and cytosol"
## [75] "Mitotic Prometaphase"
## [76] "DNA Double-Strand Break Repair"
## [77] "Mitotic Metaphase and Anaphase"
## [78] "Cell Cycle Checkpoints"
## [79] "Mitotic Anaphase"
## [80] "Cell Cycle, Mitotic"
## [81] "RHO GTPases Activate Formins"
## [82] "Homology Directed Repair"
## [83] "HCMV Early Events"
## [84] "Dual incision in TC-NER"
## [85] "HDR through Homologous Recombination (HRR) or Single Strand Annealing (SSA)"
## [86] "Cell Cycle"
## [87] "SUMOylation of DNA damage response and repair proteins"
## [88] "M Phase"
## [89] "DNA Replication"
## [90] "Synthesis of DNA"
## [91] "HCMV Late Events"
## [92] "Separation of Sister Chromatids"
## [93] "Transcription-Coupled Nucleotide Excision Repair (TC-NER)"
## [94] "SUMOylation"
## [95] "S Phase"
## [96] "APC/C:Cdc20 mediated degradation of mitotic proteins"
## [97] "DNA Repair"
## [98] "tRNA processing"
## [99] "Cellular response to heat stress"
## [100] "Regulation of APC/C activators between G1/S and early anaphase"
## [101] "Activation of APC/C and APC/C:Cdc20 mediated degradation of mitotic proteins"
## [102] "SUMO E3 ligases SUMOylate target proteins"
## [103] "Host Interactions of HIV factors"
## [104] "APC/C-mediated degradation of cell cycle proteins"
## [105] "Regulation of mitotic cell cycle"
## [106] "DNA Replication Pre-Initiation"
## [107] "Global Genome Nucleotide Excision Repair (GG-NER)"
## [108] "Chromosome Maintenance"
## [109] "Late Phase of HIV Life Cycle"
## [110] "G2/M Transition"
## [111] "HIV Life Cycle"
## [112] "Mitotic G2-G2/M phases"
## [113] "Nucleotide Excision Repair"
## [114] "Processing of Capped Intron-Containing Pre-mRNA"
## [115] "Metabolism of RNA"
## [116] "mRNA Splicing - Major Pathway"
## [117] "RHO GTPase Effectors"
## [118] "mRNA Splicing"
## [119] "HIV Infection"
## [120] "Regulation of TP53 Activity"
## [121] "Transcriptional Regulation by TP53"
## [122] "Cilium Assembly"
## [123] "rRNA processing in the nucleus and cytosol"
## [124] "rRNA processing"
## [125] "Gene expression (Transcription)"
## [126] "RNA Polymerase II Transcription"message("Common downregulated")## Common downregulatedintersect(ora_sig_dn, fcs_sig_dn)## [1] "Interleukin-20 family signaling"
## [2] "Endosomal/Vacuolar pathway"
## [3] "Interferon alpha/beta signaling"
## [4] "Interferon gamma signaling"
## [5] "Modulation of host responses by IFN-stimulated genes"
## [6] "MyD88 deficiency (TLR2/4)"
## [7] "Interleukin-10 signaling"
## [8] "Interleukin-2 family signaling"
## [9] "Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell"
## [10] "Collagen degradation"
## [11] "Interleukin-4 and Interleukin-13 signaling"
## [12] "Antigen Presentation: Folding, assembly and peptide loading of class I MHC"
## [13] "Integrin cell surface interactions"
## [14] "Regulation of Insulin-like Growth Factor (IGF) transport and uptake by Insulin-like Growth Factor Binding Proteins (IGFBPs)"
## [15] "Post-translational protein phosphorylation"
## [16] "Amyloid fiber formation"
## [17] "Interferon Signaling"
## [18] "Extracellular matrix organization"
## [19] "Class A/1 (Rhodopsin-like receptors)"
## [20] "Neutrophil degranulation"
## [21] "Degradation of the extracellular matrix"
## [22] "Cell surface interactions at the vascular wall"
## [23] "Antigen processing-Cross presentation"
## [24] "Platelet degranulation"
## [25] "ER-Phagosome pathway"
## [26] "Cytokine Signaling in Immune system"
## [27] "GPCR ligand binding"
## [28] "Response to elevated platelet cytosolic Ca2+"
## [29] "RSV-host interactions"
## [30] "Innate Immune System"
## [31] "Sensory Perception"
## [32] "Signaling by Interleukins"
## [33] "GPCR downstream signalling"
## [34] "Immune System"
## [35] "Signaling by GPCR"
## [36] "Platelet activation, signaling and aggregation"
## [37] "Hemostasis"
## [38] "Signaling by Receptor Tyrosine Kinases"
## [39] "Signal Transduction"
## [40] "Developmental Biology"
## [41] "Adaptive Immune System"
## [42] "Disease"Enrichment plot.
mypw <- gs[which(names(gs) %in% fcs_dn$pathway)]
plotGseaTable(mypw, destat, fcs_res, gseaParam = 0.5)
Session information
sessionInfo()## R version 4.5.3 (2026-03-11)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 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.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] vioplot_0.5.1 zoo_1.8-15
## [3] sm_2.2-6.0 beeswarm_0.4.0
## [5] RhpcBLASctl_0.23-42 fgsea_1.34.2
## [7] mitch_1.20.0 gplots_3.3.0
## [9] eulerr_7.0.4 kableExtra_1.4.0
## [11] DESeq2_1.48.2 SummarizedExperiment_1.38.1
## [13] Biobase_2.68.0 MatrixGenerics_1.20.0
## [15] matrixStats_1.5.0 GenomicRanges_1.60.0
## [17] GenomeInfoDb_1.44.3 IRanges_2.42.0
## [19] S4Vectors_0.46.0 BiocGenerics_0.54.1
## [21] generics_0.1.4 getDEE2_1.18.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-9 gridExtra_2.3 echarts4r_0.5.0
## [4] rlang_1.1.7 magrittr_2.0.4 otel_0.2.0
## [7] compiler_4.5.3 polylabelr_1.0.0 systemfonts_1.3.1
## [10] vctrs_0.7.1 reshape2_1.4.5 htm2txt_2.2.2
## [13] stringr_1.6.0 pkgconfig_2.0.3 crayon_1.5.3
## [16] fastmap_1.2.0 XVector_0.48.0 labeling_0.4.3
## [19] caTools_1.18.3 promises_1.5.0 rmarkdown_2.30
## [22] UCSC.utils_1.4.0 network_1.20.0 purrr_1.2.1
## [25] xfun_0.56 cachem_1.1.0 jsonlite_2.0.0
## [28] later_1.4.6 DelayedArray_0.34.1 BiocParallel_1.42.2
## [31] R6_2.6.1 bslib_0.10.0 stringi_1.8.7
## [34] RColorBrewer_1.1-3 GGally_2.4.0 jquerylib_0.1.4
## [37] Rcpp_1.1.1 knitr_1.51 httpuv_1.6.16
## [40] Matrix_1.7-4 tidyselect_1.2.1 yaml_2.3.12
## [43] rstudioapi_0.18.0 dichromat_2.0-0.1 abind_1.4-8
## [46] codetools_0.2-20 lattice_0.22-9 tibble_3.3.1
## [49] plyr_1.8.9 withr_3.0.2 shiny_1.12.1
## [52] S7_0.2.1 coda_0.19-4.1 evaluate_1.0.5
## [55] polyclip_1.10-7 ggstats_0.12.0 xml2_1.5.2
## [58] pillar_1.11.1 KernSmooth_2.23-26 ggplot2_4.0.2
## [61] scales_1.4.0 gtools_3.9.5 xtable_1.8-4
## [64] glue_1.8.0 tools_4.5.3 data.table_1.18.2.1
## [67] locfit_1.5-9.12 fastmatch_1.1-8 cowplot_1.2.0
## [70] grid_4.5.3 tidyr_1.3.2 GenomeInfoDbData_1.2.14
## [73] cli_3.6.5 textshaping_1.0.4 S4Arrays_1.8.1
## [76] viridisLite_0.4.3 svglite_2.2.2 dplyr_1.2.0
## [79] gtable_0.3.6 sass_0.4.10 digest_0.6.39
## [82] SparseArray_1.8.1 htmlwidgets_1.6.4 farver_2.1.2
## [85] htmltools_0.5.9 lifecycle_1.0.5 httr_1.4.8
## [88] statnet.common_4.13.0 mime_0.13 MASS_7.3-65save.image("data/session2.Rdata")