Introduction

suppressPackageStartupMessages({
    library("plyr")
    library("RhpcBLASctl")
    library("gplots")
    library("mitch")
    library("eulerr")
    library("limma")
    library("beeswarm")
    library("kableExtra")
})

RhpcBLASctl::blas_set_num_threads(1)

Import gene sets

We’ll fetch reactome and GO.

go <- gmt_import("c5.go.v2025.1.Hs.symbols.gmt")
names(go) <- lapply(names(go),function(n) { gsub("_"," ",n) } )

reactome <- gmt_import("ReactomePathways_2025-10-13.gmt")

head(go)
## $`GOBP 10 FORMYLTETRAHYDROFOLATE METABOLIC PROCESS`
## [1] "AASDHPPT" "ALDH1L1"  "ALDH1L2"  "MTHFD1"   "MTHFD1L"  "MTHFD2L" 
## 
## $`GOBP 2FE 2S CLUSTER ASSEMBLY`
##  [1] "BOLA2"   "BOLA2B"  "FDX2"    "FXN"     "GLRX3"   "GLRX5"   "HSCB"   
##  [8] "ISCU"    "LYRM4"   "NDUFAB1" "NFS1"   
## 
## $`GOBP 2 OXOGLUTARATE METABOLIC PROCESS`
##  [1] "AADAT"  "ADHFE1" "COL6A1" "D2HGDH" "DLD"    "DLST"   "GOT1"   "GOT2"  
##  [9] "GPT2"   "IDH1"   "IDH2"   "KGD4"   "KYAT3"  "L2HGDH" "OGDH"   "OGDHL" 
## [17] "PHYH"   "TAT"   
## 
## $`GOBP 3 UTR MEDIATED MRNA DESTABILIZATION`
##  [1] "CPEB3"   "DHX36"   "DND1"    "HNRNPD"  "KHSRP"   "MOV10"   "PLEKHN1"
##  [8] "PUM1"    "PUM2"    "RBM24"   "RC3H1"   "TARDBP"  "TRIM71"  "UPF1"   
## [15] "ZC3H12A" "ZC3H12D" "ZFP36"   "ZFP36L1" "ZFP36L2"
## 
## $`GOBP 3 UTR MEDIATED MRNA STABILIZATION`
##  [1] "ANGEL2"   "ARID5A"   "BOLL"     "DAZ1"     "DAZ2"     "DAZ3"    
##  [7] "DAZ4"     "DAZL"     "ELAVL1"   "ELAVL4"   "HNRNPA0"  "HNRNPC"  
## [13] "LARP4B"   "MAPK14"   "MAPKAPK2" "METTL16"  "MYD88"    "NICOL1"  
## [19] "QKI"      "RBM10"    "RBM24"    "RBM38"    "RBM47"    "TARDBP"  
## [25] "TENT4A"   "TENT4B"   "TIRAP"    "YBX3"     "ZFP36"   
## 
## $`GOBP 4FE 4S CLUSTER ASSEMBLY`
## [1] "FDX2"  "FXN"   "ISCU"  "LYRM4" "NFS1"
length(go)
## [1] 10480
length(reactome)
## [1] 2810

Import mouse to human gene mapping

Make a gene table file from mouse stable IDs to current gene names.

mm2hs <- read.table("mart_export.txt",fill=NA,header=TRUE,sep="\t")
gt <- unique(mm2hs[,c(1,3)])
gt <- gt[which(gt$Human.gene.name != ""),]
dim(gt)
## [1] 25102     2
head(gt)
##        Gene.stable.ID Human.gene.name
## 6  ENSMUSG00000064341          MT-ND1
## 10 ENSMUSG00000064345          MT-ND2
## 16 ENSMUSG00000064351          MT-CO1
## 19 ENSMUSG00000064354          MT-CO2
## 21 ENSMUSG00000064356         MT-ATP8
## 22 ENSMUSG00000064357         MT-ATP6

Effect of diabetes on gene expression

Characterise DBDB in these samples (DVvCV).

Distal_tubule_cells.DBDB_CL25vCV.tsv Distal_tubule_cells.DBDB_CL50vCV.tsv Distal_tubule_cells.DBDB_DL25vCL25.tsv Distal_tubule_cells.DBDB_DL25vDV.tsv Distal_tubule_cells.DBDB_DL50vCL50.tsv Distal_tubule_cells.DBDB_DL50vDV.tsv Distal_tubule_cells.DBDB_DVvCV.tsv

  1. “Distal_tubule_cells.DBDB_DVvCV.tsv”
  2. “Endothelial_cells.DBDB_DVvCV.tsv”
  3. “Erythroid-like_and_erythroid_precursor_cells.DBDB_DVvCV.tsv”
  4. “Intercalated_cells.DBDB_DVvCV.tsv”
  5. “LOH.DBDB_DVvCV.tsv”
  6. “Macrophages.DBDB_DVvCV.tsv”
  7. “Mesangial_cells.DBDB_DVvCV.tsv”
  8. “Neutrophils.DBDB_DVvCV.tsv”
  9. “NK_cells.DBDB_DVvCV.tsv”
  10. “Podocytes.DBDB_DVvCV.tsv”
  11. “Principal_cells.DBDB_DVvCV.tsv”
  12. “Proximal_tubule_cells.DBDB_DVvCV.tsv”
