Introduction

Testing whether ORA biases towards gene sets based on their basemean expression levels.

library("fgsea")
library("kableExtra")
library("beeswarm")
library("eulerr")
library("tictoc")
library("parallel")

RhpcBLASctl::blas_set_num_threads(1)

Read data and make a smear plot.

x <- readRDS("dataprep/bulkrna3.Rds")

x$rank <- rank(x$stat)

head(x)
##                         baseMean log2FoldChange     lfcSE      stat
## ENSG00000000003 TSPAN6  933.9853     -0.5835878 0.2016855 -2.893553
## ENSG00000000419 DPM1    878.3732      1.3127388 0.1567541  8.374509
## ENSG00000000457 SCYL3   160.3316      0.5569394 0.2128733  2.616296
## ENSG00000000460 FIRRM   342.8626     -0.4963699 0.1950339 -2.545044
## ENSG00000000971 CFH    1226.9676      0.3813575 0.2177729  1.751171
## ENSG00000001036 FUCA2  1587.9857     -0.8921835 0.1501429 -5.942231
##                              pvalue         padj  rank
## ENSG00000000003 TSPAN6 3.809098e-03 7.464152e-03  3797
## ENSG00000000419 DPM1   5.545453e-17 4.929044e-16 14583
## ENSG00000000457 SCYL3  8.888959e-03 1.623368e-02 11024
## ENSG00000000460 FIRRM  1.092641e-02 1.964745e-02  4138
## ENSG00000000971 CFH    7.991647e-02 1.177568e-01  9984
## ENSG00000001036 FUCA2  2.811687e-09 1.271511e-08  1737
dim(x)
## [1] 15380     7
plot(log10(x$baseMean),x$log2FoldChange,pch=19,cex=0.5,
  xlab="log10(baseMean)",ylab="log2(foldchange)",
  ylim=c(-8,8))
sig <- subset(x,padj<0.01)
points(log10(sig$baseMean),sig$log2FoldChange,pch=19,cex=0.5,col="red")

plot(x$log2FoldChange,-log10(x$pvalue),  pch=19,cex=0.5,
  xlab="log2(foldchange)",ylab="-log10(p-value)")
sig <- subset(x,padj<0.01)
points(sig$log2FoldChange,-log10(sig$pvalue),pch=19,cex=0.5,col="red")

Break the data into deciles based on basemean expression.

xs <- x[order(-x$baseMean),]

NCHUNKS = 10

CHUNKSIZE <- round(nrow(xs)/NCHUNKS)

START = 1

xl <- as.list(1:10)

for ( i in 1:NCHUNKS ) {
  START <- ( ( i - 1 )  * CHUNKSIZE ) + 1
  END <- START + CHUNKSIZE - 1
  if ( END > nrow(xs) ) { END <- nrow(xs) }
  xl[[i]] <- rownames(xs[START:END,])
}
names(xl) <- paste("level",1:10)

str(xl)
## List of 10
##  $ level 1 : chr [1:1538] "ENSG00000075624 ACTB" "ENSG00000184009 ACTG1" "ENSG00000156508 EEF1A1" "ENSG00000115414 FN1" ...
##  $ level 2 : chr [1:1538] "ENSG00000198909 MAP3K3" "ENSG00000134590 RTL8C" "ENSG00000184014 DENND5A" "ENSG00000090857 PDPR" ...
##  $ level 3 : chr [1:1538] "ENSG00000115541 HSPE1" "ENSG00000141084 RANBP10" "ENSG00000120053 GOT1" "ENSG00000071994 PDCD2" ...
##  $ level 4 : chr [1:1538] "ENSG00000049618 ARID1B" "ENSG00000083520 DIS3" "ENSG00000061794 MRPS35" "ENSG00000134882 UBAC2" ...
##  $ level 5 : chr [1:1538] "ENSG00000086544 ITPKC" "ENSG00000243716 NPIPB5" "ENSG00000167695 TLCD3A" "ENSG00000183048 SLC25A10" ...
##  $ level 6 : chr [1:1538] "ENSG00000206284 WDR46" "ENSG00000107719 PALD1" "ENSG00000115540 MOB4" "ENSG00000136404 TM6SF1" ...
##  $ level 7 : chr [1:1538] "ENSG00000162227 TAF6L" "ENSG00000100288 CHKB" "ENSG00000184226 PCDH9" "ENSG00000164603 BMT2" ...
##  $ level 8 : chr [1:1538] "ENSG00000169193 CCDC126" "ENSG00000175305 CCNE2" "ENSG00000142494 SLC47A1" "ENSG00000136492 BRIP1" ...
##  $ level 9 : chr [1:1538] "ENSG00000185220 PGBD2" "ENSG00000266728 " "ENSG00000224320 HLA-A" "ENSG00000229215 HLA-A" ...
##  $ level 10: chr [1:1538] "ENSG00000214309 MBLAC1" "ENSG00000196793 ZNF239" "ENSG00000274570 SPDYE10" "ENSG00000139173 TMEM117" ...

ORA p-value only

xp <- x[order(x$pvalue),]

up <- rownames(head(subset(xp,log2FoldChange>0),800))
dn <- rownames(head(subset(xp,log2FoldChange<0),800))
bg <- rownames(xp)
ora_up <- fora(pathways=xl, genes=up, universe=bg, minSize = 5, maxSize = Inf)
ora_dn <- fora(pathways=xl, genes=dn, universe=bg, minSize = 5, maxSize = Inf)

ora_up[,1:6] %>%
  kbl(caption="ORA of upregulated DEGs based on p-value against basemean gene categories",row.names=FALSE) %>%
  kable_paper("hover", full_width = F)
ORA of upregulated DEGs based on p-value against basemean gene categories
pathway pval padj foldEnrichment overlap size
level 7 0.0000000 0.0000000 1.7625 141 1538
level 8 0.0000000 0.0000000 1.7125 137 1538
level 6 0.0003078 0.0010261 1.3750 110 1538
level 5 0.1259363 0.3148408 1.1250 90 1538
level 2 0.3761235 0.6268725 1.0375 83 1538
level 4 0.3761235 0.6268725 1.0375 83 1538
level 1 0.6597047 0.8246309 0.9625 77 1538
level 3 0.6597047 0.8246309 0.9625 77 1538
level 9 1.0000000 1.0000000 0.0250 2 1538
level 10 1.0000000 1.0000000 0.0000 0 1538 
ora_dn[,1:6] %>%
  kbl(caption="ORA of downregulated DEGs based on p-value against basemean gene categories",row.names=FALSE) %>%
  kable_paper("hover", full_width = F)
