PADDI RNA expression analysis - Timecourse Analysis - Downstream
Mark Ziemann
2025-04-17
Introduction
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")
})Load the data
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)Overview
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.
Functions
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)Object names
Timecourse in low CRP group
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.
Timecourse in high CRP group
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.
Timecourse in low CRP group and treatment group A
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.
Timecourse in low CRP group and treatment group B
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.
Timecourse in high CRP group and treatment group A
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.
Timecourse in high CRP group and treatment group B
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.
low CRP group
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.
high CRP group
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.
Timecourse in low CRP group and treatment group A
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.
Timecourse in low CRP group and treatment group B
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.
Timecourse in high CRP group and treatment group A
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.
Timecourse in high CRP group and treatment group B
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.
Comparison of high and low CRP groups T0 vs EOS
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.
Comparison of high and low CRP groups EOS vs POD1
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.
Comparison of high and low CRP groups T0 vs POD1
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.
Comparison of treatment groups A and B in low CRP T0 vs EOS
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.
Comparison of treatment groups A and B in low CRP EOS vs POD1
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.
Comparison of treatment groups A and B in low CRP T0 vs POD1
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.
Comparison of treatment groups A and B in high CRP T0 vs EOS
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.
Comparison of treatment groups A and B in high CRP EOS vs POD1
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.
Comparison of treatment groups A and B in high CRP T0 vs POD1
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.
STOP HERE
This one is getting long, so I will continue downstream analysis in a separate script.
Session information
For reproducibility
sessionInfo()## R version 4.4.3 (2025-02-28)
## 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
##
## 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] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] kableExtra_1.4.0 eulerr_7.0.2
## [3] MASS_7.3-65 mitch_1.19.3
## [5] DESeq2_1.44.0 SummarizedExperiment_1.34.0
## [7] Biobase_2.64.0 MatrixGenerics_1.16.0
## [9] matrixStats_1.5.0 GenomicRanges_1.56.2
## [11] GenomeInfoDb_1.40.1 IRanges_2.38.1
## [13] S4Vectors_0.42.1 BiocGenerics_0.50.0
## [15] dplyr_1.1.4 WGCNA_1.73
## [17] fastcluster_1.2.6 dynamicTreeCut_1.63-1
## [19] reshape2_1.4.4 gplots_3.2.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 bitops_1.0-9 gridExtra_2.3
## [4] echarts4r_0.4.5 rlang_1.1.5 magrittr_2.0.3
## [7] compiler_4.4.3 RSQLite_2.3.9 systemfonts_1.2.1
## [10] png_0.1-8 vctrs_0.6.5 stringr_1.5.1
## [13] pkgconfig_2.0.3 crayon_1.5.3 fastmap_1.2.0
## [16] backports_1.5.0 XVector_0.44.0 caTools_1.18.3
## [19] promises_1.3.2 rmarkdown_2.29 UCSC.utils_1.0.0
## [22] preprocessCore_1.66.0 purrr_1.0.4 bit_4.6.0
## [25] xfun_0.51 zlibbioc_1.50.0 cachem_1.1.0
## [28] jsonlite_1.9.1 blob_1.2.4 later_1.4.1
## [31] DelayedArray_0.30.1 BiocParallel_1.38.0 parallel_4.4.3
## [34] cluster_2.1.8.1 R6_2.6.1 RColorBrewer_1.1-3
## [37] bslib_0.9.0 stringi_1.8.4 GGally_2.2.1
## [40] rpart_4.1.24 jquerylib_0.1.4 Rcpp_1.0.14
## [43] iterators_1.0.14 knitr_1.50 base64enc_0.1-3
## [46] httpuv_1.6.15 Matrix_1.7-3 splines_4.4.3
## [49] nnet_7.3-20 tidyselect_1.2.1 rstudioapi_0.17.1
## [52] abind_1.4-8 yaml_2.3.10 doParallel_1.0.17
## [55] codetools_0.2-20 lattice_0.22-6 tibble_3.2.1
## [58] plyr_1.8.9 shiny_1.10.0 KEGGREST_1.44.1
## [61] evaluate_1.0.3 foreign_0.8-89 survival_3.8-3
## [64] ggstats_0.9.0 xml2_1.3.8 Biostrings_2.72.1
## [67] pillar_1.10.1 KernSmooth_2.23-26 checkmate_2.3.2
## [70] foreach_1.5.2 generics_0.1.3 ggplot2_3.5.1
## [73] munsell_0.5.1 scales_1.3.0 xtable_1.8-4
## [76] gtools_3.9.5 glue_1.8.0 Hmisc_5.2-3
## [79] tools_4.4.3 data.table_1.17.0 locfit_1.5-9.12
## [82] grid_4.4.3 impute_1.80.0 tidyr_1.3.1
## [85] AnnotationDbi_1.66.0 colorspace_2.1-1 GenomeInfoDbData_1.2.12
## [88] beeswarm_0.4.0 htmlTable_2.4.3 Formula_1.2-5
## [91] cli_3.6.4 viridisLite_0.4.2 S4Arrays_1.4.1
## [94] svglite_2.1.3 gtable_0.3.6 sass_0.4.9
## [97] digest_0.6.37 SparseArray_1.4.8 htmlwidgets_1.6.4
## [100] memoise_2.0.1.9000 htmltools_0.5.8.1 lifecycle_1.0.4
## [103] httr_1.4.7 mime_0.13 GO.db_3.19.1
## [106] bit64_4.6.0-1