vignettes/countsimQC.Rmd
countsimQC.RmdThe 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.
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 GroupNext, 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”:
calculateStatistics = FALSE, only plots will
be generated. This is the fastest way of running
countsimQCReport(), and in many cases generates enough
information for the user to make a visual evaluation of the count data
set(s).calculateStatistics = TRUE and
permutationPvalues = FALSE, some quantitative pairwise
comparisons between data sets will be performed. In particular, the
Kolmogorov-Smirnov test and the Wald-Wolfowitz runs test will be used to
compare distributions, and additional statistics will be calculated to
evaluate how similar the evaluated aspects are between pairs of data
sets.calculateStatistics = TRUE and
permutationPvalues = TRUE (and giving the requested number
of permutations via the nPermutations argument),
permutation of data set labels will be used to evaluate the significance
of the statistics calculated in the previous point. Naturally, this
increases the run time of the analysis considerably.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.
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'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.
## [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)## R version 4.4.0 (2024-04-24)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.4.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: UTC
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DESeq2_1.43.5 SummarizedExperiment_1.33.3
## [3] Biobase_2.63.1 MatrixGenerics_1.15.1
## [5] matrixStats_1.3.0 GenomicRanges_1.55.4
## [7] GenomeInfoDb_1.39.14 IRanges_2.37.1
## [9] S4Vectors_0.41.7 BiocGenerics_0.49.1
## [11] countsimQC_1.23.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.2 bitops_1.0-7 rlang_1.1.3
## [4] magrittr_2.0.3 compiler_4.4.0 RSQLite_2.3.6
## [7] mgcv_1.9-1 png_0.1-8 systemfonts_1.0.6
## [10] vctrs_0.6.5 pkgconfig_2.0.3 crayon_1.5.2
## [13] fastmap_1.1.1 XVector_0.43.1 labeling_0.4.3
## [16] caTools_1.18.2 utf8_1.2.4 rmarkdown_2.26
## [19] UCSC.utils_0.99.7 ragg_1.3.0 purrr_1.0.2
## [22] bit_4.0.5 xfun_0.43 zlibbioc_1.49.3
## [25] cachem_1.0.8 jsonlite_1.8.8 blob_1.2.4
## [28] highr_0.10 DelayedArray_0.29.9 BiocParallel_1.37.1
## [31] parallel_4.4.0 R6_2.5.1 bslib_0.7.0
## [34] limma_3.59.10 genefilter_1.85.1 jquerylib_0.1.4
## [37] Rcpp_1.0.12 knitr_1.46 Matrix_1.7-0
## [40] splines_4.4.0 tidyselect_1.2.1 abind_1.4-5
## [43] yaml_2.3.8 codetools_0.2-20 randtests_1.0.2
## [46] lattice_0.22-6 tibble_3.2.1 withr_3.0.0
## [49] KEGGREST_1.43.1 evaluate_0.23 desc_1.4.3
## [52] survival_3.6-4 Biostrings_2.71.6 pillar_1.9.0
## [55] DT_0.33 generics_0.1.3 ggplot2_3.5.1
## [58] munsell_0.5.1 scales_1.3.0 xtable_1.8-4
## [61] glue_1.7.0 tools_4.4.0 annotate_1.81.2
## [64] locfit_1.5-9.9 fs_1.6.4 XML_3.99-0.16.1
## [67] grid_4.4.0 tidyr_1.3.1 crosstalk_1.2.1
## [70] AnnotationDbi_1.65.2 edgeR_4.1.31 colorspace_2.1-0
## [73] nlme_3.1-164 GenomeInfoDbData_1.2.12 cli_3.6.2
## [76] textshaping_0.3.7 fansi_1.0.6 S4Arrays_1.3.7
## [79] dplyr_1.1.4 gtable_0.3.5 sass_0.4.9
## [82] digest_0.6.35 SparseArray_1.3.7 htmlwidgets_1.6.4
## [85] farver_2.1.1 memoise_2.0.1 htmltools_0.5.8.1
## [88] pkgdown_2.0.9.9000 lifecycle_1.0.4 httr_1.4.7
## [91] statmod_1.5.0 bit64_4.0.5