ORA of downregulated DEGs based on p-value against basemean gene categories
pathway pval padj foldEnrichment overlap size
level 1 0.0000000 0.0000000 1.9125 153 1538
level 3 0.0000000 0.0000000 1.8625 149 1538
level 2 0.0000000 0.0000000 1.7375 139 1538
level 5 0.0000553 0.0001383 1.4250 114 1538
level 4 0.0002036 0.0004072 1.3875 111 1538
level 6 0.7826833 1.0000000 0.9250 74 1538
level 7 0.9999999 1.0000000 0.5125 41 1538
level 8 1.0000000 1.0000000 0.2375 19 1538
level 9 1.0000000 1.0000000 0.0000 0 1538
level 10 1.0000000 1.0000000 0.0000 0 1538

ORA p-value and foldchange threshold

up <- rownames(head(subset(xp,log2FoldChange>0.5),800))
dn <- rownames(head(subset(xp,log2FoldChange < -0.5),800))

ora2_up <- fora(pathways=xl, genes=up, universe=bg, minSize = 5, maxSize = Inf)
ora2_dn <- fora(pathways=xl, genes=dn, universe=bg, minSize = 5, maxSize = Inf)

ora2_up[,1:6] %>%
  kbl(caption="ORA of upregulated DEGs based on p-value and logfoldchange threshold against basemean gene categories",row.names=FALSE) %>%
  kable_paper("hover", full_width = F)
ORA of upregulated DEGs based on p-value and logfoldchange threshold against basemean gene categories
pathway pval padj foldEnrichment overlap size
level 7 0.0000000 0.0000000 1.7625 141 1538
level 8 0.0000000 0.0000000 1.7125 137 1538
level 6 0.0003078 0.0010261 1.3750 110 1538
level 5 0.1259363 0.3148408 1.1250 90 1538
level 2 0.3761235 0.6268725 1.0375 83 1538
level 4 0.3761235 0.6268725 1.0375 83 1538
level 1 0.6597047 0.8246309 0.9625 77 1538
level 3 0.6597047 0.8246309 0.9625 77 1538
level 9 1.0000000 1.0000000 0.0250 2 1538
level 10 1.0000000 1.0000000 0.0000 0 1538 
ora2_dn[,1:6] %>%
  kbl(caption="ORA of downregulated DEGs based on p-value and logfoldchange threshold against basemean gene categories",row.names=FALSE) %>%
  kable_paper("hover", full_width = F)
ORA of downregulated DEGs based on p-value and logfoldchange threshold against basemean gene categories
pathway pval padj foldEnrichment overlap size
level 1 0.0000000 0.0000000 1.9125 153 1538
level 3 0.0000000 0.0000000 1.8625 149 1538
level 2 0.0000000 0.0000000 1.7375 139 1538
level 5 0.0000553 0.0001383 1.4250 114 1538
level 4 0.0002036 0.0004072 1.3875 111 1538
level 6 0.7826833 1.0000000 0.9250 74 1538
level 7 0.9999999 1.0000000 0.5125 41 1538
level 8 1.0000000 1.0000000 0.2375 19 1538
level 9 1.0000000 1.0000000 0.0000 0 1538
level 10 1.0000000 1.0000000 0.0000 0 1538

FCS now.

stat <- x$stat

names(stat) <- rownames(x)

