Introduction

The differential analysis was already conducted using the qc.Rmd script. Here we are doing downstream analysis including pathway analysis and charts.

suppressPackageStartupMessages({
  library("gplots")
  library("reshape2")
  library("WGCNA")
  library("dplyr")
  library("DESeq2")
  library("mitch")
  library("MASS")
  library("kableExtra")
  library("ggplot2")
  library("eulerr")
  library("DT")
  library("xlsx")
})

load("qc.Rds")

Functions

volcanoplot <- function(x,title) {
  tot <- nrow(x)
  sig <- subset(x,padj<0.05)
  nsig <- nrow(subset(x,padj<0.05))
  dn <- nrow(subset(x,padj<0.05 & log2FoldChange<0))
  up <- nrow(subset(x,padj<0.05 & log2FoldChange>0))
  header <- paste(title,":", tot, "genes,", nsig ,"@5% FDR,", up, "up,", dn, "dn")

  plot(x$log2FoldChange,-log10(x$pvalue),bty="n",
    cex=0.5, pch=19, xlab="log2 fold change",
    ylab="-log10 p-value", main=header)

  points(sig$log2FoldChange,-log10(sig$pvalue),
    cex=0.5, pch=19, col="red")
}

smearplot <- function(x,title) {
  tot <- nrow(x)
  sig <- subset(x,padj<0.05)
  nsig <- nrow(subset(x,padj<0.05))
  dn <- nrow(subset(x,padj<0.05 & log2FoldChange<0))
  up <- nrow(subset(x,padj<0.05 & log2FoldChange>0))
  header <- paste(title,":", tot, "genes,", nsig ,"@5% FDR,", up, "up,", dn, "dn")

  plot(log10(x$baseMean),x$log2FoldChange,bty="n",
    cex=0.5, pch=19, xlab="log10 base mean",
    ylab="log2 fold change", main=header)

  points(log10(sig$baseMean),sig$log2FoldChange,
    cex=0.5, pch=19, col="red")
}

debar <- function(l) {
  nn <- sapply(l,function(x) {
    nup <- nrow(subset(x,padj<0.05 & log2FoldChange>0 ))
    ndn <- nrow(subset(x,padj<0.05 & log2FoldChange<0 )) *-1
    return(c(nup,ndn))
  })
  rownames(nn) <- c("up","dn")
  par(mar = c(8.1, 4.1, 4.1, 2.1))
  barplot(nn[1,],col="blue",beside = TRUE, border=TRUE,ylim=c(min(nn),max(nn)) , las=2)
  barplot(nn[2,],add=TRUE,col="red",names.arg = "", ylab="", yaxt="n")
  mtext("blue=up,red=down")
  par(mar = c(5.1, 4.1, 4.1, 2.1))
}

myeuler <- function(l) {
  ups <- lapply(l,function(x) { rownames(subset(x,padj<0.05 & log2FoldChange>0) )})
  dns <- lapply(l,function(x) { rownames(subset(x,padj<0.05 & log2FoldChange<0) )})
  names(ups) <- paste(names(ups),"up")
  names(dns) <- paste(names(dns),"dn")
  ll <- c(ups,dns)
  plot(euler(ll),quantities = TRUE)
}

writeTSV <- function(l) {
  lapply(1:length(l), function(i) {
    filename <- paste(names(l)[i],".tsv",sep="")
    write.table(x=l[[i]],file=filename,sep="\t",quote=FALSE)
  })
}

DE results

Unadjusted for clinical covariates

Unstratified, unadjusted

  • crp_t0
  • crp_eos
  • crp_pod1
  • dex_t0
  • dex_eos
  • dex_pod1
l1 <- list("crp_t0"=crp_t0, "crp_eos"=crp_eos, "crp_pod1"=crp_pod1,
  "dex_t0"=dex_t0,"dex_eos"=dex_eos,"dex_pod1"=dex_pod1 )

lapply(1:length(l1),function(i) {
  volcanoplot(x=l1[[i]], title=names(l1)[[i]])
  smearplot(x=l1[[i]], title=names(l1)[[i]])
  print( head(l1[[i]],50) |>
    kbl(caption = paste(names(l1)[[i]],"top 50 genes")) |>
    kable_paper("hover", full_width = F))
} )

