countsimQC - Comparing characteristic features across count data sets

Charlotte Soneson

2022-02-03

Introduction

The countsimQC package provides a simple way to compare the characteristic features of a collection of (e.g., RNA-seq) count data sets. An important application is in situations where a synthetic count data set has been generated using a real count data set as an underlying source of parameters, in which case it is often important to verify that the final synthetic data captures the main features of the original data set. However, the package can be used to create a visual overview of any collection of one or more count data sets.

Data preparation

In this vignette we will show how to generate a comparative report from a collection of two simulated data sets and the original, underlying real data set. First, we load the object containing the three data sets. The object is a named list, where each element is a DESeqDataSet object, containing the count matrix, a sample information data frame and a model formula (necessary to calculate dispersions). For more information about the DESeqDataSet class, please see the DESeq2 Bioconductor package. For speed reasons, we use only a subset of the features in each data set for the following calculations.

suppressPackageStartupMessages({
  library(countsimQC)
  library(DESeq2)
})

data(countsimExample)
countsimExample
## $Original
## class: DESeqDataSet 
## dim: 10000 11 
## metadata(1): version
## assays(1): counts
## rownames(10000): ENSMUSG00000000001.4 ENSMUSG00000000028.14 ...
##   ENSMUSG00000048027.7 ENSMUSG00000048029.10
## rowData names(0):
## colnames(11): GSM1923445 GSM1923446 ... GSM1923578 GSM1923579
## colData names(2): group sample
## 
## $Sim1
## class: DESeqDataSet 
## dim: 10000 11 
## metadata(1): version
## assays(1): counts
## rownames(10000): Gene1 Gene2 ... Gene9999 Gene10000
## rowData names(0):
## colnames(11): Cell1 Cell2 ... Cell88 Cell89
## colData names(4): Cell Batch Group ExpLibSize
## 
## $Sim2
## class: DESeqDataSet 
## dim: 10000 11 
## metadata(1): version
## assays(1): counts
## rownames(10000): Gene1 Gene2 ... Gene9999 Gene10000
## rowData names(0):
## colnames(11): Cell1 Cell2 ... Cell88 Cell89
## colData names(3): Cell CellFac Group
countsimExample <- lapply(countsimExample, function(cse) {
  cse[seq_len(1500), ]
})

Report generation

Next, we generate the report using the countsimQCReport() function. Depending on the level of detail and the type of information that are required for the final report, this function can be run in different “modes”:

Here, for the sake of speed, we calculate statistics for a small subset of the observations (subsampleSize = 25) and refrain from calculating permutation p-values.

tempDir <- tempdir()
countsimQCReport(ddsList = countsimExample, outputFile = "countsim_report.html",
                 outputDir = tempDir, outputFormat = "html_document", 
                 showCode = FALSE, forceOverwrite = TRUE,
                 savePlots = TRUE, description = "This is my test report.", 
                 maxNForCorr = 25, maxNForDisp = Inf, 
                 calculateStatistics = TRUE, subsampleSize = 25,
                 kfrac = 0.01, kmin = 5, 
                 permutationPvalues = FALSE, nPermutations = NULL)

The countsimQCReport() function can generate either an HTML file (by setting outputFormat = "html_document" or outputFormat = NULL) or a pdf file (by setting outputFormat = "pdf_document"). The description argument can be used to provide a more extensive description of the data set(s) that are included in the report.

Generation of individual figures

If the argument savePlots is set to TRUE, an .rds file containing the individual ggplot objects will be generated. These objects can be used to perform fine-tuning of the visualizations if desired. Note, however, that the .rds file can become large if the number of data sets is large, or if the individual data sets have many samples or features. The convenience function generateIndividualPlots() can be used to quickly generate individual figures for all plots included in the report, using a variety of devices. For example, to generate each plot in pdf format:

ggplots <- readRDS(file.path(tempDir, "countsim_report_ggplots.rds"))
if (!dir.exists(file.path(tempDir, "figures"))) {
  dir.create(file.path(tempDir, "figures"))
}
generateIndividualPlots(ggplots, device = "pdf", nDatasets = 3, 
                        outputDir = file.path(tempDir, "figures"))
## `geom_smooth()` using formula 'y ~ x'

Input data format