d <- read.table("Distal_tubule_cells.DBDB_DVvCV.tsv")
d <- d[order(d$p_val),]
g <- -log10(d$p_val) * sign(d$avg_logFC)
names(g) <- sapply( strsplit(rownames(d),"_"), "[[", 1  )
h <- mitch_import(x=as.data.frame(g),DEtype="preranked",geneTable=gt)
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1612
## Note: no. genes in output = 1521
## Note: estimated proportion of input genes in output = 0.944
m <- mitch_calc(x=h,genesets=reactome,minsetsize=5,cores=8,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
subset(m$enrichment_result,p.adjustANOVA<0.05)
##           set setSize       pANOVA     s.dist p.adjustANOVA
## 162 Drug ADME      17 1.214465e-07 -0.7422951  8.914171e-05

Now analyse all cell types for diabetes.

myfiles <- list.files(".",pattern="DBDB_DVvCV.tsv")

dbdb_res0 <- lapply(myfiles, function(myfile) {
  d <- read.table(myfile)
  d <- d[order(d$p_val),]
  subset(d,p_val_adj<0.05)
})
names(dbdb_res0) <- myfiles
lapply(dbdb_res0,nrow)
## $Distal_tubule_cells.DBDB_DVvCV.tsv
## [1] 10
## 
## $Endothelial_cells.DBDB_DVvCV.tsv
## [1] 34
## 
## $`Erythroid-like_and_erythroid_precursor_cells.DBDB_DVvCV.tsv`
## [1] 0
## 
## $Intercalated_cells.DBDB_DVvCV.tsv
## [1] 0
## 
## $LOH.DBDB_DVvCV.tsv
## [1] 1
## 
## $Macrophages.DBDB_DVvCV.tsv
## [1] 1
## 
## $Mesangial_cells.DBDB_DVvCV.tsv
## [1] 0
## 
## $Neutrophils.DBDB_DVvCV.tsv
## [1] 0
## 
## $NK_cells.DBDB_DVvCV.tsv
## [1] 0
## 
## $Podocytes.DBDB_DVvCV.tsv
## [1] 4
## 
## $Principal_cells.DBDB_DVvCV.tsv
## [1] 0
## 
## $Proximal_tubule_cells.DBDB_DVvCV.tsv
## [1] 1014
names(dbdb_res0) <- gsub(".DBDB_DVvCV.tsv","",names(dbdb_res0))

nsig <- unlist(lapply(dbdb_res0,nrow))
nsig <- sort(nsig)

par(mar=c(5.1,20.1,4.1,2.1)) ; barplot(nsig,horiz=TRUE,las=1,xlim=c(0,1200))
text(nsig+100,(1:length(nsig)*1.18)-0.3,labels=nsig)
mtext("No. significant genes")

lapply(dbdb_res0,head,10)
## $Distal_tubule_cells
##                                   p_val  avg_logFC pct.1 pct.2    p_val_adj
## ENSMUSG00000032758_Kap     3.229231e-29 -1.9368299 0.690 0.986 8.140247e-25
## ENSMUSG00000028713_Cyp4b1  6.832508e-14 -1.0696966 0.388 0.797 1.722339e-09
## ENSMUSG00000024866_Acy3    8.568379e-13 -1.1276846 0.302 0.703 2.159917e-08
## ENSMUSG00000041052_Slc7a13 5.484520e-10 -0.6500861 0.008 0.297 1.382538e-05
## ENSMUSG00000018339_Gpx3    8.641442e-10  0.3373351 1.000 1.000 2.178335e-05
## ENSMUSG00000031958_Ldhd    1.110887e-07 -0.7508959 0.202 0.554 2.800324e-03
## ENSMUSG00000022037_Clu     1.333308e-07 -0.9357106 0.884 0.919 3.361003e-03
## ENSMUSG00000074207_Adh1    8.403786e-07 -0.5739375 0.047 0.297 2.118426e-02
## ENSMUSG00000015656_Hspa8   1.301640e-06 -0.4860473 0.829 0.919 3.281175e-02
## ENSMUSG00000025479_Cyp2e1  1.582278e-06 -0.7130597 0.202 0.514 3.988607e-02
##                              MusID
## ENSMUSG00000032758_Kap         Kap
## ENSMUSG00000028713_Cyp4b1   Cyp4b1
## ENSMUSG00000024866_Acy3       Acy3
## ENSMUSG00000041052_Slc7a13 Slc7a13
## ENSMUSG00000018339_Gpx3       Gpx3
## ENSMUSG00000031958_Ldhd       Ldhd
## ENSMUSG00000022037_Clu         Clu
## ENSMUSG00000074207_Adh1       Adh1
## ENSMUSG00000015656_Hspa8     Hspa8
## ENSMUSG00000025479_Cyp2e1   Cyp2e1
## 
## $Endothelial_cells
##                                    p_val  avg_logFC pct.1 pct.2    p_val_adj
## ENSMUSG00000032758_Kap      1.497865e-37 -1.8739494 0.627 0.981 3.775817e-33
## ENSMUSG00000023944_Hsp90ab1 1.561473e-20 -0.9442815 0.854 0.942 3.936161e-16
## ENSMUSG00000090877_Hspa1b   1.675375e-19 -2.3889747 0.032 0.490 4.223285e-15
## ENSMUSG00000091971_Hspa1a   2.679158e-18 -2.1146852 0.006 0.423 6.753621e-14
## ENSMUSG00000028410_Dnaja1   3.719704e-14 -1.1828842 0.418 0.769 9.376631e-10
## ENSMUSG00000064339_mt-Rnr2  1.721192e-11  0.3090577 1.000 1.000 4.338781e-07
## ENSMUSG00000018339_Gpx3     4.510087e-11  0.2364702 1.000 0.990 1.136903e-06
## ENSMUSG00000029304_Spp1     1.563092e-10  0.9865497 0.437 0.058 3.940243e-06
## ENSMUSG00000063903_Klk1     2.149393e-10  0.9649652 0.570 0.173 5.418189e-06
## ENSMUSG00000059824_Dbp      4.609871e-10  0.9622835 0.405 0.048 1.162056e-05
##                                MusID
## ENSMUSG00000032758_Kap           Kap
## ENSMUSG00000023944_Hsp90ab1 Hsp90ab1
## ENSMUSG00000090877_Hspa1b     Hspa1b
## ENSMUSG00000091971_Hspa1a     Hspa1a
## ENSMUSG00000028410_Dnaja1     Dnaja1
## ENSMUSG00000064339_mt-Rnr2   mt-Rnr2
## ENSMUSG00000018339_Gpx3         Gpx3
## ENSMUSG00000029304_Spp1         Spp1
## ENSMUSG00000063903_Klk1         Klk1
## ENSMUSG00000059824_Dbp           Dbp
## 
## $`Erythroid-like_and_erythroid_precursor_cells`
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Intercalated_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $LOH
##                               p_val avg_logFC pct.1 pct.2    p_val_adj MusID
## ENSMUSG00000032758_Kap 1.238408e-10 -1.733313 0.853     1 3.121778e-06   Kap
## 
## $Macrophages
##                               p_val avg_logFC pct.1 pct.2   p_val_adj MusID
## ENSMUSG00000032758_Kap 4.149869e-08 -2.399468  0.56     1 0.001046099   Kap
## 
## $Mesangial_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Neutrophils
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $NK_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Podocytes
##                                    p_val avg_logFC pct.1 pct.2    p_val_adj
## ENSMUSG00000032758_Kap      5.013010e-12 -2.121591 0.500 1.000 1.263680e-07
## ENSMUSG00000023944_Hsp90ab1 2.736756e-09 -1.153839 1.000 0.983 6.898815e-05
## ENSMUSG00000015656_Hspa8    2.756868e-07 -1.057190 0.667 0.931 6.949512e-03
## ENSMUSG00000037573_Tob1     1.540323e-06 -1.303275 0.292 0.845 3.882847e-02
##                                MusID
## ENSMUSG00000032758_Kap           Kap
## ENSMUSG00000023944_Hsp90ab1 Hsp90ab1
## ENSMUSG00000015656_Hspa8       Hspa8
## ENSMUSG00000037573_Tob1         Tob1
## 
## $Principal_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Proximal_tubule_cells
##                                     p_val  avg_logFC pct.1 pct.2     p_val_adj
## ENSMUSG00000032758_Kap       0.000000e+00 -2.0059378 0.719 0.999  0.000000e+00
## ENSMUSG00000029304_Spp1      0.000000e+00  1.5413347 0.815 0.255  0.000000e+00
## ENSMUSG00000024866_Acy3     1.292677e-300 -1.3425344 0.560 0.886 3.258579e-296
## ENSMUSG00000064373_Selenop  1.029043e-249 -1.0008097 0.739 0.920 2.594011e-245
## ENSMUSG00000028713_Cyp4b1   4.286158e-243 -1.1748531 0.763 0.932 1.080455e-238
## ENSMUSG00000024644_Cndp2    1.223609e-207 -1.3744239 0.338 0.724 3.084473e-203
## ENSMUSG00000023944_Hsp90ab1 3.352546e-198 -0.6982340 0.847 0.939 8.451099e-194
## ENSMUSG00000030945_Acsm2    1.246829e-176 -0.4920636 0.978 0.991 3.143007e-172
## ENSMUSG00000061906_Ugt2b38  4.058400e-171 -1.0268020 0.597 0.829 1.023042e-166
## ENSMUSG00000025479_Cyp2e1   1.564145e-156 -1.0977585 0.551 0.805 3.942897e-152
##                                MusID
## ENSMUSG00000032758_Kap           Kap
## ENSMUSG00000029304_Spp1         Spp1
## ENSMUSG00000024866_Acy3         Acy3
## ENSMUSG00000064373_Selenop   Selenop
## ENSMUSG00000028713_Cyp4b1     Cyp4b1
## ENSMUSG00000024644_Cndp2       Cndp2
## ENSMUSG00000023944_Hsp90ab1 Hsp90ab1
## ENSMUSG00000030945_Acsm2       Acsm2
## ENSMUSG00000061906_Ugt2b38   Ugt2b38
## ENSMUSG00000025479_Cyp2e1     Cyp2e1

In all subsequent work, just focus on:

  1. “Distal_tubule_cells.DBDB_DVvCV.tsv”
  2. “Endothelial_cells.DBDB_DVvCV.tsv”
  3. “Podocytes.DBDB_DVvCV.tsv”
  4. “Proximal_tubule_cells.DBDB_DVvCV.tsv”

GO enrichment for diabetes

dbdb_res_go <- lapply(myfiles, function(myfile) {
  d <- read.table(myfile)
  d <- d[order(d$p_val),]
  g <- -log10(d$p_val) * sign(d$avg_logFC)
  names(g) <- sapply( strsplit(rownames(d),"_"), "[[", 1  )
  h <- mitch_import(x=as.data.frame(g),DEtype="preranked",geneTable=gt)
  m <- mitch_calc(x=h,genesets=go,minsetsize=5,cores=8,priority="effect")
  subset(m$enrichment_result,p.adjustANOVA<0.05)
})
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1612
## Note: no. genes in output = 1521
## Note: estimated proportion of input genes in output = 0.944
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1804
## Note: no. genes in output = 1746
## Note: estimated proportion of input genes in output = 0.968
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1984
## Note: no. genes in output = 1872
## Note: estimated proportion of input genes in output = 0.944
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1797
## Note: no. genes in output = 1620
## Note: estimated proportion of input genes in output = 0.902
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1756
## Note: no. genes in output = 1651
## Note: estimated proportion of input genes in output = 0.94
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1359
## Note: no. genes in output = 1310
## Note: estimated proportion of input genes in output = 0.964
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 605
## Note: no. genes in output = 574
## Note: estimated proportion of input genes in output = 0.949
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 223
## Note: no. genes in output = 218
## Note: estimated proportion of input genes in output = 0.978
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 440
## Note: no. genes in output = 401
## Note: estimated proportion of input genes in output = 0.911
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1684
## Note: no. genes in output = 1581
## Note: estimated proportion of input genes in output = 0.939
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1535
## Note: no. genes in output = 1405
## Note: estimated proportion of input genes in output = 0.915
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 3147
## Note: no. genes in output = 3005
## Note: estimated proportion of input genes in output = 0.955
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
names(dbdb_res_go) <- myfiles
names(dbdb_res_go) <- gsub(".DBDB_DVvCV.tsv","",names(dbdb_res_go))
lapply(dbdb_res_go,nrow)
## $Distal_tubule_cells
## [1] 22
## 
## $Endothelial_cells
## [1] 37
## 
## $`Erythroid-like_and_erythroid_precursor_cells`
## [1] 68
## 
## $Intercalated_cells
## [1] 1
## 
## $LOH
## [1] 0
## 
## $Macrophages
## [1] 5
## 
## $Mesangial_cells
## [1] 7
## 
## $Neutrophils
## [1] 0
## 
## $NK_cells
## [1] 3
## 
## $Podocytes
## [1] 8
## 
## $Principal_cells
## [1] 0
## 
## $Proximal_tubule_cells
## [1] 266
nsig <- unlist(lapply(dbdb_res_go,nrow))
nsig <- sort(nsig)
nsig
##                                          LOH 
##                                            0 
##                                  Neutrophils 
##                                            0 
##                              Principal_cells 
##                                            0 
##                           Intercalated_cells 
##                                            1 
##                                     NK_cells 
##                                            3 
##                                  Macrophages 
##                                            5 
##                              Mesangial_cells 
##                                            7 
##                                    Podocytes 
##                                            8 
##                          Distal_tubule_cells 
##                                           22 
##                            Endothelial_cells 
##                                           37 
## Erythroid-like_and_erythroid_precursor_cells 
##                                           68 
##                        Proximal_tubule_cells 
##                                          266
par(mar=c(5.1,20.1,4.1,2.1)) ; barplot(nsig,horiz=TRUE,las=1,xlim=c(0,300))
text(nsig+20,(1:length(nsig)*1.18)-0.3,labels=nsig)
mtext("No. significant GO terms")

x <- dbdb_res_go$Proximal_tubule_cells
up <- head(subset(x,s.dist>0),20)
ups <- up$s.dist
names(ups) <- up$set
dn <- head(subset(x,s.dist<0),20)
dns <- dn$s.dist
names(dns) <- dn$set
dns <- rev(dns)
pws <- c(ups,dns)
pws <- sort(pws)
par(mar=c(5.1,25.1,4.1,2.1))
barplot(pws,horiz=TRUE,las=1,cex.names=0.6,xlab="ES")
mtext("Proximal tubule GO terms")

mypws <- names(which(abs(pws)>0.4))

x <- dbdb_res_go$Endothelial_cells
up <- head(subset(x,s.dist>0),20)
ups <- up$s.dist
names(ups) <- up$set
dn <- head(subset(x,s.dist<0),20)
dns <- dn$s.dist
names(dns) <- dn$set
dns <- rev(dns)
pws <- c(ups,dns)
pws <- sort(pws)
par(mar=c(5.1,25.1,4.1,2.1))
barplot(pws,horiz=TRUE,las=1,cex.names=0.7,xlab="ES")
mtext("Endothelial cell GO terms")

mypws <- c(mypws,names(which(abs(pws)>0.4)))

x <- dbdb_res_go$Podocytes
up <- head(subset(x,s.dist>0),20)
ups <- up$s.dist
names(ups) <- up$set
dn <- head(subset(x,s.dist<0),20)
dns <- dn$s.dist
names(dns) <- dn$set
dns <- rev(dns)
pws <- c(ups,dns)
pws <- sort(pws)
par(mar=c(5.1,25.1,4.1,2.1))
barplot(pws,horiz=TRUE,las=1,cex.names=0.7,xlab="ES")
mtext("Podocyte GO terms")

mypws <- c(mypws,names(which(abs(pws)>0.4)))

x <- dbdb_res_go$Distal_tubule_cells
up <- head(subset(x,s.dist>0),20)
ups <- up$s.dist
names(ups) <- up$set
dn <- head(subset(x,s.dist<0),20)
dns <- dn$s.dist
names(dns) <- dn$set
dns <- rev(dns)
pws <- c(ups,dns)
pws <- sort(pws)
par(mar=c(5.1,25.1,4.1,2.1))
barplot(pws,horiz=TRUE,las=1,cex.names=0.7,xlab="ES")
mtext("Distal tubule cell GO terms")

mypws <- c(mypws,names(which(abs(pws)>0.4)))

mypws <- unique(mypws)

prox <- dbdb_res_go$Proximal_tubule_cells[dbdb_res_go$Proximal_tubule_cells$set %in% mypws,c("set","s.dist")]
colnames(prox) <- c("set","prox")
endo <- dbdb_res_go$Endothelial_cells[dbdb_res_go$Endothelial_cells$set %in% mypws,c("set","s.dist")]
colnames(endo) <- c("set","endo")
podo <- dbdb_res_go$Podocytes[dbdb_res_go$Podocytes$set %in% mypws,c("set","s.dist")]
colnames(podo) <- c("set","podo")
dstl <- dbdb_res_go$Distal_tubule_cells[dbdb_res_go$Distal_tubule_cells$set %in% mypws,c("set","s.dist")]
colnames(dstl) <- c("set","dstl")

multi3 <- join_all(list(prox,endo,podo,dstl))
## Joining by: set
## Joining by: set
## Joining by: set
rownames(multi3) <- multi3$set
multi3$set=NULL

colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(  as.matrix(multi3), col=colfunc(25),scale="none", trace="none",
  margins = c(6,25), cexRow=0.5, cexCol=0.9, main="Top GO terms")

Deeper look at proximal tubules.

d <- read.table("Proximal_tubule_cells.DBDB_DVvCV.tsv")
d <- d[order(d$p_val),]
g <- -log10(d$p_val) * sign(d$avg_logFC)
names(g) <- sapply( strsplit(rownames(d),"_"), "[[", 1  )
h <- mitch_import(x=as.data.frame(g),DEtype="preranked",geneTable=gt)
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 3147
## Note: no. genes in output = 3005
## Note: estimated proportion of input genes in output = 0.955
m <- mitch_calc(x=h,genesets=go,minsetsize=5,cores=8,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(subset(m$enrichment_result,p.adjustANOVA<0.05 & s.dist > 0 ) ,20) %>%
  kbl(caption="Up-regulated GO terms") %>%
  kable_paper("hover", full_width = F)
Up-regulated GO terms
set setSize pANOVA s.dist p.adjustANOVA
2862 GOBP RENAL SODIUM ION ABSORPTION 7 0.0000435 0.8922615 0.0030220
2863 GOBP RENAL SODIUM ION TRANSPORT 7 0.0000435 0.8922615 0.0030220
4442 GOMF SUGAR TRANSMEMBRANE TRANSPORTER ACTIVITY 5 0.0016141 0.8145333 0.0347108
699 GOBP FAT SOLUBLE VITAMIN METABOLIC PROCESS 7 0.0003331 0.7833794 0.0116497
1871 GOBP POSITIVE REGULATION OF FATTY ACID TRANSPORT 5 0.0027407 0.7738000 0.0476047
1321 GOBP NEGATIVE REGULATION OF CYTOKINE PRODUCTION INVOLVED IN INFLAMMATORY RESPONSE 7 0.0008051 0.7317259 0.0214040
1514 GOBP NEGATIVE REGULATION OF T CELL PROLIFERATION 6 0.0020719 0.7263532 0.0401481
1389 GOBP NEGATIVE REGULATION OF LIPID STORAGE 6 0.0020884 0.7257975 0.0401481
3699 GOCC POSTSYNAPTIC CYTOSOL 7 0.0013384 0.7004670 0.0309961
463 GOBP CYTOKINE PRODUCTION INVOLVED IN INFLAMMATORY RESPONSE 9 0.0004437 0.6766058 0.0141111
4108 GOMF LIGAND GATED MONOATOMIC CATION CHANNEL ACTIVITY 9 0.0004777 0.6728230 0.0145808
2078 GOBP POSITIVE REGULATION OF TRANSPORTER ACTIVITY 8 0.0022628 0.6238739 0.0426361
3737 GOCC RESPIRATORY CHAIN COMPLEX IV 12 0.0002958 0.6039648 0.0110075
2482 GOBP REGULATION OF FATTY ACID TRANSPORT 9 0.0017494 0.6029892 0.0360024
1387 GOBP NEGATIVE REGULATION OF LIPID LOCALIZATION 11 0.0007845 0.5853525 0.0211166
3484 GOCC CYTOSOLIC SMALL RIBOSOMAL SUBUNIT 19 0.0000107 0.5843410 0.0010593
2909 GOBP RESPONSE TO DIETARY EXCESS 10 0.0015451 0.5787646 0.0338730
4107 GOMF LIGAND GATED CHANNEL ACTIVITY 10 0.0019773 0.5655426 0.0391641
971 GOBP LIPID EXPORT FROM CELL 11 0.0026341 0.5243214 0.0464676
3726 GOCC PROTON TRANSPORTING ATP SYNTHASE COMPLEX 14 0.0007121 0.5232364 0.0197152
# Here are some of the top pathways upregulated in DBDB
m$detailed_sets$`GOBP RENAL SODIUM ION TRANSPORT`
##  CLCNKB  GUCA2B   KLHL3    SGK1 SLC12A3    UMOD    WNK4 
##  1362.5   943.0  1288.0  1450.0  1432.0  1350.0  1257.0
m$detailed_sets$`GOMF SUGAR TRANSMEMBRANE TRANSPORTER ACTIVITY`
## SLC23A1  SLC2A2  SLC2A5 SLC5A10  SLC5A2 
##     931    1129    1401    1028    1420
m$detailed_sets$`GOBP NEGATIVE REGULATION OF CYTOKINE PRODUCTION INVOLVED IN INFLAMMATORY RESPONSE`
## CUEDC2 NFKBIA  PDCD4  PPARA  SIRPA  STAT3  ZFP36 
##   1003   1260    761   1427    621   1026   1300
m$detailed_sets$`GOBP NEGATIVE REGULATION OF T CELL PROLIFERATION`
##    CEBPB     ITCH  MARCHF7  PRKAR1A     PRNP TNFRSF21 
##     1375     1349      946      123     1163     1339
m$detailed_sets$`GOBP CYTOKINE PRODUCTION INVOLVED IN INFLAMMATORY RESPONSE`
## CUEDC2  HIF1A MAPK14 NFKBIA  PDCD4  PPARA  SIRPA  STAT3  ZFP36 
##   1003    966    398   1260    761   1427    621   1026   1300
m$detailed_sets$`GOBP POSITIVE REGULATION OF FATTY ACID TRANSPORT`
##   ACSL1 CYP4A11   FABP3     MIF   P2RX4 
##  1101.0   996.5  1461.0  1439.0   606.0
head(subset(m$enrichment_result,p.adjustANOVA<0.05 & s.dist < 0 ) ,20) %>%
  kbl(caption="Down-regulated GO terms") %>%
  kable_paper("hover", full_width = F)
Down-regulated GO terms
set setSize pANOVA s.dist p.adjustANOVA
757 GOBP GLUTATHIONE CATABOLIC PROCESS 5 0.0003110 -0.9312000 0.0113540
3499 GOCC ENDOPLASMIC RETICULUM CHAPERONE COMPLEX 5 0.0004886 -0.9005333 0.0148091
4458 GOMF THREONINE TYPE PEPTIDASE ACTIVITY 5 0.0005591 -0.8912000 0.0166113
2154 GOBP PROTEIN FOLDING IN ENDOPLASMIC RETICULUM 8 0.0000202 -0.8702870 0.0018051
4229 GOMF OMEGA PEPTIDASE ACTIVITY 5 0.0009787 -0.8514667 0.0247289
3011 GOBP RETINOIC ACID METABOLIC PROCESS 6 0.0004778 -0.8236079 0.0145808
4332 GOMF PROTEIN DISULFIDE ISOMERASE ACTIVITY 7 0.0001618 -0.8235014 0.0071626
1575 GOBP NONRIBOSOMAL PEPTIDE BIOSYNTHETIC PROCESS 8 0.0001415 -0.7772356 0.0067248
456 GOBP CYSTEINE METABOLIC PROCESS 8 0.0002926 -0.7397814 0.0110075
3907 GOMF ATP DEPENDENT PROTEIN FOLDING CHAPERONE 21 0.0000000 -0.7317758 0.0000030
3717 GOCC PROTEASOME CORE COMPLEX ALPHA SUBUNIT COMPLEX 6 0.0028458 -0.7037902 0.0488663
3716 GOCC PROTEASOME CORE COMPLEX 9 0.0002947 -0.6972259 0.0110075
3719 GOCC PROTEASOME REGULATORY PARTICLE LID SUBCOMPLEX 7 0.0015060 -0.6930335 0.0335420
3174 GOBP SULFUR COMPOUND CATABOLIC PROCESS 19 0.0000004 -0.6751507 0.0000770
507 GOBP DE NOVO POST TRANSLATIONAL PROTEIN FOLDING 13 0.0000477 -0.6519951 0.0031749
508 GOBP DE NOVO PROTEIN FOLDING 16 0.0000094 -0.6405361 0.0009836
4026 GOMF FERROUS IRON BINDING 8 0.0017539 -0.6393060 0.0360024
3203 GOBP TELOMERASE RNA LOCALIZATION 10 0.0005484 -0.6316528 0.0164004
4091 GOMF INTRAMOLECULAR OXIDOREDUCTASE ACTIVITY 21 0.0000009 -0.6214413 0.0001672
1688 GOBP PEPTIDE BIOSYNTHETIC PROCESS 9 0.0015433 -0.6100356 0.0338730
m$detailed_sets$`GOBP GLUTATHIONE CATABOLIC PROCESS`
##   CHAC2    GGT1  GGTLC1  GGTLC2  GGTLC3 
## -1242.0 -1485.5 -1485.5 -1485.5 -1485.5
m$detailed_sets$`GOCC ENDOPLASMIC RETICULUM CHAPERONE COMPLEX`
## DNAJB11 HSP90B1   HSPA5    P4HB   PDIA6 
##   -1255   -1471   -1494   -1324   -1410

Effect of LIRA in diabetic mice

Now take a look at effect of LIRA25 in DBDB mice.

  1. “Distal_tubule_cells.DBDB_DL25vDV.tsv”
  2. “Endothelial_cells.DBDB_DL25vDV.tsv”
  3. “Erythroid-like_and_erythroid_precursor_cells.DBDB_DL25vDV.tsv”
  4. “Intercalated_cells.DBDB_DL25vDV.tsv”
  5. “LOH.DBDB_DL25vDV.tsv”
  6. “Macrophages.DBDB_DL25vDV.tsv”
  7. “Mesangial_cells.DBDB_DL25vDV.tsv”
  8. “NK_cells.DBDB_DL25vDV.tsv”
  9. “Podocytes.DBDB_DL25vDV.tsv”
  10. “Principal_cells.DBDB_DL25vDV.tsv”
  11. “Proximal_tubule_cells.DBDB_DL25vDV.tsv”

Lira25 didn’t do much in proximal tubules.

Instead, let’s look at lira50.

myfiles <- list.files(".", pattern="DBDB_DL50vDV.tsv")

dbdb_res2 <- lapply(myfiles, function(myfile) {
  d <- read.table(myfile)
  d <- d[order(d$p_val),]
  subset(d,p_val_adj<0.05)
})
names(dbdb_res2) <- myfiles
lapply(dbdb_res2,nrow)
## $Distal_tubule_cells.DBDB_DL50vDV.tsv
## [1] 0
## 
## $Endothelial_cells.DBDB_DL50vDV.tsv
## [1] 8
## 
## $`Erythroid-like_and_erythroid_precursor_cells.DBDB_DL50vDV.tsv`
## [1] 0
## 
## $Intercalated_cells.DBDB_DL50vDV.tsv
## [1] 0
## 
## $LOH.DBDB_DL50vDV.tsv
## [1] 0
## 
## $Macrophages.DBDB_DL50vDV.tsv
## [1] 0
## 
## $Mesangial_cells.DBDB_DL50vDV.tsv
## [1] 0
## 
## $NK_cells.DBDB_DL50vDV.tsv
## [1] 0
## 
## $Podocytes.DBDB_DL50vDV.tsv
## [1] 0
## 
## $Principal_cells.DBDB_DL50vDV.tsv
## [1] 0
## 
## $Proximal_tubule_cells.DBDB_DL50vDV.tsv
## [1] 47
names(dbdb_res2) <- gsub(".DBDB_DL50vDV.tsv","",names(dbdb_res2))

nsig <- unlist(lapply(dbdb_res2,nrow))
nsig <- sort(nsig)

par(mar=c(5.1,20.1,4.1,2.1)) ; barplot(nsig,horiz=TRUE,las=1,xlim=c(0,60))
text(nsig+10,(1:length(nsig)*1.18)-0.3,labels=nsig)
mtext("No. significant genes")

lapply(dbdb_res2,head,10)
## $Distal_tubule_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Endothelial_cells
##                                  p_val  avg_logFC pct.1 pct.2    p_val_adj
## ENSMUSG00000073411_H2-D1  1.808645e-09  0.3316971 0.989 0.968 4.559231e-05
## ENSMUSG00000018339_Gpx3   9.307544e-09 -0.2152107 0.993 1.000 2.346246e-04
## ENSMUSG00000052305_Hbb-bs 2.248596e-07  0.5323022 0.529 0.278 5.668260e-03
## ENSMUSG00000023034_Nr4a1  2.737746e-07 -0.7545073 0.130 0.329 6.901311e-03
## ENSMUSG00000033769_Exoc6b 2.827786e-07 -0.6182112 0.025 0.158 7.128283e-03
## ENSMUSG00000071866_Ppia   6.730229e-07  0.4125461 0.815 0.646 1.696556e-02
## ENSMUSG00000060802_B2m    1.227371e-06  0.3009533 0.964 0.930 3.093957e-02
## ENSMUSG00000001240_Ramp2  1.370890e-06  0.3950925 0.819 0.665 3.455740e-02
##                            MusID
## ENSMUSG00000073411_H2-D1   H2-D1
## ENSMUSG00000018339_Gpx3     Gpx3
## ENSMUSG00000052305_Hbb-bs Hbb-bs
## ENSMUSG00000023034_Nr4a1   Nr4a1
## ENSMUSG00000033769_Exoc6b Exoc6b
## ENSMUSG00000071866_Ppia     Ppia
## ENSMUSG00000060802_B2m       B2m
## ENSMUSG00000001240_Ramp2   Ramp2
## 
## $`Erythroid-like_and_erythroid_precursor_cells`
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Intercalated_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $LOH
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Macrophages
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Mesangial_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $NK_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Podocytes
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Principal_cells
## [1] p_val     avg_logFC pct.1     pct.2     p_val_adj MusID    
## <0 rows> (or 0-length row.names)
## 
## $Proximal_tubule_cells
##                                   p_val  avg_logFC pct.1 pct.2    p_val_adj
## ENSMUSG00000052305_Hbb-bs  1.558407e-61  0.6332160 0.608 0.371 3.928433e-57
## ENSMUSG00000027199_Gatm    4.663092e-48  0.3276040 0.958 0.906 1.175472e-43
## ENSMUSG00000069919_Hba-a1  5.873930e-28  0.5904470 0.215 0.086 1.480700e-23
## ENSMUSG00000021453_Gadd45g 6.337038e-21 -0.3533638 0.044 0.130 1.597441e-16
## ENSMUSG00000069917_Hba-a2  6.033574e-19  0.4225553 0.160 0.067 1.520943e-14
## ENSMUSG00000029304_Spp1    1.877751e-18 -0.2647368 0.731 0.815 4.733436e-14
## ENSMUSG00000098178_Gm42418 1.687602e-17 -0.3576872 0.685 0.792 4.254107e-13
## ENSMUSG00000022613_Miox    4.954426e-17  0.1767109 0.965 0.937 1.248912e-12
## ENSMUSG00000064341_mt-Nd1  2.803843e-15  0.1993856 0.974 0.964 7.067928e-11
## ENSMUSG00000024190_Dusp1   3.297884e-15 -0.3301572 0.135 0.231 8.313306e-11
##                              MusID
## ENSMUSG00000052305_Hbb-bs   Hbb-bs
## ENSMUSG00000027199_Gatm       Gatm
## ENSMUSG00000069919_Hba-a1   Hba-a1
## ENSMUSG00000021453_Gadd45g Gadd45g
## ENSMUSG00000069917_Hba-a2   Hba-a2
## ENSMUSG00000029304_Spp1       Spp1
## ENSMUSG00000098178_Gm42418 Gm42418
## ENSMUSG00000022613_Miox       Miox
## ENSMUSG00000064341_mt-Nd1   mt-Nd1
## ENSMUSG00000024190_Dusp1     Dusp1

Venn diagram

Proximal tubules:

  • “Proximal_tubule_cells.DBDB_DVvCV.tsv”
  • “Proximal_tubule_cells.DBDB_DL25vDV.tsv”
pd <- read.table("Proximal_tubule_cells.DBDB_DVvCV.tsv")
pd$g <- sapply(strsplit(rownames(pd),"_"),"[[",1)

pl <- read.table("Proximal_tubule_cells.DBDB_DL50vDV.tsv")
pl$g <- sapply(strsplit(rownames(pl),"_"),"[[",1)

pd_up <- unique(subset(pd,p_val_adj<0.05 & avg_logFC > 0 )$MusID)
pd_dn <- unique(subset(pd,p_val_adj<0.05 & avg_logFC < 0 )$MusID)

pl_up <- unique(subset(pl,p_val_adj<0.05 & avg_logFC > 0 )$MusID)
pl_dn <- unique(subset(pl,p_val_adj<0.05 & avg_logFC < 0 )$MusID)

v1 <- list("Diab up"=pd_up, "Diab dn"=pd_dn,
  "Lira up"=pl_up,"Lira dn"=pl_dn)

plot(euler(v1),quantities = TRUE,main="Proximal tubule cells")

message("Genes up-regulated in diabetes that are lowered by lira50")
## Genes up-regulated in diabetes that are lowered by lira50
itx1 <-intersect(pd_up,pl_dn)
itx1
##  [1] "Apoe"    "Dusp1"   "Eif3a"   "Errfi1"  "G6pc"    "Gadd45g" "Gm42418"
##  [8] "mt-Rnr1" "mt-Rnr2" "Sfrp1"   "Slc12a1" "Slc13a1" "Slc5a12" "Spink1" 
## [15] "Spp1"    "Umod"
message("Genes down-regulated in diabetes that are partially restored by lira50")
## Genes down-regulated in diabetes that are partially restored by lira50
itx2 <- intersect(pd_dn,pl_up)
itx2
##  [1] "Akr1a1"  "Cth"     "Eef1a1"  "Fth1"    "Gm6977"  "Gm8797"  "Miox"   
##  [8] "mt-Nd1"  "Prdx1"   "Slc25a5" "Tubb4b"

GO terms.

dbdb_res3 <- lapply(myfiles, function(myfile) {
  d <- read.table(myfile)
  d <- d[order(d$p_val),]
  g <- -log10(d$p_val) * sign(d$avg_logFC)
  names(g) <- sapply( strsplit(rownames(d),"_"), "[[", 1  )
  h <- mitch_import(x=as.data.frame(g),DEtype="preranked",geneTable=gt)
  m <- mitch_calc(x=h,genesets=go,minsetsize=5,cores=8,priority="effect")
  subset(m$enrichment_result,p.adjustANOVA<0.05)
})
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1326
## Note: no. genes in output = 1246
## Note: estimated proportion of input genes in output = 0.94
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1508
## Note: no. genes in output = 1436
## Note: estimated proportion of input genes in output = 0.952
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1419
## Note: no. genes in output = 1346
## Note: estimated proportion of input genes in output = 0.949
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1975
## Note: no. genes in output = 1697
## Note: estimated proportion of input genes in output = 0.859
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1701
## Note: no. genes in output = 1565
## Note: estimated proportion of input genes in output = 0.92
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1273
## Note: no. genes in output = 1171
## Note: estimated proportion of input genes in output = 0.92
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 2266
## Note: no. genes in output = 2036
## Note: estimated proportion of input genes in output = 0.898
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 392
## Note: no. genes in output = 353
## Note: estimated proportion of input genes in output = 0.901
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1539
## Note: no. genes in output = 1447
## Note: estimated proportion of input genes in output = 0.94
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1075
## Note: no. genes in output = 1006
## Note: estimated proportion of input genes in output = 0.936
## Warning in mitch_rank(input_profile): Warning: >60% of genes have the same score. This isn't
##             optimal for rank based enrichment analysis.
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 1648
## Note: no. genes in output = 1550
## Note: estimated proportion of input genes in output = 0.941
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
names(dbdb_res3) <- myfiles
names(dbdb_res3) <- gsub(".DBDB_DL50vDV.tsv","",names(dbdb_res3))
lapply(dbdb_res3,nrow)
## $Distal_tubule_cells
## [1] 0
## 
## $Endothelial_cells
## [1] 1
## 
## $`Erythroid-like_and_erythroid_precursor_cells`
## [1] 0
## 
## $Intercalated_cells
## [1] 1
## 
## $LOH
## [1] 1
## 
## $Macrophages
## [1] 0
## 
## $Mesangial_cells
## [1] 2
## 
## $NK_cells
## [1] 7
## 
## $Podocytes
## [1] 0
## 
## $Principal_cells
## [1] 0
## 
## $Proximal_tubule_cells
## [1] 246
nsig <- unlist(lapply(dbdb_res3,nrow))
nsig <- sort(nsig)
nsig
##                          Distal_tubule_cells 
##                                            0 
## Erythroid-like_and_erythroid_precursor_cells 
##                                            0 
##                                  Macrophages 
##                                            0 
##                                    Podocytes 
##                                            0 
##                              Principal_cells 
##                                            0 
##                            Endothelial_cells 
##                                            1 
##                           Intercalated_cells 
##                                            1 
##                                          LOH 
##                                            1 
##                              Mesangial_cells 
##                                            2 
##                                     NK_cells 
##                                            7 
##                        Proximal_tubule_cells 
##                                          246
par(mar=c(5.1,20.1,4.1,2.1)) ; barplot(nsig,horiz=TRUE,las=1,xlim=c(0,300))
text(nsig+50,(1:length(nsig)*1.18)-0.3,labels=nsig)
mtext("No. significant GO terms")

x <- dbdb_res3$Proximal_tubule_cells
up <- head(subset(x,s.dist>0),20)
ups <- up$s.dist
names(ups) <- up$set
dn <- head(subset(x,s.dist<0),20)
dns <- dn$s.dist
names(dns) <- dn$set
dns <- rev(dns)
pws <- c(ups,dns)
pws <- sort(pws)
par(mar=c(5.1,25.1,4.1,2.1))
if ( length(pws)>2 ) {
  barplot(pws,horiz=TRUE,las=1,cex.names=0.6,xlab="ES")
  mtext("Proximal tubule GO terms")
}

mypws <- names(which(abs(pws)>0.4))

x <- dbdb_res3$Endothelial_cells
up <- head(subset(x,s.dist>0),20)
ups <- up$s.dist
names(ups) <- up$set
dn <- head(subset(x,s.dist<0),20)
dns <- dn$s.dist
names(dns) <- dn$set
dns <- rev(dns)
pws <- c(ups,dns)
pws <- sort(pws)
par(mar=c(5.1,25.1,4.1,2.1))
if ( length(pws)>2 ) {
  barplot(pws,horiz=TRUE,las=1,cex.names=0.7,xlab="ES")
  mtext("Endothelial cell GO terms")
}
mypws <- c(mypws,names(which(abs(pws)>0.4)))

x <- dbdb_res3$Podocytes
up <- head(subset(x,s.dist>0),20)
ups <- up$s.dist
names(ups) <- up$set
dn <- head(subset(x,s.dist<0),20)
dns <- dn$s.dist
names(dns) <- dn$set
dns <- rev(dns)
pws <- c(ups,dns)
pws <- sort(pws)
par(mar=c(5.1,25.1,4.1,2.1))
if ( length(pws)>2 ) {
  barplot(pws,horiz=TRUE,las=1,cex.names=0.7,xlab="ES")
  mtext("Podocyte GO terms")
}
mypws <- c(mypws,names(which(abs(pws)>0.4)))

x <- dbdb_res3$Distal_tubule_cells
up <- head(subset(x,s.dist>0),20)
ups <- up$s.dist
names(ups) <- up$set
dn <- head(subset(x,s.dist<0),20)
dns <- dn$s.dist
names(dns) <- dn$set
dns <- rev(dns)
pws <- c(ups,dns)
pws <- sort(pws)
par(mar=c(5.1,25.1,4.1,2.1))
if ( length(pws)>2 ) {
  barplot(pws,horiz=TRUE,las=1,cex.names=0.7,xlab="ES")
  mtext("Distal tubule cell GO terms")
}
mypws <- c(mypws,names(which(abs(pws)>0.4)))

Now focus only on proximal tubule cells in DBDB and LIRA50 contrasts.

pd <- read.table("Proximal_tubule_cells.DBDB_DVvCV.tsv")
pd$g <- sapply(strsplit(rownames(pd),"_"),"[[",1)

pl <- read.table("Proximal_tubule_cells.DBDB_DL50vDV.tsv")
pl$g <- sapply(strsplit(rownames(pl),"_"),"[[",1)

x <- list("diab"=pd,"lira"=pl)

#mitch_import(x=x,DEtype="seurat",geneIDcol="g",geneTable=gt)

#xm <- mitch_import(x=x,DEtype="seurat",geneIDcol="g",geneTable=gt) #outer join not working for report

xm <- mitch_import(x=x,DEtype="seurat",geneIDcol="g",geneTable=gt,joinType="inner")
## Note: Mean no. genes in input = 2397.5
## Note: no. genes in output = 1229
## Note: estimated proportion of input genes in output = 0.513
dim(xm)
## [1] 1229    2
mres <- mitch_calc(x=xm,genesets=go,minsetsize=5,cores=8,priority="effect")
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
head(mres$enrichment_result,10) %>%
  kbl(caption="Top GO terms, not filtered for FDR") %>%
  kable_paper("hover", full_width = F)
Top GO terms, not filtered for FDR
set setSize pMANOVA s.diab s.lira p.diab p.lira s.dist SD p.adjustMANOVA
1252 GOBP POTASSIUM ION HOMEOSTASIS 6 0.0003445 0.8702644 -0.7552467 0.0002245 0.0013751 1.1522836 1.1494099 0.0079282
1253 GOBP POTASSIUM ION IMPORT ACROSS PLASMA MEMBRANE 5 0.0029754 0.7228758 -0.8039216 0.0051804 0.0018696 1.0811287 1.0796088 0.0392032
197 GOBP CELLULAR RESPONSE TO OXYGEN RADICAL 5 0.0056948 -0.7315359 0.7138889 0.0046665 0.0057672 1.0221459 1.0220697 0.0582641
1782 GOBP RESPONSE TO OXYGEN RADICAL 5 0.0056948 -0.7315359 0.7138889 0.0046665 0.0057672 1.0221459 1.0220697 0.0582641
218 GOBP CELL CELL ADHESION VIA PLASMA MEMBRANE ADHESION MOLECULES 5 0.0051827 0.5264706 -0.8303922 0.0419003 0.0013134 0.9832204 0.9594468 0.0566743
2423 GOMF FERROUS IRON BINDING 6 0.0048767 -0.6124285 0.7151812 0.0095176 0.0024502 0.9415693 0.9387618 0.0539831
1080 GOBP POSITIVE REGULATION OF CARBOHYDRATE METABOLIC PROCESS 9 0.0002610 0.7734062 -0.5016393 0.0000599 0.0093778 0.9218455 0.9015933 0.0061844
1715 GOBP RENAL SODIUM ION ABSORPTION 6 0.0016286 0.8412374 -0.3590897 0.0003625 0.1287786 0.9146725 0.8487594 0.0258497
1716 GOBP RENAL SODIUM ION TRANSPORT 6 0.0016286 0.8412374 -0.3590897 0.0003625 0.1287786 0.9146725 0.8487594 0.0258497
1874 GOBP SOMITE DEVELOPMENT 5 0.0122065 0.4918301 -0.7581699 0.0573755 0.0033618 0.9037248 0.8838835 0.0923523
if ( ! file.exists("prox_diab_lira50.html") ) {
  mitch_report(res=mres,outfile="prox_diab_lira50.html")
}

mres2 <- subset(mres$enrichment_result,p.adjustMANOVA<0.05)

head(mres2,30) %>%
  kbl(caption="Top GO terms after 5% FDR filtering") %>%
  kable_paper("hover", full_width = F)
Top GO terms after 5% FDR filtering
set setSize pMANOVA s.diab s.lira p.diab p.lira s.dist SD p.adjustMANOVA
1252 GOBP POTASSIUM ION HOMEOSTASIS 6 0.0003445 0.8702644 -0.7552467 0.0002245 0.0013751 1.1522836 1.1494099 0.0079282
1253 GOBP POTASSIUM ION IMPORT ACROSS PLASMA MEMBRANE 5 0.0029754 0.7228758 -0.8039216 0.0051804 0.0018696 1.0811287 1.0796088 0.0392032
1080 GOBP POSITIVE REGULATION OF CARBOHYDRATE METABOLIC PROCESS 9 0.0002610 0.7734062 -0.5016393 0.0000599 0.0093778 0.9218455 0.9015933 0.0061844
1715 GOBP RENAL SODIUM ION ABSORPTION 6 0.0016286 0.8412374 -0.3590897 0.0003625 0.1287786 0.9146725 0.8487594 0.0258497
1716 GOBP RENAL SODIUM ION TRANSPORT 6 0.0016286 0.8412374 -0.3590897 0.0003625 0.1287786 0.9146725 0.8487594 0.0258497
1285 GOBP PROTEIN FOLDING IN ENDOPLASMIC RETICULUM 5 0.0023459 -0.8581699 0.2098039 0.0008970 0.4178266 0.8834440 0.7551716 0.0333487
620 GOBP LOW DENSITY LIPOPROTEIN PARTICLE CLEARANCE 5 0.0036151 0.8431373 -0.2604575 0.0011041 0.3144738 0.8824503 0.7803593 0.0444205
1950 GOBP TRICARBOXYLIC ACID CYCLE 9 0.0015238 -0.5362477 0.6542805 0.0054737 0.0006944 0.8459578 0.8418306 0.0249644
1845 GOBP SERINE FAMILY AMINO ACID METABOLIC PROCESS 9 0.0012079 -0.4826958 0.6910747 0.0124332 0.0003390 0.8429588 0.8299811 0.0209845
2206 GOCC NUCLEOID 9 0.0018283 -0.5063752 0.6557377 0.0087271 0.0006754 0.8284973 0.8217379 0.0277653
479 GOBP GLYCEROLIPID CATABOLIC PROCESS 6 0.0037771 0.2646498 -0.7732352 0.2630470 0.0010515 0.8172712 0.7338955 0.0449426
1493 GOBP REGULATION OF GLYCOLYTIC PROCESS 8 0.0041791 0.6537674 -0.4868960 0.0013947 0.0174143 0.8151561 0.8065708 0.0486547
1951 GOBP TRICARBOXYLIC ACID METABOLIC PROCESS 6 0.0036781 -0.2300354 0.7669665 0.3306668 0.0011553 0.8007208 0.7049868 0.0448812
2520 GOMF NAD BINDING 13 0.0001660 -0.4277581 0.6621964 0.0078534 0.0000373 0.7883407 0.7707142 0.0044308
618 GOBP LONG TERM MEMORY 5 0.0039901 0.7794118 -0.0898693 0.0025696 0.7285843 0.7845758 0.6146745 0.0468571
2168 GOCC MICROBODY LUMEN 17 0.0000629 -0.5488255 0.5264997 0.0000961 0.0001835 0.7605336 0.7603697 0.0020720
2099 GOCC CYTOSOLIC SMALL RIBOSOMAL SUBUNIT 11 0.0000000 0.4036423 0.6429318 0.0209708 0.0002303 0.7591366 0.1692032 0.0000006
445 GOBP FAT SOLUBLE VITAMIN METABOLIC PROCESS 5 0.0024753 0.7535948 0.0281046 0.0035592 0.9135857 0.7541187 0.5129990 0.0350040
624 GOBP LYMPHOCYTE APOPTOTIC PROCESS 5 0.0022018 0.7509804 0.0464052 0.0036766 0.8577965 0.7524128 0.4982099 0.0316339
26 GOBP ADULT BEHAVIOR 11 0.0015231 0.4321541 -0.6139722 0.0134371 0.0004373 0.7508122 0.7397230 0.0249644
1249 GOBP POST GOLGI VESICLE MEDIATED TRANSPORT 11 0.0015250 0.4067771 -0.6203911 0.0199930 0.0003803 0.7418576 0.7263176 0.0249644
2551 GOMF OXIDOREDUCTASE ACTIVITY ACTING ON PAIRED DONORS WITH INCORPORATION OR REDUCTION OF MOLECULAR OXYGEN REDUCED FLAVIN OR FLAVOPROTEIN AS ONE DONOR AND INCORPORATION OF ONE ATOM OF OXYGEN 6 0.0001014 0.2946307 0.6786590 0.2127330 0.0040507 0.7398549 0.2715490 0.0030552
589 GOBP LEUKOCYTE APOPTOTIC PROCESS 6 0.0027178 0.7339875 -0.0801308 0.0018747 0.7347611 0.7383485 0.5756686 0.0367043
1714 GOBP RENAL ABSORPTION 10 0.0025876 0.6275636 -0.3862182 0.0006080 0.0351570 0.7368857 0.7168520 0.0360260
2460 GOMF INTRAMOLECULAR OXIDOREDUCTASE ACTIVITY 12 0.0016321 -0.5625856 0.4657628 0.0007701 0.0053948 0.7303681 0.7271521 0.0258497
1871 GOBP SODIUM ION TRANSMEMBRANE TRANSPORT 17 0.0000994 0.5788196 -0.4446709 0.0000387 0.0015954 0.7299071 0.7237171 0.0030552
2257 GOCC RESPIRATORY CHAIN COMPLEX IV 9 0.0000009 0.3693989 0.6218579 0.0558898 0.0012687 0.7233000 0.1785155 0.0000441
1905 GOBP SULFUR COMPOUND CATABOLIC PROCESS 12 0.0013084 -0.6020268 0.3900301 0.0003178 0.0198624 0.7173283 0.7014902 0.0225090
2 GOBP ACTIN FILAMENT BASED MOVEMENT 10 0.0037057 0.3474979 -0.6114848 0.0580888 0.0008384 0.7033267 0.6781032 0.0448834
687 GOBP MITOCHONDRIAL ELECTRON TRANSPORT CYTOCHROME C TO OXYGEN 10 0.0000015 0.2474159 0.6556194 0.1774147 0.0003409 0.7007506 0.2886434 0.0000736

Visualisations

colfunc <- colorRampPalette(c("blue", "white", "red"))

mres3 <- mres2[1:30,c("s.diab","s.lira")]
rownames(mres3) <- mres2$set[1:30]

heatmap.2(  as.matrix(mres3), col=colfunc(25),scale="none",
  trace="none",margins = c(6,25), cexCol=0.9,cexRow=0.6, main="Top GO terms")

itx1h <- mm2hs[match(itx1,mm2hs$Gene.name),3]
itx2h <- mm2hs[match(itx2,mm2hs$Gene.name),3]

itx1h <- itx1h[which(itx1h!="")]
itx2h <- itx2h[which(itx2h!="")]

mitx <- lapply(mres$detailed_sets, function(x) {
  gs <- rownames(x)
  itx1res <- length(which(gs %in% itx1h))
  itx2res <- length(which(gs %in% itx2h))
  out=c(itx1res,itx2res)
  out
} )

mitx[which(lapply(mitx,sum)>1)]
## $`GOBP POTASSIUM ION HOMEOSTASIS`
## [1] 2 0
## 
## $`GOMF MONOATOMIC ANION MONOATOMIC CATION SYMPORTER ACTIVITY`
## [1] 2 0
## 
## $`GOMF MONOATOMIC ANION SODIUM SYMPORTER ACTIVITY`
## [1] 2 0
## 
## $`GOBP RESPONSE TO VITAMIN D`
## [1] 2 0
mres$detailed_sets$`GOBP POTASSIUM ION HOMEOSTASIS`
##         diab lira
## ATP1A1   715 -362
## KLHL3    608 -339
## SLC12A1  609 -688
## SLC12A3  702 -626
## UMOD     646 -694
## WNK1     567 -560
mres$detailed_sets$`GOBP RESPONSE TO VITAMIN D`
##         diab lira
## CYP24A1  696 -553
## RXRA     162 -101
## SFRP1    656 -683
## SNW1     359  148
## SPP1     723 -696

Venn diagram for pathways

diab_up <- subset(dbdb_res_go$Proximal_tubule_cells,s.dist>0)$set
diab_dn <- subset(dbdb_res_go$Proximal_tubule_cells,s.dist<0)$set

lira_up <- subset(dbdb_res3$Proximal_tubule_cells,s.dist > 0)$set
lira_dn <- subset(dbdb_res3$Proximal_tubule_cells,s.dist < 0)$set

v1 <- list("Diab up"=diab_up, "Diab dn"=diab_dn,
  "Lira up"=lira_up,"Lira dn"=lira_dn)

plot(euler(v1),quantities = TRUE,main="GO terms - Proximal tubule cells")

message("GO terms upregulated in diabetes and lowered with lira50 at 5% FDR")
## GO terms upregulated in diabetes and lowered with lira50 at 5% FDR
intersect(diab_up,lira_dn)
##  [1] "GOBP SODIUM ION TRANSMEMBRANE TRANSPORT"                       
##  [2] "GOBP ACTIN FILAMENT BASED MOVEMENT"                            
##  [3] "GOBP DEVELOPMENTAL CELL GROWTH"                                
##  [4] "GOBP SODIUM ION TRANSPORT"                                     
##  [5] "GOBP DEVELOPMENTAL GROWTH INVOLVED IN MORPHOGENESIS"           
##  [6] "GOMF GTPASE ACTIVATOR ACTIVITY"                                
##  [7] "GOMF GUANYL NUCLEOTIDE EXCHANGE FACTOR ACTIVITY"               
##  [8] "GOMF ACTIN FILAMENT BINDING"                                   
##  [9] "GOBP REGULATION OF DEVELOPMENTAL GROWTH"                       
## [10] "GOMF NUCLEOSIDE TRIPHOSPHATASE REGULATOR ACTIVITY"             
## [11] "GOMF ACTIN BINDING"                                            
## [12] "GOBP METAL ION TRANSPORT"                                      
## [13] "GOBP DEVELOPMENTAL GROWTH"                                     
## [14] "GOBP CELL PROJECTION MORPHOGENESIS"                            
## [15] "GOCC BASAL PART OF CELL"                                       
## [16] "GOBP CELL MORPHOGENESIS INVOLVED IN NEURON DIFFERENTIATION"    
## [17] "GOBP REGULATION OF TRANSMEMBRANE TRANSPORT"                    
## [18] "GOBP REGULATION OF NEURON PROJECTION DEVELOPMENT"              
## [19] "GOBP REGULATION OF NERVOUS SYSTEM DEVELOPMENT"                 
## [20] "GOBP ENZYME LINKED RECEPTOR PROTEIN SIGNALING PATHWAY"         
## [21] "GOBP GENERATION OF NEURONS"                                    
## [22] "GOBP POSITIVE REGULATION OF RNA METABOLIC PROCESS"             
## [23] "GOBP CELL MORPHOGENESIS"                                       
## [24] "GOBP POSITIVE REGULATION OF TRANSCRIPTION BY RNA POLYMERASE II"
## [25] "GOBP ACTIN FILAMENT BASED PROCESS"                             
## [26] "GOBP CELL JUNCTION ORGANIZATION"                               
## [27] "GOBP REGULATION OF ANATOMICAL STRUCTURE MORPHOGENESIS"         
## [28] "GOBP CELL PROJECTION ORGANIZATION"                             
## [29] "GOBP NEURON DEVELOPMENT"                                       
## [30] "GOBP NEUROGENESIS"                                             
## [31] "GOMF CHROMATIN BINDING"                                        
## [32] "GOBP REGULATION OF CELL ADHESION"                              
## [33] "GOBP REGULATION OF CELL PROJECTION ORGANIZATION"               
## [34] "GOMF CYTOSKELETAL PROTEIN BINDING"                             
## [35] "GOBP POSITIVE REGULATION OF DEVELOPMENTAL PROCESS"             
## [36] "GOBP REGULATION OF MULTICELLULAR ORGANISMAL DEVELOPMENT"       
## [37] "GOCC CHROMATIN"                                                
## [38] "GOMF SEQUENCE SPECIFIC DNA BINDING"                            
## [39] "GOBP REGULATION OF CELL DIFFERENTIATION"                       
## [40] "GOBP CELL ADHESION"                                            
## [41] "GOBP GROWTH"                                                   
## [42] "GOMF ENZYME REGULATOR ACTIVITY"                                
## [43] "GOCC PLASMA MEMBRANE REGION"                                   
## [44] "GOBP CYTOSKELETON ORGANIZATION"
message("GO terms downregulated in diabetes and raise with lira50 at 5% FDR")
## GO terms downregulated in diabetes and raise with lira50 at 5% FDR
intersect(diab_dn,lira_up)
##  [1] "GOCC MICROBODY LUMEN"                                        
##  [2] "GOMF OXIDOREDUCTASE ACTIVITY ACTING ON CH OH GROUP OF DONORS"
##  [3] "GOBP GLUTATHIONE METABOLIC PROCESS"                          
##  [4] "GOBP SULFUR COMPOUND METABOLIC PROCESS"                      
##  [5] "GOCC FICOLIN 1 RICH GRANULE LUMEN"                           
##  [6] "GOBP AMINO ACID CATABOLIC PROCESS"                           
##  [7] "GOMF LYASE ACTIVITY"                                         
##  [8] "GOBP ALPHA AMINO ACID METABOLIC PROCESS"                     
##  [9] "GOBP AMINO ACID METABOLIC PROCESS"                           
## [10] "GOCC MITOCHONDRIAL MATRIX"                                   
## [11] "GOBP ORGANIC ACID CATABOLIC PROCESS"                         
## [12] "GOBP ORGANIC ACID METABOLIC PROCESS"                         
## [13] "GOMF OXIDOREDUCTASE ACTIVITY"

Now with FDR<0.01 for more strict filtering.

diab_up <- subset(dbdb_res_go$Proximal_tubule_cells,s.dist>0&p.adjustANOVA<0.01)$set
diab_dn <- subset(dbdb_res_go$Proximal_tubule_cells,s.dist<0&p.adjustANOVA<0.01)$set

lira_up <- subset(dbdb_res3$Proximal_tubule_cells,s.dist>0&p.adjustANOVA<0.01)$set
lira_dn <- subset(dbdb_res3$Proximal_tubule_cells,s.dist<0&p.adjustANOVA<0.01)$set

v1 <- list("Diab up"=diab_up, "Diab dn"=diab_dn,
  "Lira up"=lira_up,"Lira dn"=lira_dn)

plot(euler(v1),quantities = TRUE,main="GO terms - Proximal tubule cells")

message("GO terms upregulated in diabetes and lowered with lira50 at 1% FDR")
## GO terms upregulated in diabetes and lowered with lira50 at 1% FDR
intersect(diab_up,lira_dn)
##  [1] "GOBP DEVELOPMENTAL CELL GROWTH"                            
##  [2] "GOBP SODIUM ION TRANSPORT"                                 
##  [3] "GOMF NUCLEOSIDE TRIPHOSPHATASE REGULATOR ACTIVITY"         
##  [4] "GOMF ACTIN BINDING"                                        
##  [5] "GOBP METAL ION TRANSPORT"                                  
##  [6] "GOBP DEVELOPMENTAL GROWTH"                                 
##  [7] "GOBP CELL PROJECTION MORPHOGENESIS"                        
##  [8] "GOBP CELL MORPHOGENESIS INVOLVED IN NEURON DIFFERENTIATION"
##  [9] "GOBP ENZYME LINKED RECEPTOR PROTEIN SIGNALING PATHWAY"     
## [10] "GOBP GENERATION OF NEURONS"                                
## [11] "GOBP POSITIVE REGULATION OF RNA METABOLIC PROCESS"         
## [12] "GOBP CELL MORPHOGENESIS"                                   
## [13] "GOBP ACTIN FILAMENT BASED PROCESS"                         
## [14] "GOBP REGULATION OF ANATOMICAL STRUCTURE MORPHOGENESIS"     
## [15] "GOBP NEUROGENESIS"                                         
## [16] "GOBP POSITIVE REGULATION OF DEVELOPMENTAL PROCESS"         
## [17] "GOBP REGULATION OF MULTICELLULAR ORGANISMAL DEVELOPMENT"
message("GO terms downregulated in diabetes and raise with lira50 at 1% FDR")
## GO terms downregulated in diabetes and raise with lira50 at 1% FDR
intersect(diab_dn,lira_up)
## [1] "GOBP SULFUR COMPOUND METABOLIC PROCESS"
## [2] "GOBP AMINO ACID METABOLIC PROCESS"     
## [3] "GOCC MITOCHONDRIAL MATRIX"             
## [4] "GOBP ORGANIC ACID METABOLIC PROCESS"

Session information

sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## 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    beeswarm_0.4.0     
## [4] limma_3.65.1        eulerr_7.0.2        mitch_1.21.3       
## [7] gplots_3.2.0        RhpcBLASctl_0.23-42 plyr_1.8.9         
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6          xfun_0.52             bslib_0.9.0          
##  [4] ggplot2_3.5.2         htmlwidgets_1.6.4     caTools_1.18.3       
##  [7] GGally_2.2.1          lattice_0.22-7        vctrs_0.6.5          
## [10] tools_4.5.1           bitops_1.0-9          generics_0.1.4       
## [13] parallel_4.5.1        tibble_3.3.0          pkgconfig_2.0.3      
## [16] KernSmooth_2.23-26    polylabelr_0.3.0      RColorBrewer_1.1-3   
## [19] lifecycle_1.0.4       compiler_4.5.1        farver_2.1.2         
## [22] stringr_1.5.1         textshaping_1.0.1     statmod_1.5.0        
## [25] httpuv_1.6.16         htmltools_0.5.8.1     sass_0.4.10          
## [28] yaml_2.3.10           later_1.4.2           pillar_1.10.2        
## [31] jquerylib_0.1.4       tidyr_1.3.1           MASS_7.3-65          
## [34] cachem_1.1.0          mime_0.13             ggstats_0.9.0        
## [37] network_1.19.0        tidyselect_1.2.1      digest_0.6.37        
## [40] stringi_1.8.7         dplyr_1.1.4           reshape2_1.4.4       
## [43] purrr_1.0.4           labeling_0.4.3        polyclip_1.10-7      
## [46] fastmap_1.2.0         grid_4.5.1            cli_3.6.5            
## [49] magrittr_2.0.3        dichromat_2.0-0.1     withr_3.0.2          
## [52] scales_1.4.0          promises_1.3.3        rmarkdown_2.29       
## [55] gridExtra_2.3         coda_0.19-4.1         shiny_1.10.0         
## [58] evaluate_1.0.3        knitr_1.50            viridisLite_0.4.2    
## [61] rlang_1.1.6           Rcpp_1.0.14           xtable_1.8-4         
## [64] glue_1.8.0            echarts4r_0.4.5       xml2_1.3.8           
## [67] svglite_2.2.1         rstudioapi_0.17.1     jsonlite_2.0.0       
## [70] R6_2.6.1              statnet.common_4.12.0 systemfonts_1.2.3