## [[1]]
## 
## [[2]]
## 
## [[3]]
## 
## [[4]]
## 
## [[5]]
## 
## [[6]]
debar(l1) ; par(mar = c(5.1, 4.1, 4.1, 2.1))

myeuler(list("crp_t0"=crp_t0,"dex_t0"=dex_t0))

myeuler(list("crp_eos"=crp_eos,"dex_eos"=dex_eos))

myeuler(list("crp_pod1"=crp_pod1,"dex_pod1"=dex_pod1))

myeuler(list("crp_pod1"=crp_pod1,"dex_eos"=dex_eos))

message("crp_t0")
## crp_t0
datatable(crp_t0[1:1000,1:6])
message("crp_eos")
## crp_eos
datatable(crp_eos[1:1000,1:6])
message("crp_pod1")
## crp_pod1
datatable(crp_pod1[1:1000,1:6])
message("dex_t0")
## dex_t0
datatable(dex_t0[1:1000,1:6])
message("dex_eos")
## dex_eos
datatable(dex_eos[1:1000,1:6])
message("dex_pod1")
## dex_pod1
datatable(dex_pod1[1:1000,1:6])

Stratified, unadjusted

  • crp_t0_a
  • crp_t0_b
  • crp_eos_a
  • crp_eos_b
  • crp_pod1_a
  • crp_pod1_b
  • dex_crplo_t0
  • dex_crphi_t0
  • dex_crplo_eos
  • dex_crphi_eos
  • dex_crplo_pod1
  • dex_crphi_pod1
l2 <- list("crp_t0_a"=crp_t0_a, "crp_t0_b"=crp_t0_b,
  "crp_eos_a"=crp_eos_a, "crp_eos_b"=crp_eos_b,
  "crp_pod1_a"= crp_pod1_a, "crp_pod1_b"=crp_pod1_b,
  "dex_crplo_t0"=dex_crplo_t0, "dex_crphi_t0"= dex_crphi_t0,
  "dex_crplo_eos"=dex_crplo_eos, "dex_crphi_eos"=dex_crphi_eos,
  "dex_crplo_pod1"=dex_crplo_pod1, "dex_crphi_pod1"=dex_crphi_pod1)

lapply(1:length(l2),function(i) {
  volcanoplot(x=l2[[i]], title=names(l2)[[i]])
  smearplot(x=l2[[i]], title=names(l2)[[i]])
  print( head(l2[[i]],50)  |>
    kbl(caption = paste(names(l2)[[i]],"top 50 genes")) |>
    kable_paper("hover", full_width = F) )
} )

## [[1]]
## 
## [[2]]
## 
## [[3]]
## 
## [[4]]
## 
## [[5]]
## 
## [[6]]
## 
## [[7]]
## 
## [[8]]
## 
## [[9]]
## 
## [[10]]
## 
## [[11]]
## 
## [[12]]
debar(l2) ; par(mar = c(5.1, 4.1, 4.1, 2.1))

myeuler(list("crp_t0_a"=crp_t0_a, "crp_t0_b"=crp_t0_b,"dex_crplo_t0"=dex_crplo_t0, "dex_crphi_t0"= dex_crphi_t0))
## Warning in colSums(id & !empty) == 0 | merged_sets: longer object length is not
## a multiple of shorter object length

myeuler(list("crp_eos_a"=crp_eos_a, "crp_eos_b"=crp_eos_b,"dex_crplo_eos"=dex_crplo_eos, "dex_crphi_eos"=dex_crphi_eos))
## Warning in colSums(id & !empty) == 0 | merged_sets: longer object length is not
## a multiple of shorter object length

myeuler(list("crp_pod1_a"= crp_pod1_a, "crp_pod1_b"=crp_pod1_b,"dex_crplo_pod1"=dex_crplo_pod1, "dex_crphi_pod1"=dex_crphi_pod1))

