1 Basics

1.1 Install awst

R is an open-source statistical environment which can be easily modified to enhance its functionality via packages. awst is a R package available via the Bioconductor repository for packages. R can be installed on any operating system from CRAN after which you can install awst by using the following commands in your R session:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("awst")

## Check that you have a valid Bioconductor installation
BiocManager::valid()

1.2 Required knowledge

awst is based on many other packages and in particular in those that have implemented the infrastructure needed for dealing with RNA-seq data. That is, packages like SummarizedExperiment.

If you are asking yourself the question “Where do I start using Bioconductor?” you might be interested in this blog post.

1.3 Asking for help

As package developers, we try to explain clearly how to use our packages and in which order to use the functions. But R and Bioconductor have a steep learning curve so it is critical to learn where to ask for help. The blog post quoted above mentions some but we would like to highlight the Bioconductor support site as the main resource for getting help: remember to use the awst tag and check the older posts. Other alternatives are available such as creating GitHub issues and tweeting. However, please note that if you want to receive help you should adhere to the posting guidelines. It is particularly critical that you provide a small reproducible example and your session information so package developers can track down the source of the error.

1.4 Citing awst

We hope that awst will be useful for your research. Please use the following information to cite the package and the overall approach. Thank you!

## Citation info
citation("awst")
#> Warning in person1(given = given[[i]], family = family[[i]], middle =
#> middle[[i]], : It is recommended to use 'given' instead of 'middle'.
#> To cite package 'awst' in publications use:
#> 
#>   Risso D, Pagnotta SM (2021). "Per-sample standardization and
#>   asymmetric winsorization lead to accurate clustering of RNA-seq
#>   expression profiles." _Bioinformatics_.
#>   doi:10.1093/bioinformatics/btab091
#>   <https://doi.org/10.1093/bioinformatics/btab091>.
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Article{,
#>     title = {Per-sample standardization and asymmetric winsorization lead to accurate clustering of RNA-seq expression profiles},
#>     author = {Davide Risso and Stefano Maria Pagnotta},
#>     year = {2021},
#>     journal = {Bioinformatics},
#>     doi = {10.1093/bioinformatics/btab091},
#>   }

2 What does awst does?

AWST aims to regularize the original read counts to reduce the effect of noise on the clustering of samples. In fact, gene expression data are characterized by high levels of noise in both lowly expressed features, which suffer from background effects and low signal-to-noise ratio, and highly expressed features, which may be the result of amplification bias and other experimental artifacts. These effects are of utmost importance in highly degraded or low input material samples, such as tumor samples and single cells.

AWST comprises two main steps. In the first one, namely the standardization step, we standardize the counts by centering and scaling them, exploiting the log-normal probability distribution. We refer to the standardized counts as z-counts. The second step, namely the smoothing step, leverages a highly skewed transformation that decreases the noise while preserving the influence of genes to separate molecular subtypes. These two steps are implemented in the awst function.

A further filtering method, implemented in the gene_filter function, is suggested to remove those features that only contribute noise to the clustering.

3 Quick start

library(awst)
library(airway)
library(SummarizedExperiment)
library(EDASeq)
library(ggplot2)

Here, we will use the data in the airway package to illustrate the awst approach.

Please, see our paper (Risso and Pagnotta, 2021) and this repository for more extensive and biologically relevant examples.

data(airway)
airway
#> class: RangedSummarizedExperiment 
#> dim: 63677 8 
#> metadata(1): ''
#> assays(1): counts
#> rownames(63677): ENSG00000000003 ENSG00000000005 ... ENSG00000273492
#>   ENSG00000273493
#> rowData names(10): gene_id gene_name ... seq_coord_system symbol
#> colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
#> colData names(9): SampleName cell ... Sample BioSample

The data are stored in a RangedSummarizedExperiment, a special case of the SummarizedExperiment class, one of the central classes in Bioconductor. If you are not familiar with it, I recomment to look at its vignette available at SummarizedExperiment.

First, we filter out non-expressed genes. For simplicity, we remove those genes with fewer than 10 reads on average across samples.

filter <- rowMeans(assay(airway)) >= 10
table(filter)
#> filter
#> FALSE  TRUE 
#> 47587 16090

se <- airway[filter,]

We are left with 16090 genes. We are now ready to apply awst to the data.

se <- awst(se)
se
#> class: RangedSummarizedExperiment 
#> dim: 16090 8 
#> metadata(1): ''
#> assays(2): counts awst
#> rownames(16090): ENSG00000000003 ENSG00000000419 ... ENSG00000273472
#>   ENSG00000273486
#> rowData names(10): gene_id gene_name ... seq_coord_system symbol
#> colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
#> colData names(9): SampleName cell ... Sample BioSample
plot(density(assay(se, "awst")[,1]), main = "Sample 1")

We can see that the majority of the values have been shrunk around −2, while the others values gradually increase up to around 4. The effect of reducing the contribution of lowly expressed genes, and of the winsorization for the highly expressed ones, results in a better separation of the samples, reflecting biological differences (Risso and Pagnotta, 2021).

The other main function of the awst package is gene_filter. It can be used to remove those genes that contribute little to nothing to the distance between samples. The function uses an entropy measure to remove the uninformative genes.

filtered <- gene_filter(se)
dim(filtered)
#> [1] 10842     8

Our final dataset is made of 8 genes.

We can see how the awst transformation leads to separation between treatment (along PC1) and cell line (along PC2).

res_pca <- prcomp(t(assay(filtered, "awst")))
df <- as.data.frame(cbind(res_pca$x, colData(airway)))
ggplot(df, aes(x = PC1, y = PC2, color = dex, shape = cell)) +
    geom_point() + theme_classic()

4 Role of normalization

Although in this example awst applied to raw data works well, a prior normalization step can help. We have found that full-quantile normalization works well and has the computational advantage of allowing awst to estimate the parameters only once for all samples (Risso and Pagnotta, 2021).

Here we show the results of awst after full-quantile normalization (implemented in EDASeq).

assay(se, "fq") <- betweenLaneNormalization(assay(se), which="full")
se <- awst(se, expr_values = "fq")

res_pca <- prcomp(t(assay(se, "awst")))
df <- as.data.frame(cbind(res_pca$x, colData(airway)))
ggplot(df, aes(x = PC1, y = PC2, color = dex, shape = cell)) +
    geom_point() + theme_classic()

5 Reproducibility

The awst package (Risso and Pagnotta, 2021) was made possible thanks to:

  • R (R Core Team, 2024)
  • BiocStyle (Oleś, 2024)
  • knitr (Xie, 2024)
  • RefManageR (McLean, 2017)
  • rmarkdown (Allaire, Xie, Dervieux, McPherson, Luraschi, Ushey, Atkins, Wickham, Cheng, Chang, and Iannone, 2024)
  • sessioninfo (Wickham, Chang, Flight, Müller, and Hester, 2021)
  • testthat (Wickham, 2011)

This package was developed using biocthis.

Code for creating the vignette

## Create the vignette
library("rmarkdown")
system.time(render("awst_intro.Rmd", "BiocStyle::html_document"))

## Extract the R code
library("knitr")
knit("awst_intro.Rmd", tangle = TRUE)

Date the vignette was generated.

#> [1] "2024-04-30 22:35:33 EDT"

Wallclock time spent generating the vignette.

#> Time difference of 21.074 secs

R session information.

#> ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.4.0 beta (2024-04-15 r86425)
#>  os       Ubuntu 22.04.4 LTS
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       America/New_York
#>  date     2024-04-30
#>  pandoc   2.7.3 @ /usr/bin/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
#>  package              * version     date (UTC) lib source
#>  abind                  1.4-5       2016-07-21 [2] CRAN (R 4.4.0)
#>  airway               * 1.23.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  AnnotationDbi          1.66.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  aroma.light            3.34.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  awst                 * 1.12.0      2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
#>  backports              1.4.1       2021-12-13 [2] CRAN (R 4.4.0)
#>  bibtex                 0.5.1       2023-01-26 [2] CRAN (R 4.4.0)
#>  Biobase              * 2.64.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  BiocFileCache          2.12.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  BiocGenerics         * 0.50.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  BiocIO                 1.14.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  BiocManager            1.30.22     2023-08-08 [2] CRAN (R 4.4.0)
#>  BiocParallel         * 1.38.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  BiocStyle            * 2.32.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  biomaRt                2.60.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  Biostrings           * 2.72.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  bit                    4.0.5       2022-11-15 [2] CRAN (R 4.4.0)
#>  bit64                  4.0.5       2020-08-30 [2] CRAN (R 4.4.0)
#>  bitops                 1.0-7       2021-04-24 [2] CRAN (R 4.4.0)
#>  blob                   1.2.4       2023-03-17 [2] CRAN (R 4.4.0)
#>  bookdown               0.39        2024-04-15 [2] CRAN (R 4.4.0)
#>  bslib                  0.7.0       2024-03-29 [2] CRAN (R 4.4.0)
#>  cachem                 1.0.8       2023-05-01 [2] CRAN (R 4.4.0)
#>  cli                    3.6.2       2023-12-11 [2] CRAN (R 4.4.0)
#>  codetools              0.2-20      2024-03-31 [3] CRAN (R 4.4.0)
#>  colorspace             2.1-0       2023-01-23 [2] CRAN (R 4.4.0)
#>  crayon                 1.5.2       2022-09-29 [2] CRAN (R 4.4.0)
#>  curl                   5.2.1       2024-03-01 [2] CRAN (R 4.4.0)
#>  DBI                    1.2.2       2024-02-16 [2] CRAN (R 4.4.0)
#>  dbplyr                 2.5.0       2024-03-19 [2] CRAN (R 4.4.0)
#>  DelayedArray           0.30.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  deldir                 2.0-4       2024-02-28 [2] CRAN (R 4.4.0)
#>  digest                 0.6.35      2024-03-11 [2] CRAN (R 4.4.0)
#>  dplyr                  1.1.4       2023-11-17 [2] CRAN (R 4.4.0)
#>  EDASeq               * 2.38.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  evaluate               0.23        2023-11-01 [2] CRAN (R 4.4.0)
#>  fansi                  1.0.6       2023-12-08 [2] CRAN (R 4.4.0)
#>  farver                 2.1.1       2022-07-06 [2] CRAN (R 4.4.0)
#>  fastmap                1.1.1       2023-02-24 [2] CRAN (R 4.4.0)
#>  filelock               1.0.3       2023-12-11 [2] CRAN (R 4.4.0)
#>  generics               0.1.3       2022-07-05 [2] CRAN (R 4.4.0)
#>  GenomeInfoDb         * 1.40.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  GenomeInfoDbData       1.2.12      2024-04-16 [2] Bioconductor
#>  GenomicAlignments    * 1.40.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  GenomicFeatures        1.56.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  GenomicRanges        * 1.56.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  ggplot2              * 3.5.1       2024-04-23 [2] CRAN (R 4.4.0)
#>  glue                   1.7.0       2024-01-09 [2] CRAN (R 4.4.0)
#>  gtable                 0.3.5       2024-04-22 [2] CRAN (R 4.4.0)
#>  highr                  0.10        2022-12-22 [2] CRAN (R 4.4.0)
#>  hms                    1.1.3       2023-03-21 [2] CRAN (R 4.4.0)
#>  htmltools              0.5.8.1     2024-04-04 [2] CRAN (R 4.4.0)
#>  httr                   1.4.7       2023-08-15 [2] CRAN (R 4.4.0)
#>  httr2                  1.0.1       2024-04-01 [2] CRAN (R 4.4.0)
#>  hwriter                1.3.2.1     2022-04-08 [2] CRAN (R 4.4.0)
#>  interp                 1.1-6       2024-01-26 [2] CRAN (R 4.4.0)
#>  IRanges              * 2.38.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  jpeg                   0.1-10      2022-11-29 [2] CRAN (R 4.4.0)
#>  jquerylib              0.1.4       2021-04-26 [2] CRAN (R 4.4.0)
#>  jsonlite               1.8.8       2023-12-04 [2] CRAN (R 4.4.0)
#>  KEGGREST               1.44.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  knitr                  1.46        2024-04-06 [2] CRAN (R 4.4.0)
#>  labeling               0.4.3       2023-08-29 [2] CRAN (R 4.4.0)
#>  lattice                0.22-6      2024-03-20 [3] CRAN (R 4.4.0)
#>  latticeExtra           0.6-30      2022-07-04 [2] CRAN (R 4.4.0)
#>  lifecycle              1.0.4       2023-11-07 [2] CRAN (R 4.4.0)
#>  lubridate              1.9.3       2023-09-27 [2] CRAN (R 4.4.0)
#>  magrittr               2.0.3       2022-03-30 [2] CRAN (R 4.4.0)
#>  Matrix                 1.7-0       2024-03-22 [3] CRAN (R 4.4.0)
#>  MatrixGenerics       * 1.16.0      2024-04-30 [2] Bioconductor 3.19 (R 4.4.0)
#>  matrixStats          * 1.3.0       2024-04-11 [2] CRAN (R 4.4.0)
#>  memoise                2.0.1       20