Exploratory analysis on

1. Environment Set Up

knitr::opts_chunk$set(echo = TRUE, warning=FALSE)
library(RNASeqBulkExploratory)  #Our Package
library(DT)
library(gridExtra)
library(SummarizedExperiment)
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
library(ggplot2)
SE_Bio <- readRDS(params$SEBio)
CpmFilt <- params$CpmFilt
SampleFilt <- params$SampleFilt
OutputFolder <- params$OutputFolder

2. Data Upload

SE_Filt <- filterSE(SE_Bio, cpmTh=CpmFilt, sampleTh=SampleFilt)
## UPDATED library.size slot
  • Re-calculate library size and calculate TMM
  • Calculate cpm and logCpm (with normalization) and store in SE_Filt
SE_Filt <- normTmmSE(SE_Filt, useNormFactors=TRUE, priorCount=0.25)

3. Sample-to-sample correlation heatmap

Correlation matrix across samples calculated on the basis of the Spearman correlation. Heatmap annotation takes into consideration sample ‘Condition’ by default.

# Set heatmap size on the basis of the number of samples
Width <- 7 + (dim(SE_Filt)[2]- 8) * 0.24
Width <- ifelse(Width<=14, Width, 14)
Height <- 5 + (dim(SE_Filt)[2]- 8) * 0.22
Cols <- c('DMSO' = 'grey30', 'CTL' = 'azure3', 
                 'AhHyd_Ag'='#F8766D', 'AhHyd_Inh'='#F8766D50',
                 'Andr_Ag'='#fccb17', 'Andr_Inh'='#C49A0050',  
                 "Estr_Ag"= '#53B400', "Estr_Inh"= '#53B40050', 
                 'GC_Ag' = '#00C094', 'GC_Inh' = '#00C09450',
                 'LivX_Ag' = '#00B6EB', 'LivX_Inh' = '#00B6EB50', 
                 'Ret_Ag' = '#A58AFF', 'Ret_Inh' = '#A58AFF50', 
                 'Thyr_Ag' = '#FB61D7', 'Thyr_Inh' = '#FB61D750'
                 )

colData(SE_Filt)$SeqRun <- factor(colData(SE_Filt)$SeqRun,
                                  levels = c("20210310", "20210724", "20220422"),
                                  labels = c("batch1_run1", "batch1_run2", "batch2_run3"))

seqrun_cols <- c(
  "batch1_run1" = "#0072B2",  
  "batch1_run2" = "#D55E00",  
  "batch2_run3" = "#CC79A7"   
)
cond_order <- c("CTL", "DMSO")
other_cond <- setdiff(unique(colData(SE_Filt)$Condition), cond_order)
cond_levels <- c(cond_order, sort(other_cond))
colData(SE_Filt)$Condition <- factor(
  colData(SE_Filt)$Condition,
  levels = cond_levels
)
ord <- order(colData(SE_Filt)$Condition)

SE_Filt_ord <- SE_Filt[, ord]
sampleCorrHeatmap2 <- function(
    SE,
    plotTitle = NULL,
    annotation_col = NULL,
    annotation_colors = NULL,
    myPalette = NULL,
    display_numbers = TRUE,
    ...
) {
    if (is.null(plotTitle)) {
        plotTitle <- deparse(substitute(SE))
    }

    if (is.null(annotation_col)) {
        ann_fields <- match("Condition", names(colData(SE)))
        if (is.na(ann_fields)) {
            warning("The SE object does not contain the default 'Condition' annotation field")
            ann_fields <- integer(0)
        }
    } else {
        ann_fields <- match(annotation_col, names(colData(SE)))
        if (any(is.na(ann_fields))) {
            warning(
                paste(
                    "WARNING: missing annotation field(s):",
                    paste(annotation_col[is.na(ann_fields)], collapse = ", ")
                )
            )
            ann_fields <- ann_fields[!is.na(ann_fields)]
        }
    }

    annotation_df <- as.data.frame(colData(SE)[ann_fields])

    if (!("logcpm" %in% names(assays(SE)))) {
        stop('SE object has no "logcpm" assay, please use normTmmSE() to compute it')
    }

    SampleCorr <- stats::cor(assays(SE)$logcpm, method = "spearman")

    if (is.null(myPalette)) {
        GreenCols <- grDevices::colorRampPalette(
            RColorBrewer::brewer.pal(7, "Greens")
        )(255)
    } else {
        GreenCols <- myPalette
    }

    if (ncol(SE) > 20) {
        display_numbers <- FALSE
        message('There are many samples. For readability "display_numbers" will be set to FALSE')
    }

    pheatmap::pheatmap(
        as.matrix(SampleCorr),
        main = paste(plotTitle, "\n Spearman correlation"),
        col = GreenCols,
        annotation_col = annotation_df,
        annotation_colors = annotation_colors,
        fontsize_col = 6,
        fontsize_row = 6,
        clustering_method = "average",
        display_numbers = display_numbers,
        ...
    )
}
sampleCorrHeatmap2(
  SE_Filt_ord,
  plotTitle = "Filtered dataset",
  annotation_col = c("SeqRun", "Condition"),
  annotation_colors = list(
    Condition = Cols,
    SeqRun = seqrun_cols),
  display_numbers = TRUE,
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  show_rownames = FALSE,
  show_colnames = FALSE
)
## There are many samples. For readability "display_numbers" will be set to FALSE

ggsave(filename=paste0(OutputFolder, 'PAPER_CORRPLOTS/SampleCorr_CTL08_v2.pdf'), 
       plot=sampleCorrHeatmap2(
  SE_Filt_ord,
  plotTitle = "Filtered dataset",
  annotation_col = c("SeqRun", "Condition"),
  annotation_colors = list(
    Condition = Cols,
    SeqRun = seqrun_cols),
  display_numbers = TRUE,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  show_rownames = FALSE,
  show_colnames = FALSE
),
       device='pdf', width=5, height=5)