message("crp_t0_a")
## crp_t0_a
datatable(crp_t0_a[1:1000,1:6])
message("crp_t0_b")
## crp_t0_b
datatable(crp_t0_b[1:1000,1:6])
message("crp_eos_a")
## crp_eos_a
datatable(crp_eos_a[1:1000,1:6])
message("crp_eos_b")
## crp_eos_b
datatable(crp_eos_b[1:1000,1:6])
message("crp_pod1_a")
## crp_pod1_a
datatable(crp_pod1_a[1:1000,1:6])
message("crp_pod1_b")
## crp_pod1_b
datatable(crp_pod1_b[1:1000,1:6])
message("dex_crplo_t0")
## dex_crplo_t0
datatable(dex_crplo_t0[1:1000,1:6])
message("dex_crphi_t0")
## dex_crphi_t0
datatable(dex_crphi_t0[1:1000,1:6])
message("dex_crplo_eos")
## dex_crplo_eos
datatable(dex_crplo_eos[1:1000,1:6])
message("dex_crphi_eos")
## dex_crphi_eos
datatable(dex_crphi_eos[1:1000,1:6])
message("dex_crplo_pod1")
## dex_crplo_pod1
datatable(dex_crplo_pod1[1:1000,1:6])
message("dex_crphi_pod1")
## dex_crphi_pod1
datatable(dex_crphi_pod1[1:1000,1:6])

Adjusted for clinical covariates

Unstratified, adjusted

  • crp_t0_adj
  • crp_eos_adj
  • crp_pod1_adj
  • dex_t0_adj
  • dex_eos_adj
  • dex_pod1_adj
l3 <- list("crp_t0_adj"=crp_t0_adj, "crp_eos_adj"=crp_eos_adj, "crp_pod1_adj"=crp_pod1_adj,
  "dex_t0_adj"=dex_t0_adj,"dex_eos_adj"=dex_eos_adj,"dex_pod1_adj"=dex_pod1_adj )

lapply(1:length(l3),function(i) {
  volcanoplot(x=l3[[i]], title=names(l3)[[i]])
  smearplot(x=l3[[i]], title=names(l3)[[i]])
  print( head(l1[[i]],50) |>
    kbl(caption = paste(names(l3)[[i]],"top 50 genes")) |>
    kable_paper("hover", full_width = F) )
} )

## [[1]]
## 
## [[2]]
## 
## [[3]]
## 
## [[4]]
## 
## [[5]]
## 
## [[6]]
debar(l3) ; par(mar = c(5.1, 4.1, 4.1, 2.1))

myeuler(list("crp_t0_adj"=crp_t0_adj,"dex_t0_adj"=dex_t0_adj))

myeuler(list("crp_eos_adj"=crp_eos_adj,"dex_eos_adj"=dex_eos_adj))

myeuler(list("crp_pod1_adj"=crp_pod1_adj,"dex_pod1_adj"=dex_pod1_adj))

message("crp_t0_adj")
## crp_t0_adj
datatable(crp_t0_adj[1:1000,1:6])
message("crp_eos_adj")
## crp_eos_adj
datatable(crp_eos_adj[1:1000,1:6])
message("crp_pod1_adj")
## crp_pod1_adj
datatable(crp_pod1_adj[1:1000,1:6])
message("dex_t0_adj")
## dex_t0_adj
datatable(dex_t0_adj[1:1000,1:6])
message("dex_eos_adj")
## dex_eos_adj
datatable(dex_eos_adj[1:1000,1:6])
message("dex_pod1_adj")
## dex_pod1_adj
datatable(dex_pod1_adj[1:1000,1:6])

Stratified, adjusted

  • crp_t0_a_adj
  • crp_t0_b_adj
  • crp_eos_a_adj
  • crp_eos_b_adj
  • crp_pod1_a_adj
  • crp_pod1_b_adj
  • dex_crplo_t0_adj
  • dex_crphi_t0_adj
  • dex_crplo_eos_adj
  • dex_crphi_eos_adj
  • dex_crplo_pod1_adj
  • dex_crphi_pod1_adj
l4 <- list("crp_t0_a_adj"=crp_t0_a_adj, "crp_t0_b_adj"=crp_t0_b_adj,
  "crp_eos_a_adj"=crp_eos_a_adj, "crp_eos_b_adj"=crp_eos_b_adj,
  "crp_pod1_a_adj"= crp_pod1_a_adj, "crp_pod1_b_adj"=crp_pod1_b_adj,
  "dex_crplo_t0_adj"=dex_crplo_t0_adj, "dex_crphi_t0_adj"=dex_crphi_t0_adj,
  "dex_crplo_eos_adj"=dex_crplo_eos_adj, "dex_crphi_eos_adj"=dex_crphi_eos_adj,
  "dex_crplo_pod1_adj"=dex_crplo_pod1_adj, "dex_crphi_pod1_adj"=dex_crphi_pod1_adj)

