HighlyReplicatedRNASeq
2023-10-26
The HighlyReplicatedRNASeq package provides functions to access the count matrix from bulk RNA-seq studies with many replicates. For example,the study from Schurch et al. (2016) has data on 86 samples of S. cerevisiae in two conditions.
1 Load Data
To load the dataset, call the Schurch16() function. It returns a SummarizedExperiment:
schurch_se <- HighlyReplicatedRNASeq::Schurch16()
#> see ?HighlyReplicatedRNASeq and browseVignettes('HighlyReplicatedRNASeq') for documentation
#> loading from cache
#> see ?HighlyReplicatedRNASeq and browseVignettes('HighlyReplicatedRNASeq') for documentation
#> loading from cache
schurch_se
#> class: SummarizedExperiment
#> dim: 7126 86
#> metadata(0):
#> assays(1): counts
#> rownames(7126): 15S_rRNA 21S_rRNA ... tY(GUA)O tY(GUA)Q
#> rowData names(0):
#> colnames(86): wildtype_01 wildtype_02 ... knockout_47 knockout_48
#> colData names(4): file_name condition replicate nameAn alternative approach that achieves exactly the same is to load the data directly from ExperimentHub
library(ExperimentHub)
eh <- ExperimentHub()
records <- query(eh, "HighlyReplicatedRNASeq")
records[1] ## display the metadata for the first resource
#> ExperimentHub with 1 record
#> # snapshotDate(): 2023-10-24
#> # names(): EH3315
#> # package(): HighlyReplicatedRNASeq
#> # $dataprovider: Geoff Barton's group on GitHub
#> # $species: Saccharomyces cerevisiae BY4741
#> # $rdataclass: matrix
#> # $rdatadateadded: 2020-04-03
#> # $title: Schurch S. cerevesiae Highly Replicated Bulk RNA-Seq Counts
#> # $description: Count matrix for bulk RNA-sequencing dataset from 86 S. cere...
#> # $taxonomyid: 1247190
#> # $genome: Ensembl release 68
#> # $sourcetype: tar.gz
#> # $sourceurl: https://github.com/bartongroup/profDGE48
#> # $sourcesize: NA
#> # $tags: c("ExperimentHub", "ExperimentData", "ExpressionData",
#> # "SequencingData", "RNASeqData")
#> # retrieve record with 'object[["EH3315"]]'
count_matrix <- records[["EH3315"]] ## load the count matrix by ID
#> see ?HighlyReplicatedRNASeq and browseVignettes('HighlyReplicatedRNASeq') for documentation
#> loading from cache
count_matrix[1:10, 1:5]
#> wildtype_01 wildtype_02 wildtype_03 wildtype_04 wildtype_05
#> 15S_rRNA 2 12 31 8 21
#> 21S_rRNA 20 76 101 99 128
#> HRA1 3 2 2 2 3
#> ICR1 75 123 107 157 98
#> LSR1 60 163 233 163 193
#> NME1 13 14 23 13 29
#> PWR1 0 0 0 0 0
#> Q0010 0 0 0 0 0
#> Q0017 0 0 0 0 0
#> Q0032 0 0 0 0 02 Explore Data
It has 7126 genes and 86 samples. The counts are between 0 and 467,000.
summary(c(assay(schurch_se, "counts")))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0 89 386 1229 924 467550To make the data easier to work with, I will “normalize” the data. First I divide it by the mean of each sample to account for the differential sequencing depth. Then, I apply the log() transformation and add a small number to avoid taking the logarithm of 0.
norm_counts <- assay(schurch_se, "counts")
norm_counts <- log(norm_counts / colMeans(norm_counts) + 0.001)The histogram of the transformed data looks very smooth:
hist(norm_counts, breaks = 100)
Finally, let us take a look at the MA-plot of the data and the volcano plot:
wt_mean <- rowMeans(norm_counts[, schurch_se$condition == "wildtype"])
ko_mean <- rowMeans(norm_counts[, schurch_se$condition == "knockout"])
plot((wt_mean+ ko_mean) / 2, wt_mean - ko_mean,
pch = 16, cex = 0.4, col = "#00000050", frame.plot = FALSE)
abline(h = 0)
pvalues <- sapply(seq_len(nrow(norm_counts)), function(idx){
tryCatch(
t.test(norm_counts[idx, schurch_se$condition == "wildtype"],
norm_counts[idx, schurch_se$condition == "knockout"])$p.value,
error = function(err) NA)
})
plot(wt_mean - ko_mean, - log10(pvalues),
pch = 16, cex = 0.5, col = "#00000050", frame.plot = FALSE)
3 Reference
- Schurch, N. J., Schofield, P., Gierliński, M., Cole, C., Sherstnev, A., Singh, V., … Barton, G. J. (2016). How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? Rna, 22(6), 839–851. https://doi.org/10.1261/rna.053959.115
4 Session Info
sessionInfo()
#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
#> 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_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
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] HighlyReplicatedRNASeq_1.14.0 ExperimentHub_2.10.0
#> [3] AnnotationHub_3.10.0 BiocFileCache_2.10.0
#> [5] dbplyr_2.3.4 SummarizedExperiment_1.32.0
#> [7] Biobase_2.62.0 GenomicRanges_1.54.0
#> [9] GenomeInfoDb_1.38.0 IRanges_2.36.0
#> [11] S4Vectors_0.40.0 BiocGenerics_0.48.0
#> [13] MatrixGenerics_1.14.0 matrixStats_1.0.0
#> [15] BiocStyle_2.30.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.0 dplyr_1.1.3
#> [3] blob_1.2.4 filelock_1.0.2
#> [5] Biostrings_2.70.1 bitops_1.0-7
#> [7] fastmap_1.1.1 RCurl_1.98-1.12
#> [9] promises_1.2.1 digest_0.6.33
#> [11] mime_0.12 lifecycle_1.0.3
#> [13] ellipsis_0.3.2 KEGGREST_1.42.0
#> [15] interactiveDisplayBase_1.40.0 RSQLite_2.3.1
#> [17] magrittr_2.0.3 compiler_4.3.1
#> [19] rlang_1.1.1 sass_0.4.7
#> [21] tools_4.3.1 utf8_1.2.4
#> [23] yaml_2.3.7 knitr_1.44
#> [25] S4Arrays_1.2.0 bit_4.0.5
#> [27] curl_5.1.0 DelayedArray_0.28.0
#> [29] abind_1.4-5 withr_2.5.1
#> [31] purrr_1.0.2 grid_4.3.1
#> [33] fansi_1.0.5 xtable_1.8-4
#> [35] cli_3.6.1 rmarkdown_2.25
#> [37] crayon_1.5.2 generics_0.1.3
#> [39] httr_1.4.7 DBI_1.1.3
#> [41] cachem_1.0.8 zlibbioc_1.48.0
#> [43] AnnotationDbi_1.64.0 BiocManager_1.30.22
#> [45] XVector_0.42.0 vctrs_0.6.4
#> [47] Matrix_1.6-1.1 jsonlite_1.8.7
#> [49] bookdown_0.36 bit64_4.0.5
#> [51] magick_2.8.1 jquerylib_0.1.4
#> [53] glue_1.6.2 BiocVersion_3.18.0
#> [55] later_1.3.1 tibble_3.2.1
#> [57] pillar_1.9.0 rappdirs_0.3.3
#> [59] htmltools_0.5.6.1 GenomeInfoDbData_1.2.11
#> [61] R6_2.5.1 evaluate_0.22
#> [63] shiny_1.7.5.1 lattice_0.22-5
#> [65] png_0.1-8 memoise_2.0.1
#> [67] httpuv_1.6.12 bslib_0.5.1
#> [69] Rcpp_1.0.11 SparseArray_1.2.0
#> [71] xfun_0.40 pkgconfig_2.0.3