fcs1 <- fgsea(pathways=xl, stats=stat,minSize=5,
  maxSize=100000,nPermSimple = 10000)
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (2.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
fcs1 <- fcs1[order(fcs1$pval),]

fcs1[,1:7] %>%
  kbl(caption="FCS of based on Walt test stat against basemean gene categories",row.names=FALSE) %>%
  kable_paper("hover", full_width = F)
FCS of based on Walt test stat against basemean gene categories
pathway pval padj log2err ES NES size
level 10 0.0000000 0.0000000 1.6904724 0.3983051 2.132217 1538
level 9 0.0000000 0.0000000 1.4099514 0.3604833 1.929748 1538
level 1 0.0000000 0.0000000 1.3267161 -0.3663123 -1.886429 1538
level 8 0.0000000 0.0000000 1.2462328 0.3382543 1.810752 1538
level 3 0.0000000 0.0000000 1.0672100 -0.3295711 -1.697220 1538
level 2 0.0000000 0.0000000 0.9436322 -0.3109811 -1.601485 1538
level 5 0.0000001 0.0000001 0.7049757 -0.2761087 -1.421900 1538
level 7 0.0000001 0.0000001 0.7049757 0.2586031 1.384360 1538
level 4 0.0001412 0.0001569 0.5188481 -0.2492018 -1.283335 1538
level 6 0.4204583 0.4204583 0.0379022 0.1877465 1.005049 1538

Idea for calculating amount of significant terms:

ora1bias <- sum(-log10(ora_up$pval)) + sum(-log10(ora_dn$pval))
ora2bias <- sum(-log10(ora2_up$pval)) + sum(-log10(ora2_dn$pval))
fcs1bias <- sum(-log10(fcs1$pval))

barplot(c("ORA p"=ora1bias,"ORA p FC"=ora2bias,"FCS"=fcs1bias),ylab="Bias measure")

Visualising deciles in a gene rank

START = 1

xr <- as.list(1:10)

for ( i in 1:NCHUNKS ) {
  START <- ( ( i - 1 )  * CHUNKSIZE ) + 1
  END <- START + CHUNKSIZE - 1
  if ( END > nrow(xs) ) { END <- nrow(xs) }
  xr[[i]] <- xs[START:END,"rank"]
}
names(xr) <- paste("level",1:10)

str(xr)
## List of 10
##  $ level 1 : num [1:1538] 13531 14705 2548 2965 206 ...
##  $ level 2 : num [1:1538] 8976 9108 5467 11302 3769 ...
##  $ level 3 : num [1:1538] 3500 14976 13800 1313 11592 ...
##  $ level 4 : num [1:1538] 66 12717 2431 3807 11589 ...
##  $ level 5 : num [1:1538] 1091 10590 236 14029 3231 ...
##  $ level 6 : num [1:1538] 3276 2435 9254 557 6055 ...
##  $ level 7 : num [1:1538] 4453 2173 9458 12930 4082 ...
##  $ level 8 : num [1:1538] 10206 6708 8220 4802 5406 ...
##  $ level 9 : num [1:1538] 2483 3782 7522 7522 9734 ...
##  $ level 10: num [1:1538] 3773 7653 8154 2456 6794 ...
boxplot(xr,cex=0,col="white",ylab="Position in gene rank")
beeswarm(xr,cex=0.25,add=TRUE,pch=19,col="darkgray")

Now quantify the data skew.

Pearson’s Coefficient of Skewness.

lapply(xr,summary)
## $`level 1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1    2211    5090    6239   10046   15380 
## 
## $`level 2`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       2    2582    6402    6919   11126   15379 
## 
## $`level 3`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       3    2682    6419    6906   11142   15372 
## 
## $`level 4`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      13    3084    7182    7333   11525   15367 
## 
## $`level 5`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      34    2866    6624    7018   10753   15360 
## 
## $`level 6`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      50    3376    7368    7504   11479   15336 
## 
## $`level 7`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     113    4280    7852    7958   11722   15273 
## 
## $`level 8`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     309    5258    8786    8720   12487   15136 
## 
## $`level 9`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     827    5965    9212    9074   12541   14609 
## 
## $`level 10`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1949    7088    9440    9234   11754   13720
pcs <- unlist(lapply(xr, function(x) { ( mean(x) - median(x) ) / sd(x) } ) )

pcs
##     level 1     level 2     level 3     level 4     level 5     level 6 
##  0.24877283  0.10997825  0.10406573  0.03252517  0.08545128  0.02985701 
##     level 7     level 8     level 9    level 10 
##  0.02385122 -0.01579670 -0.03758411 -0.07431847
barplot(pcs)

sum(abs(pcs))
## [1] 0.7622008

Consider running Rank ANOVA test for enrichment. Optionally can correct for basemean bias.

df <- x[,"rank",drop=FALSE]
df$baseMean <- rank(x$baseMean)

plot(df)

i=1
gs <- xl[[i]]
df$inset <- as.numeric(rownames(df) %in% gs)
mdl0 <- aov(rank ~ inset , data = df)
p0 <- summary(mdl0)[[1]][1,5]
mdl1 <- aov(rank ~ baseMean + inset, data = df)
p1 <- summary(mdl1)[[1]][2,5]
summary(mdl0)
##                Df    Sum Sq   Mean Sq F value Pr(>F)    
## inset           1 3.599e+09 3.599e+09   184.7 <2e-16 ***
## Residuals   15378 2.996e+11 1.948e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mdl1)
##                Df    Sum Sq   Mean Sq F value Pr(>F)    
## baseMean        1 1.359e+10 1.359e+10 721.717 <2e-16 ***
## inset           1 4.705e+05 4.705e+05   0.025  0.874    
## Residuals   15377 2.896e+11 1.883e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p0
## [1] 7.769292e-42
p1
## [1] 0.8744022
myinset <- subset(df,inset==1)
points(myinset[,1:2],col="red")

df <- x[,"rank",drop=FALSE]
df$baseMean <- scale(log10(x$baseMean))

mclapply(1:length(xl), function(i) {
  gs <- xl[[i]]
  df$inset <- as.numeric(rownames(df) %in% gs)
  mdl0 <- aov(rank ~ inset , data = df)
  p0 <- summary(mdl0)[[1]][1,5]
  mdl1 <- aov(rank ~ baseMean + inset, data = df)
  p1 <- summary(mdl1)[[1]][2,5]
  res <- c("p0"=p0,"p1"=p1)
  return(res)
},mc.cores=10)
## [[1]]
##           p0           p1 
## 7.769292e-42 1.603824e-01 
## 
## [[2]]
##           p0           p1 
## 6.611814e-13 6.207328e-02 
## 
## [[3]]
##           p0           p1 
## 2.754843e-13 2.819922e-01 
## 
## [[4]]
##           p0           p1 
## 0.0008779276 0.4076083238 
## 
## [[5]]
##           p0           p1 
## 3.647856e-10 1.735879e-05 
## 
## [[6]]
##         p0         p1 
## 0.08244166 0.05490743 
## 
## [[7]]
##         p0         p1 
## 0.01282455 0.61413676 
## 
## [[8]]
##           p0           p1 
## 7.934177e-22 9.629952e-04 
## 
## [[9]]
##           p0           p1 
## 3.773612e-38 2.586777e-02 
## 
## [[10]]
##           p0           p1 
## 4.204555e-47 1.274065e-01

Run enrichment test on Reactome gene sets with and without correction.

Relate gene symbols to gene IDs.

x$derank <- rank(x$stat)
x$bmrank <- rank(x$baseMean)
x$symbol <- lapply(strsplit(rownames(x)," "), function(x) { x[length(x) ]} )

Plot DE rank vs basemean rank.

Testing code.

df <- x

plot(df[,c("derank","bmrank")])

reactome <- gmtPathways("c2.cp.reactome.v2023.2.Hs.symbols.gmt")

i=1
gs <- reactome[[i]]
df$inset <- as.numeric(df$symbol %in% gs)
mdl0 <- aov(derank ~ inset , data = df)
p0 <- summary(mdl0)[[1]][1,5]

# basic es calc
mrset <- mean(subset(df,inset==1)$derank)
mrbg <- mean(subset(df,inset!=1)$derank)
tot <- nrow(df)
es0 <- (mrset - mrbg) / (tot/2)

mdl1 <- aov(derank ~ bmrank + inset, data = df)
p1 <- summary(mdl1)[[1]][2,5]
summary(mdl0)
##                Df    Sum Sq  Mean Sq F value Pr(>F)
## inset           1 1.544e+07 15439147   0.783  0.376
## Residuals   15378 3.032e+11 19713593
summary(mdl1)
##                Df    Sum Sq   Mean Sq F value Pr(>F)    
## bmrank          1 1.359e+10 1.359e+10 721.724 <2e-16 ***
## inset           1 3.059e+06 3.059e+06   0.162  0.687    
## Residuals   15377 2.896e+11 1.883e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p0
## [1] 0.3761864
p1
## [1] 0.6869117
myinset <- subset(df,inset==1)
points(myinset[,c("derank","bmrank")],col="red",pch=19)

# corrected es calc
bmranks <- mean(subset(df,inset==1)$bmrank)
df2 <- df[order(df$bmrank),]

idx <- which(df2$inset==1)
block <- round(nrow(df)/100)

mysubsets <- lapply(idx, function(i) {
  df3 <- df2[(i-block):(i+block),]
  df3 <- subset(df3,inset!=1)
  return(df3)
} )

# background based on genes in bmrank vicinity rather than all genes
mrbg1 <- mean(unlist(lapply(idx, function(i) {
  df3 <- df2[(i-block):(i+block),]
  df3 <- subset(df3,inset!=1)
  return(df3$derank)
} )))
tot <- nrow(df)
es1 <- (mrset - mrbg1) / (tot/2)

Now run for all Reactome gene sets.

RhpcBLASctl::blas_set_num_threads(1)
tic()
myres <- mclapply(1:length(reactome), function(i) {
  gs <- reactome[[i]]
  gsname <- names(reactome)[i]
  df$inset <- as.numeric(df$symbol %in% gs)
  NGENES <- length(unique(unlist(subset( df,inset==1 )$symbol)))
  mdl0 <- aov(derank ~ inset , data = df)
  p0 <- summary(mdl0)[[1]][1,5]
  mdl1 <- aov(derank ~ bmrank + inset, data = df)
  p1 <- summary(mdl1)[[1]][2,5]
  mrset <- mean(subset(df,inset==1)$derank)
  mrbg <- mean(subset(df,inset!=1)$derank)
  tot <- nrow(df)
  es0 <- (mrset - mrbg) / (tot/2)
  # corrected es
  df2 <- df[order(df$bmrank),]
  idx <- which(df2$inset==1)
  block <- round(nrow(df)/100)
  mrbg1 <- mean(unlist(lapply(idx, function(j) {
    START=j-block
    if (START < 1) { START=1 }
    END=j+block
    if (END > tot) { END=tot }
    df3 <- df2[START:END,]
    df3 <- subset(df3,inset!=1)
    return(df3$derank)
  } )))
  tot <- nrow(df)
  es1 <- (mrset - mrbg1) / (tot/2)
  res <- c("GeneSet"=gsname,"Ngenes"=NGENES,"ES0"=es0,"ES1"=es1,"p0"=p0,"p1"=p1)
  return(res)
},mc.cores=8)
toc()
## 10.473 sec elapsed
myres <- do.call(rbind,myres)
myrownames <- myres[,1]
myres <- myres[,-1]
myres <- apply(myres,2,as.numeric)
rownames(myres) <- myrownames
myres <- as.data.frame(myres)
myres <- subset(myres,Ngenes>=5)

# FDR correction is vital
myres$fdr0 <- p.adjust(myres$p0,method="fdr")
myres$fdr1 <- p.adjust(myres$p1,method="fdr")

# show top results
myres <- myres[order(myres[,"p0"]),]
head(myres,10) %>%
  kbl(caption="Top results for normal ES calc",row.names=TRUE) %>%
  kable_paper("hover", full_width = F)
Top results for normal ES calc
Ngenes ES0 ES1 p0 p1 fdr0 fdr1
REACTOME_METABOLISM_OF_RNA 706 -0.2480424 -0.1660465 0 0.0e+00 0 0.0000000
REACTOME_TRANSLATION 290 -0.3455556 -0.2322516 0 0.0e+00 0 0.0000000
REACTOME_CELL_CYCLE 625 -0.2240389 -0.1673730 0 0.0e+00 0 0.0000000
REACTOME_CELL_CYCLE_MITOTIC 507 -0.2472756 -0.1853153 0 0.0e+00 0 0.0000000
REACTOME_RRNA_PROCESSING 200 -0.3364327 -0.2168823 0 0.0e+00 0 0.0000028
REACTOME_CELL_CYCLE_CHECKPOINTS 267 -0.2942184 -0.2427673 0 0.0e+00 0 0.0000000
REACTOME_NONSENSE_MEDIATED_DECAY_NMD 113 -0.4325272 -0.2575996 0 2.0e-07 0 0.0000173
REACTOME_SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE 110 -0.4357951 -0.2533666 0 4.0e-07 0 0.0000337
REACTOME_EUKARYOTIC_TRANSLATION_INITIATION 117 -0.4156588 -0.2325140 0 1.8e-06 0 0.0001048
REACTOME_MITOTIC_METAPHASE_AND_ANAPHASE 226 -0.3110235 -0.2344685 0 0.0e+00 0 0.0000002 
myres <- myres[order(myres[,"p1"]),]
head(myres,10) %>%
  kbl(caption="Top results for adjusted ES calc",row.names=TRUE) %>%
  kable_paper("hover", full_width = F)
Top results for adjusted ES calc
Ngenes ES0 ES1 p0 p1 fdr0 fdr1
REACTOME_METABOLISM_OF_RNA 706 -0.2480424 -0.1660465 0.00e+00 0 0.0000000 0.0e+00
REACTOME_CELL_CYCLE 625 -0.2240389 -0.1673730 0.00e+00 0 0.0000000 0.0e+00
REACTOME_CELL_CYCLE_MITOTIC 507 -0.2472756 -0.1853153 0.00e+00 0 0.0000000 0.0e+00
REACTOME_TRANSLATION 290 -0.3455556 -0.2322516 0.00e+00 0 0.0000000 0.0e+00
REACTOME_CELL_CYCLE_CHECKPOINTS 267 -0.2942184 -0.2427673 0.00e+00 0 0.0000000 0.0e+00
REACTOME_MITOTIC_METAPHASE_AND_ANAPHASE 226 -0.3110235 -0.2344685 0.00e+00 0 0.0000000 2.0e-07
REACTOME_M_PHASE 365 -0.2421764 -0.1776868 0.00e+00 0 0.0000000 6.0e-07
REACTOME_RRNA_PROCESSING 200 -0.3364327 -0.2168823 0.00e+00 0 0.0000000 2.8e-06
REACTOME_NEURONAL_SYSTEM 254 0.2038831 0.1919859 0.00e+00 0 0.0000002 3.5e-06
REACTOME_METABOLISM_OF_LIPIDS 597 0.0972610 0.1253196 2.55e-05 0 0.0004205 3.5e-06

Now make an Euler plot.

fdr0set <- rownames(subset(myres,fdr0<0.01))
fdr1set <- rownames(subset(myres,fdr1<0.01))

FDR0LEN <- length(fdr0set)
FDR1LEN <- length(fdr1set)
JACCARD <- signif(length(intersect(fdr0set,fdr1set)) / length(union(fdr0set,fdr1set)),3)
HEADER <- paste("N sig (uncorr)=",FDR0LEN,", N sig (corr)",FDR1LEN,", Jaccard=",JACCARD)
HEADER
## [1] "N sig (uncorr)= 172 , N sig (corr) 87 , Jaccard= 0.392"
v1 <- list("fdr0set"=fdr0set, "fdr1set"=fdr1set)

plot(euler(v1),quantities = list(cex = 2), labels = list(cex = 2))

#mtext(HEADER,cex=0.5)

Now separate up and down.

fdr0up <- rownames(subset(myres,fdr0<0.01&ES0>0))
fdr0dn <- rownames(subset(myres,fdr0<0.01&ES0<0))

fdr1up <- rownames(subset(myres,fdr1<0.01&ES1>0))
fdr1dn <- rownames(subset(myres,fdr1<0.01&ES1<0))

v1 <- list("fdr0up"=fdr0up, "fdr0dn"=fdr0dn, "fdr1up"=fdr1up,"fdr1dn"=fdr1dn)

plot(euler(v1),quantities = list(cex = 2), labels = list(cex = 2))

#mtext(HEADER,cex=0.5)

setdiff(fdr0up,fdr1up)
## [1] "REACTOME_CARDIAC_CONDUCTION"                             
## [2] "REACTOME_GABA_SYNTHESIS_RELEASE_REUPTAKE_AND_DEGRADATION"
## [3] "REACTOME_POTASSIUM_CHANNELS"
setdiff(fdr1up,fdr0up)
##  [1] "REACTOME_SIGNALING_BY_RECEPTOR_TYROSINE_KINASES"                   
##  [2] "REACTOME_VESICLE_MEDIATED_TRANSPORT"                               
##  [3] "REACTOME_SIGNALING_BY_NTRKS"                                       
##  [4] "REACTOME_MEMBRANE_TRAFFICKING"                                     
##  [5] "REACTOME_CLATHRIN_MEDIATED_ENDOCYTOSIS"                            
##  [6] "REACTOME_PHOSPHOLIPID_METABOLISM"                                  
##  [7] "REACTOME_NUCLEAR_EVENTS_KINASE_AND_TRANSCRIPTION_FACTOR_ACTIVATION"
##  [8] "REACTOME_CILIUM_ASSEMBLY"                                          
##  [9] "REACTOME_BRANCHED_CHAIN_AMINO_ACID_CATABOLISM"                     
## [10] "REACTOME_HEPARAN_SULFATE_HEPARIN_HS_GAG_METABOLISM"                
## [11] "REACTOME_INOSITOL_PHOSPHATE_METABOLISM"                            
## [12] "REACTOME_NGF_STIMULATED_TRANSCRIPTION"                             
## [13] "REACTOME_PLASMA_LIPOPROTEIN_ASSEMBLY_REMODELING_AND_CLEARANCE"     
## [14] "REACTOME_HS_GAG_BIOSYNTHESIS"
setdiff(fdr0dn,fdr1dn)
##  [1] "REACTOME_HOMOLOGY_DIRECTED_REPAIR"                                                                                        
##  [2] "REACTOME_VIRAL_INFECTION_PATHWAYS"                                                                                        
##  [3] "REACTOME_ACTIVATION_OF_ATR_IN_RESPONSE_TO_REPLICATION_STRESS"                                                             
##  [4] "REACTOME_MITOCHONDRIAL_PROTEIN_IMPORT"                                                                                    
##  [5] "REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2"                                                                             
##  [6] "REACTOME_STABILIZATION_OF_P53"                                                                                            
##  [7] "REACTOME_MITOTIC_PROPHASE"                                                                                                
##  [8] "REACTOME_SARS_COV_1_MODULATES_HOST_TRANSLATION_MACHINERY"                                                                 
##  [9] "REACTOME_CHK1_CHK2_CDS1_MEDIATED_INACTIVATION_OF_CYCLIN_B_CDK1_COMPLEX"                                                   
## [10] "REACTOME_DNA_DOUBLE_STRAND_BREAK_REPAIR"                                                                                  
## [11] "REACTOME_SARS_COV_2_MODULATES_HOST_TRANSLATION_MACHINERY"                                                                 
## [12] "REACTOME_RHO_GTPASES_ACTIVATE_FORMINS"                                                                                    
## [13] "REACTOME_MITOTIC_G2_G2_M_PHASES"                                                                                          
## [14] "REACTOME_G1_S_DNA_DAMAGE_CHECKPOINTS"                                                                                     
## [15] "REACTOME_TRNA_AMINOACYLATION"                                                                                             
## [16] "REACTOME_MRNA_SPLICING"                                                                                                   
## [17] "REACTOME_ANTIGEN_PROCESSING_UBIQUITINATION_PROTEASOME_DEGRADATION"                                                        
## [18] "REACTOME_DNA_STRAND_ELONGATION"                                                                                           
## [19] "REACTOME_ACTIVATION_OF_THE_MRNA_UPON_BINDING_OF_THE_CAP_BINDING_COMPLEX_AND_EIFS_AND_SUBSEQUENT_BINDING_TO_43S"           
## [20] "REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY"                                                                     
## [21] "REACTOME_CROSS_PRESENTATION_OF_SOLUBLE_EXOGENOUS_ANTIGENS_ENDOSOMES"                                                      
## [22] "REACTOME_HOST_INTERACTIONS_OF_HIV_FACTORS"                                                                                
## [23] "REACTOME_REGULATION_OF_PTEN_STABILITY_AND_ACTIVITY"                                                                       
## [24] "REACTOME_SCF_SKP2_MEDIATED_DEGRADATION_OF_P27_P21"                                                                        
## [25] "REACTOME_CYCLIN_A_CDK2_ASSOCIATED_EVENTS_AT_S_PHASE_ENTRY"                                                                
## [26] "REACTOME_DNA_REPAIR"                                                                                                      
## [27] "REACTOME_TCR_SIGNALING"                                                                                                   
## [28] "REACTOME_NUCLEAR_ENVELOPE_BREAKDOWN"                                                                                      
## [29] "REACTOME_PROGRAMMED_CELL_DEATH"                                                                                           
## [30] "REACTOME_FORMATION_OF_PARAXIAL_MESODERM"                                                                                  
## [31] "REACTOME_THE_ROLE_OF_GTSE1_IN_G2_M_PROGRESSION_AFTER_G2_CHECKPOINT"                                                       
## [32] "REACTOME_TRNA_PROCESSING"                                                                                                 
## [33] "REACTOME_CYTOKINE_SIGNALING_IN_IMMUNE_SYSTEM"                                                                             
## [34] "REACTOME_SCF_BETA_TRCP_MEDIATED_DEGRADATION_OF_EMI1"                                                                      
## [35] "REACTOME_PHOSPHORYLATION_OF_THE_APC_C"                                                                                    
## [36] "REACTOME_INFECTIOUS_DISEASE"                                                                                              
## [37] "REACTOME_REGULATED_NECROSIS"                                                                                              
## [38] "REACTOME_DEGRADATION_OF_BETA_CATENIN_BY_THE_DESTRUCTION_COMPLEX"                                                          
## [39] "REACTOME_NEDDYLATION"                                                                                                     
## [40] "REACTOME_INTERLEUKIN_1_FAMILY_SIGNALING"                                                                                  
## [41] "REACTOME_RUNX1_REGULATES_TRANSCRIPTION_OF_GENES_INVOLVED_IN_DIFFERENTIATION_OF_HSCS"                                      
## [42] "REACTOME_DEGRADATION_OF_DVL"                                                                                              
## [43] "REACTOME_INTERACTIONS_OF_REV_WITH_HOST_CELLULAR_PROTEINS"                                                                 
## [44] "REACTOME_MAPK6_MAPK4_SIGNALING"                                                                                           
## [45] "REACTOME_GOLGI_CISTERNAE_PERICENTRIOLAR_STACK_REORGANIZATION"                                                             
## [46] "REACTOME_ABC_TRANSPORTER_DISORDERS"                                                                                       
## [47] "REACTOME_ADAPTIVE_IMMUNE_SYSTEM"                                                                                          
## [48] "REACTOME_KEAP1_NFE2L2_PATHWAY"                                                                                            
## [49] "REACTOME_DEGRADATION_OF_AXIN"                                                                                             
## [50] "REACTOME_FCERI_MEDIATED_NF_KB_ACTIVATION"                                                                                 
## [51] "REACTOME_SUMOYLATION"                                                                                                     
## [52] "REACTOME_HEDGEHOG_ON_STATE"                                                                                               
## [53] "REACTOME_SARS_COV_INFECTIONS"                                                                                             
## [54] "REACTOME_INITIATION_OF_NUCLEAR_ENVELOPE_NE_REFORMATION"                                                                   
## [55] "REACTOME_INTERLEUKIN_1_SIGNALING"                                                                                         
## [56] "REACTOME_PTEN_REGULATION"                                                                                                 
## [57] "REACTOME_REGULATION_OF_RAS_BY_GAPS"                                                                                       
## [58] "REACTOME_UCH_PROTEINASES"                                                                                                 
## [59] "REACTOME_TRANSCRIPTIONAL_REGULATION_BY_TP53"                                                                              
## [60] "REACTOME_SARS_COV_2_HOST_INTERACTIONS"                                                                                    
## [61] "REACTOME_DOWNSTREAM_SIGNALING_EVENTS_OF_B_CELL_RECEPTOR_BCR"                                                              
## [62] "REACTOME_TNFR2_NON_CANONICAL_NF_KB_PATHWAY"                                                                               
## [63] "REACTOME_DEFECTIVE_CFTR_CAUSES_CYSTIC_FIBROSIS"                                                                           
## [64] "REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1"                                                                             
## [65] "REACTOME_DISORDERS_OF_TRANSMEMBRANE_TRANSPORTERS"                                                                         
## [66] "REACTOME_AUF1_HNRNP_D0_BINDS_AND_DESTABILIZES_MRNA"                                                                       
## [67] "REACTOME_SOMITOGENESIS"                                                                                                   
## [68] "REACTOME_APOPTOSIS"                                                                                                       
## [69] "REACTOME_DEGRADATION_OF_GLI1_BY_THE_PROTEASOME"                                                                           
## [70] "REACTOME_RHOBTB2_GTPASE_CYCLE"                                                                                            
## [71] "REACTOME_SARS_COV_1_TARGETS_HOST_INTRACELLULAR_SIGNALLING_AND_REGULATORY_PATHWAYS"                                        
## [72] "REACTOME_CRISTAE_FORMATION"                                                                                               
## [73] "REACTOME_NUCLEAR_PORE_COMPLEX_NPC_DISASSEMBLY"                                                                            
## [74] "REACTOME_CELLULAR_RESPONSE_TO_HYPOXIA"                                                                                    
## [75] "REACTOME_NUCLEAR_IMPORT_OF_REV_PROTEIN"                                                                                   
## [76] "REACTOME_NUCLEAR_EVENTS_MEDIATED_BY_NFE2L2"                                                                               
## [77] "REACTOME_REGULATION_OF_MRNA_STABILITY_BY_PROTEINS_THAT_BIND_AU_RICH_ELEMENTS"                                             
## [78] "REACTOME_CELLULAR_RESPONSE_TO_STARVATION"                                                                                 
## [79] "REACTOME_RRNA_MODIFICATION_IN_THE_NUCLEUS_AND_CYTOSOL"                                                                    
## [80] "REACTOME_DEUBIQUITINATION"                                                                                                
## [81] "REACTOME_NUCLEAR_ENVELOPE_NE_REASSEMBLY"                                                                                  
## [82] "REACTOME_NS1_MEDIATED_EFFECTS_ON_HOST_PATHWAYS"                                                                           
## [83] "REACTOME_SUMOYLATION_OF_DNA_REPLICATION_PROTEINS"                                                                         
## [84] "REACTOME_DECTIN_1_MEDIATED_NONCANONICAL_NF_KB_SIGNALING"                                                                  
## [85] "REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY"                                                                     
## [86] "REACTOME_RESPIRATORY_ELECTRON_TRANSPORT_ATP_SYNTHESIS_BY_CHEMIOSMOTIC_COUPLING_AND_HEAT_PRODUCTION_BY_UNCOUPLING_PROTEINS"
## [87] "REACTOME_SIGNALING_BY_INTERLEUKINS"                                                                                       
## [88] "REACTOME_UB_SPECIFIC_PROCESSING_PROTEASES"                                                                                
## [89] "REACTOME_SARS_COV_2_INFECTION"                                                                                            
## [90] "REACTOME_HIV_INFECTION"                                                                                                   
## [91] "REACTOME_CELLULAR_RESPONSE_TO_CHEMICAL_STRESS"                                                                            
## [92] "REACTOME_CELLULAR_RESPONSES_TO_STIMULI"                                                                                   
## [93] "REACTOME_SIGNALING_BY_RHO_GTPASES_MIRO_GTPASES_AND_RHOBTB3"                                                               
## [94] "REACTOME_POST_TRANSLATIONAL_PROTEIN_MODIFICATION"                                                                         
## [95] "REACTOME_DEVELOPMENTAL_BIOLOGY"                                                                                           
## [96] "REACTOME_NERVOUS_SYSTEM_DEVELOPMENT"
setdiff(fdr1dn,fdr0dn)
## character(0)

Now take a look at the top results.

myres <- myres[order(myres[,"p0"]),]

top0dn <- rownames(head(subset(myres,ES0<0),10))
top0up <- rownames(head(subset(myres,ES0>0),10))

myres <- myres[order(myres[,"p1"]),]

top1dn <- rownames(head(subset(myres,ES1<0),10))
top1up <- rownames(head(subset(myres,ES1>0),10))

v1 <- list("up0"=top0up, "up1"=top1up, "dn0"=top0dn, "dn1"=top1dn)

plot(euler(v1),quantities = list(cex = 2), labels = list(cex = 2))

Now look at the top pathway results by ES after removing non-significant results (FDR>0.01).

top <- subset(myres,fdr0<0.01)
top <- top[order(top$ES0),]
dn0 <- head(top,15)$ES0
names(dn0) <- rownames(head(top,15))
up0 <- tail(top,15)$ES0
names(up0) <- rownames(tail(top,15))
top0 <- c(dn0,up0)
names(top0) <- gsub("REACTOME_","",names(top0))
names(top0) <- gsub("_"," ",names(top0))

par(mar=c(5.1,20.1,4.1,2.1))
barplot(top0,horiz=TRUE,las=1,cex.names=0.6,main="Uncorr",xlab="ES0")

top <- subset(myres,fdr1<0.01)
top <- top[order(top$ES1),]
dn0 <- head(top,15)$ES1
names(dn0) <- rownames(head(top,15))
up0 <- tail(top,15)$ES1
names(up0) <- rownames(tail(top,15))
top0 <- c(dn0,up0)
names(top0) <- gsub("REACTOME_","",names(top0))
names(top0) <- gsub("_"," ",names(top0))

barplot(top0,horiz=TRUE,las=1,cex.names=0.6,main="Corr",xlab="ES1")

par(mar=c(5.1,4.1,4.1,2.1))

Now plot ES values with both approaches.

plot(myres$ES0,myres$ES1,pch=19,col="gray",cex=0.5)
sig0 <- subset(myres,fdr0<0.01)
sig1 <- subset(myres,fdr1<0.01)
sig2 <- subset(myres,fdr0<0.01 & fdr1<0.01)

points(sig0$ES0,sig0$ES1,pch=19,col="blue",cex=0.5)
points(sig1$ES0,sig1$ES1,pch=19,col="orange",cex=0.5)
points(sig2$ES0,sig2$ES1,pch=19,col="black",cex=0.5)

abline(0,1)

Go back andcheck whether the adjustment reduces bias on basemean gene sets.

RhpcBLASctl::blas_set_num_threads(1)
tic()
myres <- mclapply(1:length(xl), function(i) {
  gs <- xl[[i]]
  gsname <- names(xl)[i]
  df$inset <- as.numeric(rownames(df) %in% gs)
  NGENES <- length(unique(unlist(subset( df,inset==1 )$symbol)))
  mdl0 <- aov(derank ~ inset , data = df)
  p0 <- summary(mdl0)[[1]][1,5]
  mdl1 <- aov(derank ~ bmrank + inset, data = df)
  p1 <- summary(mdl1)[[1]][2,5]
  mrset <- mean(subset(df,inset==1)$derank)
  mrbg <- mean(subset(df,inset!=1)$derank)
  tot <- nrow(df)
  es0 <- (mrset - mrbg) / (tot/2)
  # corrected es
  df2 <- df[order(df$bmrank),]
  idx <- which(df2$inset==1)
  block <- round(nrow(df)/100)
  mrbg1 <- mean(unlist(lapply(idx, function(j) {
    START=j-block
    if (START < 1) { START=1 }
    END=j+block
    if (END > tot) { END=tot }
    df3 <- df2[START:END,]
    df3 <- subset(df3,inset!=1)
    return(df3$derank)
  } )))
  tot <- nrow(df)
  es1 <- (mrset - mrbg1) / (tot/2)
  res <- c("GeneSet"=gsname,"Ngenes"=NGENES,"ES0"=es0,"ES1"=es1,"p0"=p0,"p1"=p1)
  return(res)
},mc.cores=8)
toc()
## 1.313 sec elapsed
myres <- do.call(rbind,myres)
myrownames <- myres[,1]
myres <- myres[,-1]
myres <- apply(myres,2,as.numeric)
rownames(myres) <- myrownames
myres <- as.data.frame(myres)
myres <- subset(myres,Ngenes>=5)

# FDR correction is vital
myres$fdr0 <- p.adjust(myres$p0,method="fdr")
myres$fdr1 <- p.adjust(myres$p1,method="fdr")

# show top results
myres <- myres[order(myres[,"p0"]),]
head(myres,10) %>%
  kbl(caption="Top results for normal ES calc",row.names=TRUE) %>%
  kable_paper("hover", full_width = F)
Top results for normal ES calc
Ngenes ES0 ES1 p0 p1 fdr0 fdr1
level 10 1460 0.2229485 -0.0279211 0.0000000 0.3867808 0.0000000 0.4834760
level 1 1518 -0.2096840 -0.1293504 0.0000000 0.8744022 0.0000000 0.8744022
level 9 1439 0.1998690 0.0165792 0.0000000 0.0111965 0.0000000 0.0279912
level 8 1412 0.1487755 0.0279545 0.0000000 0.0319623 0.0000000 0.0539870
level 3 1511 -0.1132844 -0.0320161 0.0000000 0.7648047 0.0000000 0.8497830
level 2 1512 -0.1114480 -0.0670285 0.0000000 0.0001242 0.0000000 0.0006208
level 5 1498 -0.0972197 0.0065321 0.0000000 0.0000011 0.0000000 0.0000112
level 4 1507 -0.0516235 0.1038482 0.0008779 0.2045079 0.0010974 0.2921542
level 7 1465 0.0386178 -0.0051133 0.0128245 0.0323922 0.0142495 0.0539870
level 6 1477 -0.0269511 -0.0719110 0.0824417 0.0008560 0.0824417 0.0028534 
myres <- myres[order(myres[,"p1"]),]
head(myres,10) %>%
  kbl(caption="Top results for adjusted ES calc",row.names=TRUE) %>%
  kable_paper("hover", full_width = F)
Top results for adjusted ES calc
Ngenes ES0 ES1 p0 p1 fdr0 fdr1
level 5 1498 -0.0972197 0.0065321 0.0000000 0.0000011 0.0000000 0.0000112
level 2 1512 -0.1114480 -0.0670285 0.0000000 0.0001242 0.0000000 0.0006208
level 6 1477 -0.0269511 -0.0719110 0.0824417 0.0008560 0.0824417 0.0028534
level 9 1439 0.1998690 0.0165792 0.0000000 0.0111965 0.0000000 0.0279912
level 8 1412 0.1487755 0.0279545 0.0000000 0.0319623 0.0000000 0.0539870
level 7 1465 0.0386178 -0.0051133 0.0128245 0.0323922 0.0142495 0.0539870
level 4 1507 -0.0516235 0.1038482 0.0008779 0.2045079 0.0010974 0.2921542
level 10 1460 0.2229485 -0.0279211 0.0000000 0.3867808 0.0000000 0.4834760
level 3 1511 -0.1132844 -0.0320161 0.0000000 0.7648047 0.0000000 0.8497830
level 1 1518 -0.2096840 -0.1293504 0.0000000 0.8744022 0.0000000 0.8744022

Plot the result.

ora1bias <- sum(-log10(ora_up$pval)) + sum(-log10(ora_dn$pval))
ora2bias <- sum(-log10(ora2_up$pval)) + sum(-log10(ora2_dn$pval))
fcs1bias <- sum(-log10(fcs1$pval))
anova0bias <- sum(-log10(myres$p0))
anova1bias <- sum(-log10(myres$p1))

res_p <- c("ORA p"=ora1bias,"ORA p FC"=ora2bias,"fgsea"=fcs1bias,"ANOVA"=anova0bias, "ANOVA adj"=anova1bias)
res_p
##     ORA p  ORA p FC     fgsea     ANOVA ANOVA adj 
##  74.02953  74.02953 165.56530 186.21954  19.13495
barplot(res_p, ylab="Bias measure (p)")

ora1bias <- sum(-log10(ora_up$padj)) + sum(-log10(ora_dn$padj))
ora2bias <- sum(-log10(ora2_up$padj)) + sum(-log10(ora2_dn$padj))
fcs1bias <- sum(-log10(fcs1$padj))
anova0bias <- sum(-log10(myres$fdr0))
anova1bias <- sum(-log10(myres$fdr1))

res_fdr <- c("ORA p"=ora1bias,"ORA p FC"=ora2bias,"fgsea"=fcs1bias,"ANOVA"=anova0bias, "ANOVA adj"=anova1bias)
res_fdr
##     ORA p  ORA p FC     fgsea     ANOVA ANOVA adj 
##  67.74499  67.74499 162.12506 182.77930  15.76809
barplot(res_fdr, ylab="Bias measure (FDR)")

Session info

For reproducibility.

sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 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_AU.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8    
##  [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8   
##  [7] LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] tictoc_1.2.1     eulerr_7.0.4     beeswarm_0.4.0   kableExtra_1.4.0
## [5] fgsea_1.34.2    
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10         generics_0.1.4      xml2_1.5.2         
##  [4] polylabelr_1.0.0    stringi_1.8.7       lattice_0.22-7     
##  [7] digest_0.6.39       magrittr_2.0.4      evaluate_1.0.5     
## [10] grid_4.5.2          RColorBrewer_1.1-3  fastmap_1.2.0      
## [13] jsonlite_2.0.0      Matrix_1.7-4        viridisLite_0.4.2  
## [16] scales_1.4.0        RhpcBLASctl_0.23-42 codetools_0.2-20   
## [19] textshaping_1.0.4   jquerylib_0.1.4     cli_3.6.5          
## [22] rlang_1.1.7         polyclip_1.10-7     cowplot_1.2.0      
## [25] cachem_1.1.0        yaml_2.3.12         otel_0.2.0         
## [28] tools_4.5.2         BiocParallel_1.42.2 dplyr_1.1.4        
## [31] ggplot2_4.0.1       fastmatch_1.1-8     vctrs_0.7.0        
## [34] R6_2.6.1            lifecycle_1.0.5     stringr_1.6.0      
## [37] pkgconfig_2.0.3     pillar_1.11.1       bslib_0.9.0        
## [40] gtable_0.3.6        glue_1.8.0          data.table_1.18.0  
## [43] Rcpp_1.1.1          systemfonts_1.3.1   xfun_0.56          
## [46] tibble_3.3.1        tidyselect_1.2.1    rstudioapi_0.18.0  
## [49] knitr_1.51          dichromat_2.0-0.1   farver_2.1.2       
## [52] htmltools_0.5.9     rmarkdown_2.30      svglite_2.2.2      
## [55] compiler_4.5.2      S7_0.2.1