lapply(1:length(l4),function(i) {
  volcanoplot(x=l2[[i]], title=names(l2)[[i]])
  smearplot(x=l2[[i]], title=names(l2)[[i]])
  print( head(l2[[i]],50) |>
    kbl(caption = paste(names(l2)[[i]],"top 50 genes")) |>
    kable_paper("hover", full_width = F) )
} )

## [[1]]
## 
## [[2]]
## 
## [[3]]
## 
## [[4]]
## 
## [[5]]
## 
## [[6]]
## 
## [[7]]
## 
## [[8]]
## 
## [[9]]
## 
## [[10]]
## 
## [[11]]
## 
## [[12]]
debar(l4) ; par(mar = c(5.1, 4.1, 4.1, 2.1))

myeuler(list("crp_t0_a_adj"=crp_t0_a_adj, "crp_t0_b_adj"=crp_t0_b_adj,
  "dex_crplo_t0_adj"=dex_crplo_t0_adj, "dex_crphi_t0_adj"= dex_crphi_t0_adj))
## Warning in colSums(id & !empty) == 0 | merged_sets: longer object length is not
## a multiple of shorter object length

myeuler(list("crp_eos_a_adj"=crp_eos_a_adj, "crp_eos_b_adj"=crp_eos_b_adj,
  "dex_crplo_eos_adj"=dex_crplo_eos_adj, "dex_crphi_eos_adj"=dex_crphi_eos_adj))
## Warning in colSums(id & !empty) == 0 | merged_sets: longer object length is not
## a multiple of shorter object length

myeuler(list("crp_pod1_a_adj"= crp_pod1_a_adj, "crp_pod1_b_adj"=crp_pod1_b_adj,
  "dex_crplo_pod1_adj"=dex_crplo_pod1_adj, "dex_crphi_pod1_adj"=dex_crphi_pod1_adj))

message("crp_t0_a_adj")
## crp_t0_a_adj
datatable(crp_t0_a_adj[1:1000,1:6])
message("crp_t0_b_adj")
## crp_t0_b_adj
datatable(crp_t0_b_adj[1:1000,1:6])
message("crp_eos_a_adj")
## crp_eos_a_adj
datatable(crp_eos_a_adj[1:1000,1:6])
message("crp_eos_b_adj")
## crp_eos_b_adj
datatable(crp_eos_b_adj[1:1000,1:6])
message("crp_pod1_a_adj")
## crp_pod1_a_adj
datatable(crp_pod1_a_adj[1:1000,1:6])
message("crp_pod1_b_adj")
## crp_pod1_b_adj
datatable(crp_pod1_b_adj[1:1000,1:6])
message("dex_crplo_t0_adj")
## dex_crplo_t0_adj
datatable(dex_crplo_t0_adj[1:1000,1:6])
message("dex_crphi_t0_adj")
## dex_crphi_t0_adj
datatable(dex_crphi_t0_adj[1:1000,1:6])
message("dex_crplo_eos_adj")
## dex_crplo_eos_adj
datatable(dex_crplo_eos_adj[1:1000,1:6])
message("dex_crphi_eos_adj")
## dex_crphi_eos_adj
datatable(dex_crphi_eos_adj[1:1000,1:6])
message("dex_crplo_pod1_adj")
## dex_crplo_pod1_adj
datatable(dex_crplo_pod1_adj[1:1000,1:6])
message("dex_crphi_pod1_adj")
## dex_crphi_pod1_adj
datatable(dex_crphi_pod1_adj[1:1000,1:6])

Save data

writeTSV(l1)
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
writeTSV(l2)
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
## 
## [[7]]
## NULL
## 
## [[8]]
## NULL
## 
## [[9]]
## NULL
## 
## [[10]]
## NULL
## 
## [[11]]
## NULL
## 
## [[12]]
## NULL
writeTSV(l3)
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
writeTSV(l4)
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
## 
## [[7]]
## NULL
## 
## [[8]]
## NULL
## 
## [[9]]
## NULL
## 
## [[10]]
## NULL
## 
## [[11]]
## NULL
## 
## [[12]]
## NULL
system('sed -i "s#^baseMean#GeneID\tbaseMean#" *tsv')

