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
  • avb_t0
  • avb_eos
  • avb_pod1
l1 <- list("crp_t0"=crp_t0, "crp_eos"=crp_eos, "crp_pod1"=crp_pod1,
  "avb_t0"=avb_t0,"avb_eos"=avb_eos,"avb_pod1"=avb_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,"avb_t0"=avb_t0))

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

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

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("avb_t0")
## avb_t0
datatable(avb_t0[1:1000,1:6])
message("avb_eos")
## avb_eos
datatable(avb_eos[1:1000,1:6])
message("avb_pod1")
## avb_pod1
datatable(avb_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
  • avb_crplo_t0
  • avb_crphi_t0
  • avb_crplo_eos
  • avb_crphi_eos
  • avb_crplo_pod1
  • avb_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,
  "avb_crplo_t0"=avb_crplo_t0, "avb_crphi_t0"= avb_crphi_t0,
  "avb_crplo_eos"=avb_crplo_eos, "avb_crphi_eos"=avb_crphi_eos,
  "avb_crplo_pod1"=avb_crplo_pod1, "avb_crphi_pod1"=avb_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,"avb_crplo_t0"=avb_crplo_t0, "avb_crphi_t0"= avb_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,"avb_crplo_eos"=avb_crplo_eos, "avb_crphi_eos"=avb_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,"avb_crplo_pod1"=avb_crplo_pod1, "avb_crphi_pod1"=avb_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("avb_crplo_t0")
## avb_crplo_t0
datatable(avb_crplo_t0[1:1000,1:6])
message("avb_crphi_t0")
## avb_crphi_t0
datatable(avb_crphi_t0[1:1000,1:6])
message("avb_crplo_eos")
## avb_crplo_eos
datatable(avb_crplo_eos[1:1000,1:6])
message("avb_crphi_eos")
## avb_crphi_eos
datatable(avb_crphi_eos[1:1000,1:6])
message("avb_crplo_pod1")
## avb_crplo_pod1
datatable(avb_crplo_pod1[1:1000,1:6])
message("avb_crphi_pod1")
## avb_crphi_pod1
datatable(avb_crphi_pod1[1:1000,1:6])

Adjusted for clinical covariates

Unstratified, adjusted

  • crp_t0_adj
  • crp_eos_adj
  • crp_pod1_adj
  • avb_t0_adj
  • avb_eos_adj
  • avb_pod1_adj
l3 <- list("crp_t0_adj"=crp_t0_adj, "crp_eos_adj"=crp_eos_adj, "crp_pod1_adj"=crp_pod1_adj,
  "avb_t0_adj"=avb_t0_adj,"avb_eos_adj"=avb_eos_adj,"avb_pod1_adj"=avb_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,"avb_t0_adj"=avb_t0_adj))

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

