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 name

An 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(): 2024-04-29
#> # 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           0

2 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  467550

To 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.4.0 beta (2024-04-15 r86425)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-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.16.0 ExperimentHub_2.12.0         
#>  [3] AnnotationHub_3.12.0          BiocFileCache_2.12.0         
#>  [5] dbplyr_2.5.0                  SummarizedExperiment_1.34.0  
#>  [7] Biobase_2.64.0                GenomicRanges_1.56.0         
#>  [9] GenomeInfoDb_1.40.0           IRanges_2.38.0               
#> [11] S4Vectors_0.42.0              BiocGenerics_0.50.0          
#> [13] MatrixGenerics_1.16.0         matrixStats_1.3.0            
#> [15] BiocStyle_2.32.0             
#> 
#> loaded via a namespace (and not attached):
#>  [1] KEGGREST_1.44.0         xfun_0.43               bslib_0.7.0            
#>  [4] lattice_0.22-6          vctrs_0.6.5             tools_4.4.0            
#>  [7] generics_0.1.3          curl_5.2.1              tibble_3.2.1           
#> [10] fansi_1.0.6             AnnotationDbi_1.66.0    RSQLite_2.3.6          
#> [13] highr_0.10              blob_1.2.4              pkgconfig_2.0.3        
#> [16] Matrix_1.7-0            lifecycle_1.0.4         GenomeInfoDbData_1.2.12
#> [19] compiler_4.4.0          Biostrings_2.72.0       tinytex_0.50           
#> [22] htmltools_0.5.8.1       sass_0.4.9              yaml_2.3.8             
#> [25] pillar_1.9.0            crayon_1.5.2            jquerylib_0.1.4        
#> [28] DelayedArray_0.30.0     cachem_1.0.8            magick_2.8.3           
#> [31] abind_1.4-5             mime_0.12               tidyselect_1.2.1       
#> [34] digest_0.6.35           purrr_1.0.2             dplyr_1.1.4            
#> [37] bookdown_0.39           BiocVersion_3.19.1      grid_4.4.0             
#> [40] fastmap_1.1.1           SparseArray_1.4.0       cli_3.6.2              
#> [43] magrittr_2.0.3          S4Arrays_1.4.0          utf8_1.2.4             
#> [46] withr_3.0.0             filelock_1.0.3          UCSC.utils_1.0.0       
#> [49] rappdirs_0.3.3          bit64_4.0.5             rmarkdown_2.26         
#> [52] XVector_0.44.0          httr_1.4.7              bit_4.0.5              
#> [55] png_0.1-8               memoise_2.0.1           evaluate_0.23          
#> [58] knitr_1.46              rlang_1.1.3             Rcpp_1.0.12            
#> [61] glue_1.7.0              DBI_1.2.2               BiocManager_1.30.22    
#> [64] jsonlite_1.8.8          R6_2.5.1                zlibbioc_1.50.0