Source: TBA

Introduction

Here we analyse the effect of VLCAD mutation.

Pathway enrichment analysis is conducted with mitch.

Gene sets were obtained from Reactome (14th March 2025) and GO Biological process sets from MSigDB version 2024.1. Here are the file names: gobp.msigdb.v2024.1.Hs.symbols.gmt and ReactomePathways_2025-03-14.gmt.

suppressPackageStartupMessages({
    library("mitch")
    library("kableExtra")
})

knitr::opts_chunk$set(dev = 'svg') # set output device to svg

Import proteomics data

dat <- read.table("Adjusted_proteomics_mito_repeat.tsv",sep="\t",header=TRUE,fill=NA)

head(dat) %>% kbl(caption = "Input data") %>% kable_paper("hover", full_width = F)
Input data
X AVG.Group.Quantity.Denominator AVG.Group.Quantity.Numerator AVG.Log2.Ratio Fold.change Absolute.AVG.Log2.Ratio X..Change Ratio Pvalue Qvalue
GLS2;>NP_001267726 7316.844 10324.63 0.4967971 0.4967971 41.1077394 1.411077 0.0140768 0.0034774 4.0851808
FASTK;>NP_006703 10227.585 10968.56 0.1009082 0.1009082 7.2448352 1.072448 0.1801507 0.0329853 5.5311232
MSRB3;>NP_001180389 5250.715 12661.83 1.2698997 1.2698997 141.1448069 2.411448 0.0009228 0.0003594 4.4669652
TRIAP1 14461.727 13509.39 -0.0982774 0.9341477 0.0982774 -6.585229 0.9341477 0.7604954 0.1214302
D2HGDH 50164.922 13601.22 -1.8829428 0.2711301 1.8829428 -72.886989 0.2711301 0.0000195 0.0000212
AKR1B10 5461.549 14165.92 1.3750426 2.5937557 1.3750426 159.375573 2.5937557 0.0014805 0.0005132

Extract the key data for scoring

gene <- sapply(strsplit(dat$X,";"),"[[",1)

lfc <- log2(dat$Fold.change)

logp <- -log10(dat$Pvalue)

score <- sign(lfc) / logp

df <- data.frame("gene"=gene,"score"=score)

head(df) %>% kbl(caption = "Cleaned data") %>% kable_paper("hover", full_width = F)
Cleaned data
gene score
GLS2 -0.4067125
FASTK -0.6749100
MSRB3 0.2903246
TRIAP1 -8.4101869
D2HGDH -0.2123158
AKR1B10 0.3534066

Data checks

Running some checks of the data.

There are some duplicated genes.

str(df)
## 'data.frame':    1015 obs. of  2 variables:
##  $ gene : chr  "GLS2" "FASTK" "MSRB3" "TRIAP1" ...
##  $ score: num  -0.407 -0.675 0.29 -8.41 -0.212 ...
dfa <- aggregate(score ~ gene,df,mean)

rownames(dfa) <- dfa$gene

dfa$gene = NULL

length(unique(df$gene))
## [1] 985
df$gene[which(duplicated(df$gene))]
##  [1] "D2HGDH"   "PNKD"     "PTPMT1"   "UNG"      "MFF"      "HTATIP2" 
##  [7] "DMAC2L"   "NDUFV3"   "ECHDC1"   "COX11"    "MTX3"     "TMEM126A"
## [13] "MRPS35"   "SDHC"     "HTRA2"    "NDUFB5"   "IDH3G"    "ALDH7A1" 
## [19] "KARS1"    "YME1L1"   "DUT"      "PGAM5"    "GLS"      "IMMT"    
## [25] "DMAC2"    "GPX4"     "GSTK1"    "MRPL43"   "SLC25A10" "TMEM126B"
str(dfa)
## 'data.frame':    982 obs. of  1 variable:
##  $ score: num  -7.246 -0.226 -0.782 0.185 0.234 ...

Histogram of p-values looks okay.

hist(dat$Pvalue)

Histogram of lfcs looks okay. Long left tail. Slightly more genes in the down direction.

hist(lfc)

table(sign(lfc))
## 
##  -1   1 
## 528 485

Pathway analysis

Reactome is good for getting a general idea about the canonical pathways.

reactome <- gmt_import("ReactomePathways_2025-03-14.gmt")

mr <- mitch_calc(dfa, reactome, priority="significance",cores=8,minsetsize=5)
## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.
mrtop <- subset(mr$enrichment_result,p.adjustANOVA < 0.05 )

head(mrtop,20) %>%
  kbl(caption = "Top gene pathway differences caused by ALVCAD mutation") %>%
  kable_paper("hover", full_width = F)
