Time course analysis of PADDI genomics data.
suppressPackageStartupMessages({
library("gplots")
library("reshape2")
library("WGCNA")
library("dplyr")
library("DESeq2")
library("mitch")
library("MASS")
library("eulerr")
library("kableExtra")
library("gplots")
library("RhpcBLASctl")
})
RhpcBLASctl::blas_set_num_threads(1)
load("tca_pairwise.Rdata")
go <- mitch::gmt_import("c5.go.v2024.1.Hs.symbols.gmt")
names(go) <- gsub("_"," ",names(go))
gt <- read.table("../ref/gencode.v38.genetable.tsv")
rownames(gt) <- paste(gt[,1],gt[,2])
gt[,1] <- rownames(gt)
We will run timecourse analysis with a simple pairwise approach. For each of the groups below, this will involve three comparisons:
T0 vs EOS
EOS vs POD1
T0 vs POD1
The groups/subgroups we will look at are:
Timecourse in low CRP group
Timecourse in high CRP group
Timecourse in low CRP group and treatment group A
Timecourse in low CRP group and treatment group B
Timecourse in high CRP group and treatment group A
Timecourse in high CRP group and treatment group B
In this report I will be loading in the DESeq2 objects and conducting downstream analysis.
mitch1d <- function(de,dename,gs,gt=gt) {
m <- mitch_import(x=de, DEtype="deseq2", geneTable=gt )
mres <- mitch_calc(x=m,genesets=gs,minsetsize=5,cores=8,priority="effect")
mtop_up <- head(subset (mres$enrichment_result,p.adjustANOVA<0.05 & s.dist>0),15)
mtop_up |> kbl(caption=paste(dename,"GO up")) |> kable_paper("hover", full_width = F)
mtop_dn <- head(subset (mres$enrichment_result,p.adjustANOVA<0.05 & s.dist<0),15)
mtop_dn |> kbl(caption=paste(dename,"GO dn")) |> kable_paper("hover", full_width = F)
mtop <- rbind(mtop_up, mtop_dn[nrow(mtop_dn):1,])
vec <- mtop$s.dist
names(vec) <- mtop$set
vec <- sort(vec)
par( mar = c(5.1, 25.1, 4.1, 2.1) )
barplot(vec,horiz=TRUE,las=1,cex.names=0.7,xlab="ES",main=dename)
par( mar = c(5.1, 4.1, 4.1, 2.1) )
# mitch_report(res=mres,outfile=paste(dename,"_mitchreport.html",sep=""),overwrite=TRUE)
}
#mitch1d(de=tc_lo_t0veos_adj, dename="tc_lo_t0veos_adj", gs=go, gt=gt)
mitch2d <- function(de1, de2, de1name, de2name, gs, gt=gt, ngenes=40) {
l <-list(de1,de2)
names(l) <- c(de1name, de2name)
m <- mitch_import(x=l, DEtype="deseq2", geneTable=gt )
mres <- mitch_calc(x=m,genesets=gs, minsetsize=5, cores=8, priority="effect")
mtop <- head( subset(mres$enrichment_result, p.adjustMANOVA < 0.05), ngenes)
mtop |> kbl(caption=paste("GO")) |> kable_paper("hover", full_width = F)
mx <- as.matrix(mtop[,4:5])
rownames(mx) <- mtop$set
mx <- rbind(mx,t(as.matrix(c(-1,1))))
rownames(mx)[nrow(mx)] <- "minmax"
colfunc <- colorRampPalette(c("darkblue","blue","lightblue", "white", "pink","red","darkred"))
heatmap.2(mx,col=colfunc(50),scale="none",trace="none",mar=c(8,20),cexRow=0.6, cexCol=0.8)
}
#mitch2d(de=tc_lo_t0veos_adj, de2=tc_lo_eosvpod1_adj,
# de1name="tc_lo_t0veos_adj", de2name="tc_lo_eosvpod1_adj", gs=go, gt=gt)
tc_lo_t0veos_adj tc_lo_eosvpod1_adj tc_lo_t0vpod1_adj
mitch1d(de=tc_lo_t0veos_adj, dename="tc_lo_t0veos_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 22144
## Note: no. genes in output = 22075
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_lo_eosvpod1_adj, dename="tc_lo_eosvpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21793
## Note: no. genes in output = 21726
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_lo_t0vpod1_adj, dename="tc_lo_t0vpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21672
## Note: no. genes in output = 21608
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
tc_hi_t0veos_adj tc_hi_eosvpod1_adj tc_hi_t0vpod1_adj
mitch1d(de=tc_hi_t0veos_adj, dename="tc_hi_t0veos_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21772
## Note: no. genes in output = 21707
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_hi_eosvpod1_adj, dename="tc_hi_eosvpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21487
## Note: no. genes in output = 21427
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_hi_t0vpod1_adj, dename="tc_hi_t0vpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21614
## Note: no. genes in output = 21550
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
tc_lo_a_t0veos_adj tc_lo_a_eosvpod1_adj tc_lo_a_t0vpod1_adj
mitch1d(de=tc_lo_a_t0veos_adj, dename="tc_lo_a_t0veos_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 22167
## Note: no. genes in output = 22099
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_lo_a_eosvpod1_adj, dename="tc_lo_a_eosvpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21814
## Note: no. genes in output = 21746
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_lo_a_t0vpod1_adj, dename="tc_lo_a_t0vpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21696
## Note: no. genes in output = 21632
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
tc_lo_b_t0veos_adj tc_lo_b_eosvpod1_adj tc_lo_b_t0vpod1_adj
mitch1d(de=tc_lo_b_t0veos_adj, dename="tc_lo_b_t0veos_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 22091
## Note: no. genes in output = 22024
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_lo_b_eosvpod1_adj, dename="tc_lo_b_eosvpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21736
## Note: no. genes in output = 21672
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_lo_b_t0vpod1_adj, dename="tc_lo_b_t0vpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21715
## Note: no. genes in output = 21651
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
tc_hi_a_t0veos_adj tc_hi_a_eosvpod1_adj tc_hi_a_t0vpod1_adj
mitch1d(de=tc_hi_a_t0veos_adj, dename="tc_hi_a_t0veos_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 22227
## Note: no. genes in output = 22154
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_hi_a_eosvpod1_adj, dename="tc_hi_a_eosvpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21654
## Note: no. genes in output = 21586
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_hi_a_t0vpod1_adj, dename="tc_hi_a_t0vpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21860
## Note: no. genes in output = 21794
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
tc_hi_b_t0veos_adj tc_hi_b_eosvpod1_adj tc_hi_b_t0vpod1_adj
mitch1d(de=tc_hi_b_t0veos_adj, dename="tc_hi_b_t0veos_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21828
## Note: no. genes in output = 21765
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_hi_b_eosvpod1_adj, dename="tc_hi_b_eosvpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21618
## Note: no. genes in output = 21556
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch1d(de=tc_hi_b_t0vpod1_adj, dename="tc_hi_b_t0vpod1_adj", gs=go, gt=gt)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 21635
## Note: no. genes in output = 21573
## Note: estimated proportion of input genes in output = 0.997
## Note: Enrichments with large effect sizes may not be
## statistically significant.
Few changes in T0 vs EOS and more in EOS vs POD1.
TODO: * Euler diagram * scatter plot of stat values * pathway analysis * bar charts and meaning of the DGE results
mitch2d(de1=tc_lo_t0veos_adj, de2=tc_lo_eosvpod1_adj,
de1name="tc_lo_t0veos_adj", de2name="tc_lo_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21968.5
## Note: no. genes in output = 21661
## Note: estimated proportion of input genes in output = 0.986
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_hi_t0veos_adj, de2=tc_hi_eosvpod1_adj,
de1name="tc_hi_t0veos_adj", de2name="tc_hi_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21629.5
## Note: no. genes in output = 21395
## Note: estimated proportion of input genes in output = 0.989
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_lo_a_t0veos_adj, de2=tc_lo_a_eosvpod1_adj,
de1name="tc_lo_a_t0veos_adj", de2name="tc_lo_a_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21990.5
## Note: no. genes in output = 21652
## Note: estimated proportion of input genes in output = 0.985
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de=tc_lo_b_t0veos_adj, de2=tc_lo_b_eosvpod1_adj,
de1name="tc_lo_b_t0veos_adj", de2name="tc_lo_b_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21913.5
## Note: no. genes in output = 21594
## Note: estimated proportion of input genes in output = 0.985
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_hi_a_t0veos_adj, de2=tc_hi_a_eosvpod1_adj,
de1name="tc_hi_a_t0veos_adj", de2name="tc_hi_a_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21940.5
## Note: no. genes in output = 21470
## Note: estimated proportion of input genes in output = 0.979
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_hi_b_t0veos_adj, de2=tc_hi_b_eosvpod1_adj,
de1name="tc_hi_b_t0veos_adj", de2name="tc_hi_b_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21723
## Note: no. genes in output = 21372
## Note: estimated proportion of input genes in output = 0.984
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_lo_t0veos_adj, de2=tc_hi_t0veos_adj,
de1name="tc_lo_t0veos_adj", de2name="tc_hi_t0veos_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21958
## Note: no. genes in output = 21623
## Note: estimated proportion of input genes in output = 0.985
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_lo_eosvpod1_adj, de2=tc_hi_eosvpod1_adj ,
de1name="tc_lo_eosvpod1_adj", de2name="tc_hi_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21640
## Note: no. genes in output = 21313
## Note: estimated proportion of input genes in output = 0.985
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_lo_t0vpod1_adj, de2=tc_hi_t0vpod1_adj,
de1name="tc_lo_t0vpod1_adj", de2name="tc_hi_t0vpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21643
## Note: no. genes in output = 21355
## Note: estimated proportion of input genes in output = 0.987
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_lo_a_t0veos_adj, de2=tc_lo_b_t0veos_adj,
de1name="tc_lo_a_t0veos_adj", de2name="tc_lo_b_t0veos_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 22129
## Note: no. genes in output = 21730
## Note: estimated proportion of input genes in output = 0.982
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_lo_a_eosvpod1_adj, de2=tc_lo_b_eosvpod1_adj,
de1name="tc_lo_a_eosvpod1_adj", de2name="tc_lo_b_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21775
## Note: no. genes in output = 21335
## Note: estimated proportion of input genes in output = 0.98
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_lo_a_t0vpod1_adj, de2=tc_lo_b_t0vpod1_adj,
de1name="tc_lo_a_t0vpod1_adj", de2name="tc_lo_b_t0vpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21705.5
## Note: no. genes in output = 21356
## Note: estimated proportion of input genes in output = 0.984
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_hi_a_t0veos_adj, de2=tc_hi_b_t0veos_adj,
de1name="tc_hi_a_t0veos_adj", de2name="tc_hi_b_t0veos_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 22027.5
## Note: no. genes in output = 21477
## Note: estimated proportion of input genes in output = 0.975
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_hi_a_eosvpod1_adj, de2=tc_hi_b_eosvpod1_adj,
de1name="tc_hi_a_eosvpod1_adj", de2name="tc_hi_b_eosvpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21636
## Note: no. genes in output = 21148
## Note: estimated proportion of input genes in output = 0.977
## Note: Enrichments with large effect sizes may not be
## statistically significant.
mitch2d(de1=tc_hi_a_t0vpod1_adj, de2=tc_hi_b_t0vpod1_adj,
de1name="tc_hi_a_t0vpod1_adj", de2name="tc_hi_b_t0vpod1_adj", gs=go, gt=gt)
## Note: Mean no. genes in input = 21747.5
## Note: no. genes in output = 21294
## Note: estimated proportion of input genes in output = 0.979
## Note: Enrichments with large effect sizes may not be
## statistically significant.
This one is getting long, so I will continue downstream analysis in a separate script.
For reproducibility
sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=en_US.UTF-8
## [9] LC_ADDRESS=en_US.UTF-8 LC_TELEPHONE=en_US.UTF-8
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=en_US.UTF-8
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] RhpcBLASctl_0.23-42 gtools_3.9.5
## [3] xlsx_0.6.5 DT_0.34.0
## [5] ggplot2_4.0.3 kableExtra_1.4.0
## [7] beeswarm_0.4.0 eulerr_7.1.0
## [9] MASS_7.3-65 mitch_1.22.1
## [11] DESeq2_1.50.2 SummarizedExperiment_1.40.0
## [13] Biobase_2.70.0 MatrixGenerics_1.22.0
## [15] matrixStats_1.5.0 GenomicRanges_1.62.1
## [17] Seqinfo_1.0.0 IRanges_2.44.0
## [19] S4Vectors_0.48.1 BiocGenerics_0.56.0
## [21] generics_0.1.4 dplyr_1.2.1
## [23] WGCNA_1.74 fastcluster_1.3.0
## [25] dynamicTreeCut_1.63-1 reshape2_1.4.5
## [27] gplots_3.3.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.18.0 jsonlite_2.0.0
## [4] magrittr_2.0.5 farver_2.1.2 rmarkdown_2.31
## [7] vctrs_0.7.3 base64enc_0.1-6 htmltools_0.5.9
## [10] S4Arrays_1.10.1 progress_1.2.3 SparseArray_1.10.10
## [13] Formula_1.2-5 sass_0.4.10 KernSmooth_2.23-26
## [16] bslib_0.10.0 htmlwidgets_1.6.4 plyr_1.8.9
## [19] echarts4r_0.5.0 impute_1.84.0 cachem_1.1.0
## [22] mime_0.13 lifecycle_1.0.5 iterators_1.0.14
## [25] pkgconfig_2.0.3 Matrix_1.7-5 R6_2.6.1
## [28] fastmap_1.2.0 shiny_1.13.0 digest_0.6.39
## [31] colorspace_2.1-2 GGally_2.4.0 textshaping_1.0.5
## [34] crosstalk_1.2.2 Hmisc_5.2-5 labeling_0.4.3
## [37] polyclip_1.10-7 abind_1.4-8 compiler_4.5.2
## [40] withr_3.0.2 doParallel_1.0.17 htmlTable_2.5.0
## [43] S7_0.2.2 backports_1.5.1 BiocParallel_1.44.0
## [46] ggstats_0.13.0 DelayedArray_0.36.1 caTools_1.18.3
## [49] tools_4.5.2 foreign_0.8-91 otel_0.2.0
## [52] httpuv_1.6.17 nnet_7.3-20 glue_1.8.1
## [55] promises_1.5.0 grid_4.5.2 polylabelr_1.0.0
## [58] checkmate_2.3.4 cluster_2.1.8.2 gtable_0.3.6
## [61] preprocessCore_1.72.0 tidyr_1.3.2 hms_1.1.4
## [64] data.table_1.18.2.1 xml2_1.5.2 XVector_0.50.0
## [67] foreach_1.5.2 pillar_1.11.1 stringr_1.6.0
## [70] limma_3.66.0 later_1.4.8 rJava_1.0-18
## [73] splines_4.5.2 lattice_0.22-9 survival_3.8-6
## [76] tidyselect_1.2.1 locfit_1.5-9.12 knitr_1.51
## [79] gridExtra_2.3 svglite_2.2.2 xfun_0.57
## [82] statmod_1.5.1 stringi_1.8.7 statnet.common_4.13.0
## [85] yaml_2.3.12 evaluate_1.0.5 codetools_0.2-20
## [88] xlsxjars_0.9.0 tibble_3.3.1 BiocManager_1.30.27
## [91] cli_3.6.6 rpart_4.1.27 xtable_1.8-8
## [94] systemfonts_1.3.2 jquerylib_0.1.4 network_1.20.0
## [97] Rcpp_1.1.1-1.1 coda_0.19-4.1 parallel_4.5.2
## [100] prettyunits_1.2.0 bitops_1.0-9 viridisLite_0.4.3
## [103] scales_1.4.0 purrr_1.2.2 crayon_1.5.3
## [106] rlang_1.2.0