## There are many samples. For readability "display_numbers" will be set to FALSE

4. Principal Component Analysis Plot

pcaSE(SE_Filt, PlotTitle='First-Second Component', components=c(1,2), Interactive=TRUE, colVec = Cols)
#SE_Filt@colData@listData[["SeqRun"]] <- as.character(SE_Filt@colData@listData[["SeqRun"]])
pcaSE(SE_Filt, PlotTitle='First-Second Component', components=c(1,2), condition = "SeqRun", colVec = seqrun_cols)
pcaSE_new <- function(SE, PlotTitle = NULL, colVec = NULL, 
                  components = c(1, 2), condition = "Condition", 
                  Interactive = TRUE) 
{
    PCA_list <- calculatePCA(SE, condition = condition)

    if (is.null(PlotTitle)) {
        PlotTitle = deparse(substitute(SE))
    }

    PCData <- PCA_list$PCAData[, c(components, 5, 6)]
    names(PCData) <- c("x_PC", "y_PC", "Sample", "Condition")

    sub <- "Principal Component Analysis"
    pcx <- paste("PC", components[1], "- VarExp: ", round(PCA_list$Var[components[1]], 2))
    pcy <- paste("PC", components[2], "- VarExp: ", round(PCA_list$Var[components[2]], 2))

    if (is.null(colVec)) {
        colVec = (scales::hue_pal())(length(unique(colData(SE)[[condition]])))
    }

    PCAPl <- ggplot(PCData, aes(x = .data$x_PC, y = .data$y_PC)) +
        geom_point(aes(col = .data$Condition), size = 3) +
        labs(x = pcx, y = pcy, title = PlotTitle, subtitle = sub) +
        scale_colour_manual(values = colVec, name = NULL) +
        theme_bw() +
        theme(plot.title = element_text(face = "bold", colour = "darkred", size = 14, hjust = 0.5),
              plot.subtitle = element_text(colour = "darkred", size = 12, hjust = 0.5),
              legend.text = element_text(size = 8),
              legend.position = "bottom")

    if (Interactive == TRUE) {
        PCAPl <- plotly::ggplotly(PCAPl, tooltip = c("Sample", "Condition")) %>%
            plotly::layout(title = list(text = paste0(PlotTitle, "<br>", "<sup>", sub)))
        return(PCAPl)
    } else {
        return(PCAPl) 
    }
}
ggsave(filename=paste0(OutputFolder, 'PAPER_CORRPLOTS/PCA_run_CTL08.pdf'),
       plot=pcaSE_new(SE_Filt, PlotTitle='First-Second Component', components=c(1,2), condition = "SeqRun", colVec = seqrun_cols, Interactive = FALSE), device='pdf', width=6, height=5)
SessionInfo <- sessionInfo()
Date <- date()
SessionInfo
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.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       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggplot2_3.4.1               SummarizedExperiment_1.28.0
##  [3] Biobase_2.58.0              GenomicRanges_1.50.2       
##  [5] GenomeInfoDb_1.34.9         IRanges_2.32.0             
##  [7] S4Vectors_0.36.1            BiocGenerics_0.44.0        
##  [9] MatrixGenerics_1.10.0       matrixStats_0.63.0         
## [11] gridExtra_2.3               DT_0.27                    
## [13] RNASeqBulkExploratory_0.2.1
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.10            locfit_1.5-9.7         lattice_0.20-45       
##  [4] tidyr_1.3.0            digest_0.6.31          utf8_1.2.3            
##  [7] R6_2.5.1               evaluate_0.20          httr_1.4.5            
## [10] highr_0.10             pillar_1.8.1           zlibbioc_1.44.0       
## [13] rlang_1.1.1            lazyeval_0.2.2         data.table_1.14.8     
## [16] rstudioapi_0.14        jquerylib_0.1.4        Matrix_1.5-3          
## [19] rmarkdown_2.20         labeling_0.4.2         textshaping_0.3.6     
## [22] htmlwidgets_1.6.1      pheatmap_1.0.12        RCurl_1.98-1.10       
## [25] munsell_0.5.0          DelayedArray_0.24.0    compiler_4.2.1        
## [28] xfun_0.37              pkgconfig_2.0.3        systemfonts_1.0.4     
## [31] htmltools_0.5.4        tidyselect_1.2.0       tibble_3.2.1          
## [34] GenomeInfoDbData_1.2.9 edgeR_3.40.2           viridisLite_0.4.1     
## [37] fansi_1.0.4            dplyr_1.1.0            withr_2.5.0           
## [40] bitops_1.0-7           grid_4.2.1             jsonlite_1.8.4        
## [43] gtable_0.3.1           lifecycle_1.0.3        magrittr_2.0.3        
## [46] scales_1.2.1           cli_3.6.1              cachem_1.0.7          
## [49] farver_2.1.1           XVector_0.38.0         limma_3.54.1          
## [52] bslib_0.4.2            ellipsis_0.3.2         ragg_1.2.3            
## [55] generics_0.1.3         vctrs_0.6.2            RColorBrewer_1.1-3    
## [58] tools_4.2.1            glue_1.6.2             purrr_1.0.1           
## [61] crosstalk_1.2.0        fastmap_1.1.1          yaml_2.3.7            
## [64] colorspace_2.1-0       plotly_4.10.1          knitr_1.42            
## [67] sass_0.4.5
Date
## [1] "Tue Apr 28 13:10:47 2026"