In the example above, all data sets were provided as DESeqDataSet objects. The advantage of this is that it allows the specification of the experimental design, which is used in the dispersion calculations. countsimQC also allows a data set to be provided as either a data.frame or a matrix. However, in these situations, it will be assumed that all samples are replicates (i.e., a design ~1). An example is provided in the countsimExample_dfmat data set, provided with the package.

data(countsimExample_dfmat)
names(countsimExample_dfmat)
## [1] "Original" "Sim1"     "Sim2"
lapply(countsimExample_dfmat, class)
## $Original
## [1] "DESeqDataSet"
## attr(,"package")
## [1] "DESeq2"
## 
## $Sim1
## [1] "matrix" "array" 
## 
## $Sim2
## [1] "data.frame"
tempDir <- tempdir()
countsimQCReport(ddsList = countsimExample_dfmat, 
                 outputFile = "countsim_report_dfmat.html",
                 outputDir = tempDir, outputFormat = "html_document", 
                 showCode = FALSE, forceOverwrite = TRUE,
                 savePlots = TRUE, description = "This is my test report.", 
                 maxNForCorr = 25, maxNForDisp = Inf, 
                 calculateStatistics = TRUE, subsampleSize = 25,
                 kfrac = 0.01, kmin = 5, 
                 permutationPvalues = FALSE, nPermutations = NULL)

Session info

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [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] DESeq2_1.34.0               SummarizedExperiment_1.24.0
##  [3] Biobase_2.54.0              MatrixGenerics_1.6.0       
##  [5] matrixStats_0.61.0          GenomicRanges_1.46.1       
##  [7] GenomeInfoDb_1.30.1         IRanges_2.28.0             
##  [9] S4Vectors_0.32.3            BiocGenerics_0.40.0        
## [11] countsimQC_1.12.1          
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-155           bitops_1.0-7           bit64_4.0.5           
##  [4] RColorBrewer_1.1-2     httr_1.4.2             tools_4.1.2           
##  [7] bslib_0.3.1            utf8_1.2.2             R6_2.5.1              
## [10] DT_0.20                mgcv_1.8-38            DBI_1.1.2             
## [13] colorspace_2.0-2       withr_2.4.3            tidyselect_1.1.1      
## [16] bit_4.0.4              compiler_4.1.2         textshaping_0.3.6     
## [19] cli_3.1.1              randtests_1.0          DelayedArray_0.20.0   
## [22] labeling_0.4.2         sass_0.4.0             caTools_1.18.2        
## [25] scales_1.1.1           genefilter_1.76.0      systemfonts_1.0.3     
## [28] stringr_1.4.0          digest_0.6.29          rmarkdown_2.11        
## [31] XVector_0.34.0         pkgconfig_2.0.3        htmltools_0.5.2       
## [34] highr_0.9              fastmap_1.1.0          limma_3.50.0          
## [37] htmlwidgets_1.5.4      rlang_1.0.1            RSQLite_2.2.9         
## [40] jquerylib_0.1.4        generics_0.1.2         farver_2.1.0          
## [43] jsonlite_1.7.3         crosstalk_1.2.0        BiocParallel_1.28.3   
## [46] dplyr_1.0.7            RCurl_1.98-1.5         magrittr_2.0.2        
## [49] GenomeInfoDbData_1.2.7 Matrix_1.4-0           Rcpp_1.0.8            
## [52] munsell_0.5.0          fansi_1.0.2            lifecycle_1.0.1       
## [55] stringi_1.7.6          yaml_2.2.2             edgeR_3.36.0          
## [58] zlibbioc_1.40.0        grid_4.1.2             blob_1.2.2            
## [61] parallel_4.1.2         crayon_1.4.2           lattice_0.20-45       
## [64] Biostrings_2.62.0      splines_4.1.2          annotate_1.72.0       
## [67] KEGGREST_1.34.0        locfit_1.5-9.4         knitr_1.37            
## [70] pillar_1.7.0           geneplotter_1.72.0     XML_3.99-0.8          
## [73] glue_1.6.1             evaluate_0.14          png_0.1-7             
## [76] vctrs_0.3.8            gtable_0.3.0           purrr_0.3.4           
## [79] tidyr_1.2.0            assertthat_0.2.1       cachem_1.0.6          
## [82] ggplot2_3.3.5          xfun_0.29              xtable_1.8-4          
## [85] ragg_1.2.1             survival_3.2-13        tibble_3.1.6          
## [88] AnnotationDbi_1.56.2   memoise_2.0.1          ellipsis_0.3.2