suppressPackageStartupMessages({
library("plyr")
library("RhpcBLASctl")
library("gplots")
library("mitch")
library("eulerr")
library("limma")
library("beeswarm")
library("kableExtra")
})
RhpcBLASctl::blas_set_num_threads(1)
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
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
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
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:
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)
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)
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
Now take a look at effect of LIRA25 in DBDB mice.
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
Proximal tubules:
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)
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)
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
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"
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