Top gene pathway differences caused by ALVCAD mutation
set setSize pANOVA s.dist p.adjustANOVA
93 Mitochondrial translation 93 0.0000214 -0.2665433 0.0025969
96 Mitochondrial translation termination 87 0.0000431 -0.2644192 0.0025969
95 Mitochondrial translation initiation 87 0.0000464 -0.2633147 0.0025969
94 Mitochondrial translation elongation 87 0.0000701 -0.2570731 0.0029455
148 Translation 113 0.0002111 -0.2134587 0.0070916
80 Metabolism of proteins 234 0.0006763 -0.1467446 0.0189367
168 tRNA processing in the mitochondrion 17 0.0013364 0.4525450 0.0320730 
write.table(mr$enrichment_result,file="mitch_reactome_june2026.tsv",quote=FALSE,sep='\t')

par(mar=c(5,27,3,3))
top <- mr$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)
## [1] 7
up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.65,main="Ctl vs Mut",xlab="ES")
grid()

if ( ! file.exists("mitch_reactome_june2026.html") ) {
  mitch_report(res=mr,outfile="mitch_reactome_june2026.html",overwrite=TRUE)
}

mitch::networkplot(mr)
## Can't plot upregulated sets. Fewer than 5 found.

Gene ontology is a bigger database and sometimes it yields more detailed results.

gobp <- gmt_import("gobp.msigdb.v2024.1.Hs.symbols.gmt")
names(gobp) <- gsub("GOBP_","",names(gobp))
names(gobp) <- gsub("_"," ",names(gobp))

#download.file("https://ziemann-lab.net/public/tmp/go_2024-11.gmt",destfile="go_2024-11.gmt")
#gobp <- gmt_import("go_2024-11.gmt")

mgo <- mitch_calc(dfa, gobp, priority="significance",cores=8, minsetsize=5)
## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.
mgo_top <- subset(mgo$enrichment_result,p.adjustANOVA < 0.05 )

head(mgo_top,20) %>%
  kbl(caption = "Top GOBP differences caused by ALVCAD mutation") %>%
  kable_paper("hover", full_width = F)
Top GOBP differences caused by ALVCAD mutation
set setSize pANOVA s.dist p.adjustANOVA
403 MITOCHONDRIAL GENE EXPRESSION 156 2.94e-05 -0.2099010 0.0325504
419 MITOCHONDRIAL TRANSLATION 128 7.89e-05 -0.2154384 0.0436519 
write.table(mgo$enrichment_result,file="mitch_gobp_june2026.tsv",quote=FALSE,sep='\t')

par(mar=c(5,27,3,3))
top <- mgo$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)
## [1] 2
up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.65,main="Ctl vs Mut",xlab="ES")
grid()

if ( ! file.exists("mitch_gobp_june2026.html") ) {
  mitch_report(res=mgo,outfile="mitch_gobp_june2026.html")
}

Session information

sessionInfo()
## R version 4.6.0 (2026-04-24)
## 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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] gtools_3.9.5     kableExtra_1.4.0 mitch_1.24.0    
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6          beeswarm_0.4.0        bslib_0.11.0         
##  [4] xfun_0.57             ggplot2_4.0.3         htmlwidgets_1.6.4    
##  [7] caTools_1.18.3        GGally_2.4.0          lattice_0.22-9       
## [10] vctrs_0.7.3           tools_4.6.0           bitops_1.0-9         
## [13] generics_0.1.4        parallel_4.6.0        tibble_3.3.1         
## [16] pkgconfig_2.0.3       KernSmooth_2.23-26    RColorBrewer_1.1-3   
## [19] S7_0.2.2              lifecycle_1.0.5       compiler_4.6.0       
## [22] farver_2.1.2          stringr_1.6.0         textshaping_1.0.5    
## [25] gplots_3.3.0          httpuv_1.6.17         sass_0.4.10          
## [28] htmltools_0.5.9       yaml_2.3.12           jquerylib_0.1.4      
## [31] later_1.4.8           pillar_1.11.1         tidyr_1.3.2          
## [34] MASS_7.3-65           cachem_1.1.0          mime_0.13            
## [37] ggstats_0.13.0        network_1.20.0        tidyselect_1.2.1     
## [40] digest_0.6.39         stringi_1.8.7         dplyr_1.2.1          
## [43] reshape2_1.4.5        purrr_1.2.2           fastmap_1.2.0        
## [46] grid_4.6.0            cli_3.6.6             magrittr_2.0.5       
## [49] dichromat_2.0-0.1     withr_3.0.2           scales_1.4.0         
## [52] promises_1.5.0        rmarkdown_2.31        otel_0.2.0           
## [55] gridExtra_2.3         coda_0.19-4.1         shiny_1.13.0         
## [58] evaluate_1.0.5        knitr_1.51            viridisLite_0.4.3    
## [61] rlang_1.2.0           Rcpp_1.1.1-1.1        xtable_1.8-8         
## [64] glue_1.8.1            echarts4r_0.5.0       xml2_1.5.2           
## [67] jsonlite_2.0.0        svglite_2.2.2         rstudioapi_0.18.0    
## [70] R6_2.6.1              plyr_1.8.9            statnet.common_4.13.0
## [73] systemfonts_1.3.2
save.image("pathway_analysis.Rdata")