Session information

For reproducibility

sessionInfo()
## R version 4.5.3 (2026-03-11)
## 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] 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] xlsx_0.6.5                  DT_0.34.0                  
##  [3] ggplot2_4.0.2               kableExtra_1.4.0           
##  [5] beeswarm_0.4.0              eulerr_7.0.4               
##  [7] MASS_7.3-65                 mitch_1.21.3               
##  [9] DESeq2_1.50.2               SummarizedExperiment_1.40.0
## [11] Biobase_2.70.0              MatrixGenerics_1.22.0      
## [13] matrixStats_1.5.0           GenomicRanges_1.62.1       
## [15] Seqinfo_1.0.0               IRanges_2.44.0             
## [17] S4Vectors_0.48.0            BiocGenerics_0.56.0        
## [19] generics_0.1.4              dplyr_1.2.0                
## [21] WGCNA_1.73                  fastcluster_1.3.0          
## [23] dynamicTreeCut_1.63-1       reshape2_1.4.5             
## [25] gplots_3.2.0               
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3    rstudioapi_0.17.1     jsonlite_2.0.0       
##   [4] magrittr_2.0.4        farver_2.1.2          rmarkdown_2.30       
##   [7] vctrs_0.7.1           memoise_2.0.1.9000    base64enc_0.1-6      
##  [10] htmltools_0.5.9       S4Arrays_1.10.1       SparseArray_1.10.8   
##  [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.4.6       impute_1.83.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          AnnotationDbi_1.71.0 
##  [34] crosstalk_1.2.2       textshaping_1.0.4     Hmisc_5.2-4          
##  [37] RSQLite_2.4.4         polyclip_1.10-7       httr_1.4.7           
##  [40] abind_1.4-8           compiler_4.5.3        withr_3.0.2          
##  [43] bit64_4.6.0-1         doParallel_1.0.17     htmlTable_2.4.3      
##  [46] S7_0.2.1              backports_1.5.0       BiocParallel_1.44.0  
##  [49] DBI_1.2.3             ggstats_0.11.0        DelayedArray_0.36.0  
##  [52] gtools_3.9.5          caTools_1.18.3        tools_4.5.3          
##  [55] foreign_0.8-91        otel_0.2.0            httpuv_1.6.16        
##  [58] nnet_7.3-20           glue_1.8.0            promises_1.5.0       
##  [61] polylabelr_0.3.0      grid_4.5.3            checkmate_2.3.3      
##  [64] cluster_2.1.8.2       gtable_0.3.6          preprocessCore_1.71.0
##  [67] tidyr_1.3.2           data.table_1.18.2.1   xml2_1.5.0           
##  [70] XVector_0.50.0        foreach_1.5.2         pillar_1.11.1        
##  [73] stringr_1.6.0         later_1.4.8           rJava_1.0-18         
##  [76] splines_4.5.3         lattice_0.22-7        survival_3.8-6       
##  [79] bit_4.6.0             tidyselect_1.2.1      GO.db_3.21.0         
##  [82] locfit_1.5-9.12       Biostrings_2.77.1     knitr_1.51           
##  [85] gridExtra_2.3         svglite_2.2.2         xfun_0.56            
##  [88] stringi_1.8.7         UCSC.utils_1.5.0      statnet.common_4.12.0
##  [91] yaml_2.3.12           xlsxjars_0.9.0        evaluate_1.0.5       
##  [94] codetools_0.2-20      tibble_3.3.1          cli_3.6.5            
##  [97] rpart_4.1.27          xtable_1.8-4          systemfonts_1.3.1    
## [100] jquerylib_0.1.4       network_1.19.0        dichromat_2.0-0.1    
## [103] Rcpp_1.1.1            GenomeInfoDb_1.45.4   coda_0.19-4.1        
## [106] png_0.1-8             parallel_4.5.3        blob_1.2.4           
## [109] bitops_1.0-9          viridisLite_0.4.3     scales_1.4.0         
## [112] purrr_1.2.1           crayon_1.5.3          rlang_1.1.7          
## [115] KEGGREST_1.49.0