Introduction to awst
21 September 2021
Package
awst 1.1.1
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")
#> It is recommended to use 'given' instead of 'middle'.
#>
#> 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 (URL:
#> 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: 64102 8
#> metadata(1): ''
#> assays(1): counts
#> rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
#> rowData names(0):
#> colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
#> colData names(9): SampleName cell ... Sample BioSampleThe 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
#> 48012 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(0):
#> 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 8Our 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, 2021)
- BiocStyle
- knitr (Xie, 2021)
- RefManageR (McLean, 2017)
- rmarkdown (Allaire, Xie, McPherson, Luraschi, Ushey, Atkins, Wickham, Cheng, Chang, and Iannone, 2021)
- sessioninfo (Csárdi, core, Wickham, Chang, Flight, Müller, and Hester, 2018)
- 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] "2021-09-21 19:26:40 EDT"Wallclock time spent generating the vignette.
#> Time difference of 14.06 secsR session information.
#> - Session info -------------------------------------------------------------------------------------------------------
#> setting value
#> version R version 4.1.1 (2021-08-10)
#> os Windows Server x64
#> system x86_64, mingw32
#> ui RTerm
#> language (EN)
#> collate C
#> ctype English_United States.1252
#> tz America/New_York
#> date 2021-09-21
#>
#> - Packages -----------------------------------------------------------------------------------------------------------
#> package * version date lib source
#> airway * 1.13.0 2021-05-22 [2] Bioconductor
#> AnnotationDbi 1.55.1 2021-09-21 [2] Bioconductor
#> aroma.light 3.23.1 2021-09-21 [2] Bioconductor
#> assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.1.0)
#> awst * 1.1.1 2021-09-21 [1] Bioconductor
#> Biobase * 2.53.0 2021-09-21 [2] Bioconductor
#> BiocFileCache 2.1.1 2021-09-21 [2] Bioconductor
#> BiocGenerics * 0.39.2 2021-09-21 [2] Bioconductor
#> BiocIO 1.3.0 2021-09-21 [2] Bioconductor
#> BiocManager 1.30.16 2021-06-15 [2] CRAN (R 4.1.0)
#> BiocParallel * 1.27.9 2021-09-21 [2] Bioconductor
#> BiocStyle * 2.21.3 2021-09-21 [2] Bioconductor
#> biomaRt 2.49.4 2021-09-21 [2] Bioconductor
#> Biostrings * 2.61.2 2021-09-21 [2] Bioconductor
#> bit 4.0.4 2020-08-04 [2] CRAN (R 4.1.0)
#> bit64 4.0.5 2020-08-30 [2] CRAN (R 4.1.0)
#> bitops 1.0-7 2021-04-24 [2] CRAN (R 4.1.0)
#> blob 1.2.2 2021-07-23 [2] CRAN (R 4.1.0)
#> bookdown 0.24 2021-09-02 [2] CRAN (R 4.1.1)
#> bslib 0.3.0 2021-09-02 [2] CRAN (R 4.1.1)
#> cachem 1.0.6 2021-08-19 [2] CRAN (R 4.1.1)
#> cli 3.0.1 2021-07-17 [2] CRAN (R 4.1.0)
#> colorspace 2.0-2 2021-06-24 [2] CRAN (R 4.1.0)
#> crayon 1.4.1 2021-02-08 [2] CRAN (R 4.1.0)
#> curl 4.3.2 2021-06-23 [2] CRAN (R 4.1.0)
#> DBI 1.1.1 2021-01-15 [2] CRAN (R 4.1.0)
#> dbplyr 2.1.1 2021-04-06 [2] CRAN (R 4.1.0)
#> DelayedArray 0.19.3 2021-09-21 [2] Bioconductor
#> digest 0.6.27 2020-10-24 [2] CRAN (R 4.1.0)
#> dplyr 1.0.7 2021-06-18 [2] CRAN (R 4.1.0)
#> EDASeq * 2.27.1 2021-09-21 [2] Bioconductor
#> ellipsis 0.3.2 2021-04-29 [2] CRAN (R 4.1.0)
#> evaluate 0.14 2019-05-28 [2] CRAN (R 4.1.0)
#> fansi 0.5.0 2021-05-25 [2] CRAN (R 4.1.0)
#> farver 2.1.0 2021-02-28 [2] CRAN (R 4.1.0)
#> fastmap 1.1.0 2021-01-25 [2] CRAN (R 4.1.0)
#> filelock 1.0.2 2018-10-05 [2] CRAN (R 4.1.0)
#> generics 0.1.0 2020-10-31 [2] CRAN (R 4.1.0)
#> GenomeInfoDb * 1.29.8 2021-09-21 [2] Bioconductor
#> GenomeInfoDbData 1.2.6 2021-08-25 [2] Bioconductor
#> GenomicAlignments * 1.29.0 2021-09-21 [2] Bioconductor
#> GenomicFeatures 1.45.2 2021-09-21 [2] Bioconductor
#> GenomicRanges * 1.45.0 2021-09-21 [2] Bioconductor
#> ggplot2 * 3.3.5 2021-06-25 [2] CRAN (R 4.1.0)
#> glue 1.4.2 2020-08-27 [2] CRAN (R 4.1.0)
#> gtable 0.3.0 2019-03-25 [2] CRAN (R 4.1.0)
#> highr 0.9 2021-04-16 [2] CRAN (R 4.1.0)
#> hms 1.1.0 2021-05-17 [2] CRAN (R 4.1.0)
#> htmltools 0.5.2 2021-08-25 [2] CRAN (R 4.1.1)
#> httr 1.4.2 2020-07-20 [2] CRAN (R 4.1.0)
#> hwriter 1.3.2 2014-09-10 [2] CRAN (R 4.1.0)
#> IRanges * 2.27.2 2021-09-21 [2] Bioconductor
#> jpeg 0.1-9 2021-07-24 [2] CRAN (R 4.1.0)
#> jquerylib 0.1.4 2021-04-26 [2] CRAN (R 4.1.0)
#> jsonlite 1.7.2 2020-12-09 [2] CRAN (R 4.1.0)
#> KEGGREST 1.33.0 2021-09-21 [2] Bioconductor
#> knitr 1.34 2021-09-09 [2] CRAN (R 4.1.1)
#> labeling 0.4.2 2020-10-20 [2] CRAN (R 4.1.0)
#> lattice 0.20-44 2021-05-02 [2] CRAN (R 4.1.1)
#> latticeExtra 0.6-29 2019-12-19 [2] CRAN (R 4.1.0)
#> lifecycle 1.0.0 2021-02-15 [2] CRAN (R 4.1.0)
#> lubridate 1.7.10 2021-02-26 [2] CRAN (R 4.1.0)
#> magrittr 2.0.1 2020-11-17 [2] CRAN (R 4.1.0)
#> Matrix 1.3-4 2021-06-01 [2] CRAN (R 4.1.1)
#> MatrixGenerics * 1.5.4 2021-09-21 [2] Bioconductor
#> matrixStats * 0.61.0 2021-09-17 [2] CRAN (R 4.1.1)
#> memoise 2.0.0 2021-01-26 [2] CRAN (R 4.1.0)
#> munsell 0.5.0 2018-06-12 [2] CRAN (R 4.1.0)
#> pillar 1.6.2 2021-07-29 [2] CRAN (R 4.1.0)
#> pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.1.0)
#> plyr 1.8.6 2020-03-03 [2] CRAN (R 4.1.0)
#> png 0.1-7 2013-12-03 [2] CRAN (R 4.1.0)
#> prettyunits 1.1.1 2020-01-24 [2] CRAN (R 4.1.0)
#> progress 1.2.2 2019-05-16 [2] CRAN (R 4.1.0)
#> purrr 0.3.4 2020-04-17 [2] CRAN (R 4.1.0)
#> R.methodsS3 1.8.1 2020-08-26 [2] CRAN (R 4.1.0)
#> R.oo 1.24.0 2020-08-26 [2] CRAN (R 4.1.0)
#> R.utils 2.10.1 2020-08-26 [2] CRAN (R 4.1.0)
#> R6 2.5.1 2021-08-19 [2] CRAN (R 4.1.1)
#> rappdirs 0.3.3 2021-01-31 [2] CRAN (R 4.1.0)
#> RColorBrewer 1.1-2 2014-12-07 [2] CRAN (R 4.1.0)
#> Rcpp 1.0.7 2021-07-07 [2] CRAN (R 4.1.0)
#> RCurl 1.98-1.5 2021-09-17 [2] CRAN (R 4.1.1)
#> RefManageR * 1.3.0 2020-11-13 [2] CRAN (R 4.1.0)
#> restfulr 0.0.13 2017-08-06 [2] CRAN (R 4.1.0)
#> rjson 0.2.20 2018-06-08 [2] CRAN (R 4.1.0)
#> rlang 0.4.11 2021-04-30 [2] CRAN (R 4.1.0)
#> rmarkdown 2.11 2021-09-14 [2] CRAN (R 4.1.1)
#> Rsamtools * 2.9.1 2021-09-21 [2] Bioconductor
#> RSQLite 2.2.8 2021-08-21 [2] CRAN (R 4.1.1)
#> rtracklayer 1.53.1 2021-09-21 [2] Bioconductor
#> S4Vectors * 0.31.3 2021-09-21 [2] Bioconductor
#> sass 0.4.0 2021-05-12 [2] CRAN (R 4.1.0)
#> scales 1.1.1 2020-05-11 [2] CRAN (R 4.1.0)
#> sessioninfo * 1.1.1 2018-11-05 [2] CRAN (R 4.1.0)
#> ShortRead * 1.51.0 2021-09-21 [2] Bioconductor
#> stringi 1.7.4 2021-08-25 [2] CRAN (R 4.1.1)
#> stringr 1.4.0 2019-02-10 [2] CRAN (R 4.1.0)
#> SummarizedExperiment * 1.23.4 2021-09-21 [2] Bioconductor
#> tibble 3.1.4 2021-08-25 [2] CRAN (R 4.1.1)
#> tidyselect 1.1.1 2021-04-30 [2] CRAN (R 4.1.0)
#> utf8 1.2.2 2021-07-24 [2] CRAN (R 4.1.0)
#> vctrs 0.3.8 2021-04-29 [2] CRAN (R 4.1.0)
#> withr 2.4.2 2021-04-18 [2] CRAN (R 4.1.0)
#> xfun 0.26 2021-09-14 [2] CRAN (R 4.1.1)
#> XML 3.99-0.8 2021-09-17 [2] CRAN (R 4.1.1)
#> xml2 1.3.2 2020-04-23 [2] CRAN (R 4.1.0)
#> XVector * 0.33.0 2021-09-21 [2] Bioconductor
#> yaml 2.2.1 2020-02-01 [2] CRAN (R 4.1.0)
#> zlibbioc 1.39.0 2021-09-21 [2] Bioconductor
#>
#> [1] D:/biocbuild/bbs-3.14-bioc/tmpdir/RtmpCSXNXs/Rinst2560626752
#> [2] D:/biocbuild/bbs-3.14-bioc/R/library6 Bibliography
This vignette was generated using BiocStyle
with knitr (Xie, 2021) and rmarkdown (Allaire, Xie, McPherson, et al., 2021) running behind the scenes.
Citations made with RefManageR (McLean, 2017).
[1] J. Allaire, Y. Xie, J. McPherson, et al. rmarkdown: Dynamic Documents for R. R package version 2.11. 2021. URL: https://github.com/rstudio/rmarkdown.
[2] G. Csárdi, R. core, H. Wickham, et al. sessioninfo: R Session Information. R package version 1.1.1. 2018. URL: https://CRAN.R-project.org/package=sessioninfo.
[3] M. W. McLean. “RefManageR: Import and Manage BibTeX and BibLaTeX References in R”. In: The Journal of Open Source Software (2017). DOI: 10.21105/joss.00338.
[4] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria, 2021. URL: https://www.R-project.org/.
[5] D. Risso and S. M. Pagnotta. “Per-sample standardization and asymmetric winsorization lead to accurate clustering of RNA-seq expression profiles”. In: Bioinformatics (2021). DOI: 10.1093/bioinformatics/btab091.
[6] H. Wickham. “testthat: Get Started with Testing”. In: The R Journal 3 (2011), pp. 5–10. URL: https://journal.r-project.org/archive/2011-1/RJournal_2011-1_Wickham.pdf.
[7] Y. Xie. knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.34. 2021. URL: https://yihui.org/knitr/.