myeuler(list("crp_pod1_adj"=crp_pod1_adj,"avb_pod1_adj"=avb_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("avb_t0_adj")
## avb_t0_adj
datatable(avb_t0_adj[1:1000,1:6])
message("avb_eos_adj")
## avb_eos_adj
datatable(avb_eos_adj[1:1000,1:6])
message("avb_pod1_adj")
## avb_pod1_adj
datatable(avb_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
  • avb_crplo_t0_adj
  • avb_crphi_t0_adj
  • avb_crplo_eos_adj
  • avb_crphi_eos_adj
  • avb_crplo_pod1_adj
  • avb_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,
  "avb_crplo_t0_adj"=avb_crplo_t0_adj, "avb_crphi_t0_adj"= avb_crphi_t0_adj,
  "avb_crplo_eos_adj"=avb_crplo_eos_adj, "avb_crphi_eos_adj"=avb_crphi_eos_adj,
  "avb_crplo_pod1_adj"=avb_crplo_pod1_adj, "avb_crphi_pod1_adj"=avb_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,
  "avb_crplo_t0_adj"=avb_crplo_t0_adj, "avb_crphi_t0_adj"= avb_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,
  "avb_crplo_eos_adj"=avb_crplo_eos_adj, "avb_crphi_eos_adj"=avb_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,
  "avb_crplo_pod1_adj"=avb_crplo_pod1_adj, "avb_crphi_pod1_adj"=avb_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("avb_crplo_t0_adj")
## avb_crplo_t0_adj
datatable(avb_crplo_t0_adj[1:1000,1:6])
message("avb_crphi_t0_adj")
## avb_crphi_t0_adj
datatable(avb_crphi_t0_adj[1:1000,1:6])
message("avb_crplo_eos_adj")
## avb_crplo_eos_adj
datatable(avb_crplo_eos_adj[1:1000,1:6])
message("avb_crphi_eos_adj")
## avb_crphi_eos_adj
datatable(avb_crphi_eos_adj[1:1000,1:6])
message("avb_crplo_pod1_adj")
## avb_crplo_pod1_adj
datatable(avb_crplo_pod1_adj[1:1000,1:6])
message("avb_crphi_pod1_adj")
## avb_crphi_pod1_adj
datatable(avb_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.4.1 (2024-06-14)
## 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=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: 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.33                    
##  [3] eulerr_7.0.2                ggplot2_3.5.1              
##  [5] kableExtra_1.4.0            MASS_7.3-61                
##  [7] mitch_1.17.4                DESeq2_1.44.0              
##  [9] SummarizedExperiment_1.34.0 Biobase_2.64.0             
## [11] MatrixGenerics_1.16.0       matrixStats_1.4.1          
## [13] GenomicRanges_1.56.1        GenomeInfoDb_1.40.1        
## [15] IRanges_2.38.1              S4Vectors_0.42.1           
## [17] BiocGenerics_0.50.0         dplyr_1.1.4                
## [19] WGCNA_1.73                  fastcluster_1.2.6          
## [21] dynamicTreeCut_1.63-1       reshape2_1.4.4             
## [23] gplots_3.2.0               
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3      rstudioapi_0.17.1       jsonlite_1.8.9         
##   [4] magrittr_2.0.3          rmarkdown_2.28          zlibbioc_1.50.0        
##   [7] vctrs_0.6.5             memoise_2.0.1.9000      base64enc_0.1-3        
##  [10] htmltools_0.5.8.1       S4Arrays_1.4.1          SparseArray_1.4.8      
##  [13] Formula_1.2-5           sass_0.4.9              bslib_0.8.0            
##  [16] KernSmooth_2.23-24      htmlwidgets_1.6.4       plyr_1.8.9             
##  [19] echarts4r_0.4.5         impute_1.78.0           cachem_1.1.0           
##  [22] mime_0.12               lifecycle_1.0.4         iterators_1.0.14       
##  [25] pkgconfig_2.0.3         Matrix_1.7-0            R6_2.5.1               
##  [28] fastmap_1.2.0           GenomeInfoDbData_1.2.12 shiny_1.9.1            
##  [31] digest_0.6.37           colorspace_2.1-1        GGally_2.2.1           
##  [34] AnnotationDbi_1.66.0    crosstalk_1.2.1         Hmisc_5.1-3            
##  [37] RSQLite_2.3.7           fansi_1.0.6             polyclip_1.10-7        
##  [40] httr_1.4.7              abind_1.4-8             compiler_4.4.1         
##  [43] withr_3.0.1             bit64_4.5.2             doParallel_1.0.17      
##  [46] htmlTable_2.4.3         backports_1.5.0         BiocParallel_1.38.0    
##  [49] DBI_1.2.3               ggstats_0.7.0           highr_0.11             
##  [52] DelayedArray_0.30.1     gtools_3.9.5            caTools_1.18.3         
##  [55] tools_4.4.1             foreign_0.8-87          beeswarm_0.4.0         
##  [58] httpuv_1.6.15           nnet_7.3-19             glue_1.8.0             
##  [61] promises_1.3.0          polylabelr_0.2.0        grid_4.4.1             
##  [64] checkmate_2.3.2         cluster_2.1.6           generics_0.1.3         
##  [67] gtable_0.3.5            preprocessCore_1.66.0   tidyr_1.3.1            
##  [70] data.table_1.16.0       xml2_1.3.6              utf8_1.2.4             
##  [73] XVector_0.44.0          foreach_1.5.2           pillar_1.9.0           
##  [76] stringr_1.5.1           later_1.3.2             rJava_1.0-11           
##  [79] splines_4.4.1           lattice_0.22-6          survival_3.7-0         
##  [82] bit_4.5.0               tidyselect_1.2.1        GO.db_3.19.1           
##  [85] locfit_1.5-9.10         Biostrings_2.72.1       knitr_1.48             
##  [88] gridExtra_2.3           svglite_2.1.3           xfun_0.48              
##  [91] stringi_1.8.4           UCSC.utils_1.0.0        yaml_2.3.10            
##  [94] xlsxjars_0.6.1          evaluate_1.0.1          codetools_0.2-20       
##  [97] tibble_3.2.1            cli_3.6.3               rpart_4.1.23           
## [100] xtable_1.8-4            systemfonts_1.1.0       jquerylib_0.1.4        
## [103] munsell_0.5.1           Rcpp_1.0.13             png_0.1-8              
## [106] parallel_4.4.1          blob_1.2.4              bitops_1.0-9           
## [109] viridisLite_0.4.2       scales_1.3.0            purrr_1.0.2            
## [112] crayon_1.5.3            rlang_1.1.4             KEGGREST_1.44.1