diff --git a/.Rbuildignore b/.Rbuildignore index c46ab48e..0f56b95e 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,3 +5,5 @@ .RData .pdf inst/extras/ +^doc$ +^Meta$ diff --git a/.gitignore b/.gitignore index b424f209..2bab5502 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,6 @@ inst/doc docs .Rhistory *.RData +.DS_Store +/doc/ +/Meta/ diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 00000000..75245367 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,349 @@ +# Contributing to this R package + +## Etiquette and conduct + +Contributors and maintainers are expected to abide by the +[Bioconductor Code of Conduct](https://bioconductor.org/about/code-of-conduct/) + +This instructions have been adapted from + +- [(Bioconductor contributor guidelines)](https://github.com/Bioconductor/Contributions/blob/master/CONTRIBUTING.md) +- [(BioJulia contributor guidelines)](https://github.com/BioJulia/Contributing) + + +## How can I contribute? + + +### Reporting Bugs + +Following these guidelines will help us to better understand your +report, replicate the issues, and track down the source problems. + +#### Before creating a bug report: + +Please do the following: + +1. Check the GitHub issues for this package + +2. If you find an open issue that describes the same problem, kindly add a comment to let + others know that you are experiencing the same issue. + +3. If no **open** issue already exists for your problem + then kindly create a new issue. + + > **Note:** If you find a **Closed** issue that seems like it is the same thing + > that you're experiencing, open a new issue and include a link to the original + > issue in the body of the new one. + +#### How to create a (good) new bug report: + +Bugs are tracked as [GitHub issues](https://guides.github.com/features/issues/). + +When you are creating a bug report, please do the following: + +1. **Explain the problem** + + - *Use a clear and descriptive title* for the issue to identify the problem. + - *Describe the exact steps which reproduce the problem* in as many details as possible. + - *Provide a specific example*. (Includes links to pastebin, gists and so on.) + If you're providing snippets in the issue, use + [Markdown code blocks](https://help.github.com/articles/markdown-basics/#multiple-lines). + + - *Describe the behaviour you observed after following the steps* + - Point out what exactly is the problem with that behaviour. + - *Explain which behaviour you expected to see instead and why.* + - *OPTIONALLY: Include screenshots and animated GIFs* which show you + following the described steps. + +2. **Provide additional context for the problem (some of these may not always apply)** + + - *Did the problem start happening recently* (e.g. after updating to a new version)? + - If the problem started recently, *can you reproduce the problem in older versions?* + - Do you know the most recent package version in which the problem doesn't happen? + + - *Can you reliably reproduce the issue?* If not... + - Provide details about how often the problem happens. + - Provide details about under which conditions it normally happens. + + - Is the problem is related to *working with files*? If so.... + - Does the problem happen for all files and projects or only some? + - Does the problem happen only when working with local or remote files? + - Does the problem happen for files of a specific type, size, or encoding? + - Is there anything else special about the files you are using? + +3. **Include details about your configuration and environment** + +- Packages and versions you are using? (see R command sessionInfo()) + +- Name and version of the OS you're using? + + +### Suggesting enhancements + +Enhancements include new features as well as minor improvements to +existing functionality. Following these suggestions will help +maintainers and the community understand your suggestion and find +related suggestions. + +#### Before Submitting An Enhancement Proposal + +1. **Perform a cursory issue search** to see if the enhancement has already been suggested. +2. If it has not, open a new issue as per the guidance below. +3. If it has... + 1. Add a comment to the existing issue instead of opening a new one. + 2. If it was closed, take the time to understand why this was so (it's ok to + ask! :) ), and consider whether anything has changed that makes the reason + outdated. If you can think of a convincing reason to reconsider the + enhancement, feel free to open a new issue as per the guidance below. + +#### How to submit a (good) new enhancement proposal + +Enhancement proposals are tracked as +[GitHub issues](https://guides.github.com/features/issues/). + +1. **Explain the enhancement** + - *Use a clear and descriptive title* for the issue to identify the suggestion. + - *Provide a step-by-step description of the suggested enhancement* in as many details as possible. + - *Provide specific examples to demonstrate the steps*. + Include copy/pasteable snippets which you use in those examples, as + [Markdown code blocks](https://help.github.com/articles/markdown-basics/#multiple-lines). + + - If you want to change current behaviour... + - Describe the *current* behaviour. + - *Explain which behaviour you expected* to see instead and *why*. + - *Will the proposed change alter APIs or existing exposed methods/types?* + If so, this may cause dependency issues and breakages, so the maintainer + will need to consider this when versioning the next release. + - *OPTIONALLY: Include screenshots and animated GIFs*. + +2. **Provide additional context for the enhancement** + + - *Explain why this enhancement would be useful* to users and + isn't something that can or should be implemented as a separate package. + + - *Do you know of other projects where this enhancement exists?* + +3. **Include details about your configuration and environment** + + - Package versions (see R sessionInfo()). + + - Name and version of the OS + + +### Making Pull Requests + + +All packages can be developed locally, and contributions suggested via pull requests. + +Let's use a git flow kind of approach. Development version should be done +against the `master` branch and then merged to `release` for release. +(https://guides.github.com/introduction/flow/) + +Before you start working on code, it is often a good idea to open an enhancement +[suggestion](#suggest-an-enhancement) + +Hit the 'Fork' button on the repositories page to create a forked copy +of the target package for your own Github account. This will ensure +your work and experiments won't hinder other users of the released and +stable package. + +From there you can clone your fork of the package and work on it on your +machine using git. +Here's an example of cloning, assuming you already forked `miaViz`: + +```sh +git clone https://github.com//miaViz.git +``` + +Git will download or "clone" your fork and put it in your local folder. + +It is beyond the scope of this document to describe good git and github use in +more specific detail, as the folks at Git and GitHub have already done that wonderfully +on their own sites. + +If you have additional questions, +see the contact details at [microbiome.github.io](microbiome.github.io). + +#### How to make (good) code contributions and new Pull-Requests + +1. **In your code changes** + + - **Branch properly!** + - If you are making a bug-fix, then you need to checkout your bug-fix branch + from the last release tag. + - If you are making a feature addition or other enhancement, checkout your + branch from master. + - See [here](#a-suggested-branching-model) for more information (or ask a package maintainer :smile:). + + - Please comment liberally for complex pieces of internal code to facilitate comprehension. + +2. **In your pull request** + + - *Describe* the changes in the pull request + + - Provide a *clear, simple, descriptive title*. + + - Do not include issue numbers in the PR title. + + - If you have implemented *new features* or behaviour + - *Provide a description of the addition* in as many details as possible. + - *Provide justification of the addition*. + - *Provide a runnable example of use of your addition*. This lets reviewers + and others try out the feature before it is merged or makes it's way to release. + + - If you have *changed current behaviour*... + - *Describe the behaviour prior to you changes* + - *Describe the behaviour after your changes* and justify why you have made the changes. + - *Does your change alter APIs or existing exposed methods/types?* + If so, this may cause dependency issues and breakages, so the maintainer + will need to consider this when versioning the next release. + - If you are implementing changes that are intended to increase performance, you + should provide the results of a simple performance benchmark exercise + demonstrating the improvement. Especially if the changes make code less legible. + +#### Reviews and merging + +You can open a pull request early on and push changes to it until it is ready, +or you can do all your editing locally and make a pull request only when it is +finished - it is up to you. + +When your pull request is ready on Github, +mention one of the maintainers of the repo in a comment e.g. `@antagomir` +and ask them to review it. +You can also use Github's review feature. +They will review the code and documentation in the pull request, +and will assess it. + +Your pull request will be accepted and merged if: + +1. The dedicated package maintainers approve the pull request for merging. +2. The automated build system confirms that all unit tests pass without any issues. + +It may also be that the reviewers or package maintainers will want to +you to make changes to your pull request before they will merge it. +Take the time to understand why any such request has been made, and +freely discuss it with the reviewers. + +## Styleguides + +### Git Commit messages + +* Use the present tense ("Add feature" not "Added feature"). +* Use the imperative mood ("Move cursor to..." not "Moves cursor to..."). +* Limit the first line to 72 characters or less. +* Reference issues and pull requests liberally after the first line. +* Consider starting the commit message with an applicable emoji: + * :art: `:art:` when improving the format/structure of the code + * :racehorse: `:racehorse:` when improving performance + * :memo: `:memo:` when writing docs + * :penguin: `:penguin:` when fixing something on Linux + * :apple: `:apple:` when fixing something on macOS + * :checkered_flag: `:checkered_flag:` when fixing something on Windows + * :bug: `:bug:` when fixing a bug + * :fire: `:fire:` when removing code or files + * :green_heart: `:green_heart:` when fixing the CI build + * :white_check_mark: `:white_check_mark:` when adding tests + * :arrow_up: `:arrow_up:` when upgrading dependencies + * :arrow_down: `:arrow_down:` when downgrading dependencies + * :exclamation: `:exclamation:` when removing warnings or depreciations + +### Additional style suggestions + +- Indent with 4 spaces. + +- Separate logical blocks of code with one blank line. Although it is common + and acceptable for short single-line functions to be defined together on + consecutive lines with no blank lines between them. + +- Function names, apart from constructors, are all camelCase. + +- Generally try to keep lines below 80-columns, unless splitting a long line + onto multiple lines makes it harder to read. + + +## How to develop the package + +A minimal workflow for package development is as follows. + +1. _Fork_ (with Git) miaViz R package from [https://github.com/microbiome/miaViz] + +1. _Clone_ (with Git) your own fork of miaViz; or if you have cloned it + earlier, then make sure that your branch is in sync with the + _upstream_ repository (see next section) + +1. Test package (in R): `devtools::test()` (and fix possible issues) + +1. Test examples: `devtools::run_examples()` + +1. Update documentation (convert R files to man/ files): `devtools::document()` + +1. Build package: `devtools::build()` + +1. Run Bioconductor checks: `BiocCheck::BiocCheck()` (for more details on Bioconductor formatting guidelines and automated checks, see [Bioconductor package guidelines](https://bioconductor.org/packages/devel/bioc/vignettes/BiocCheck/inst/doc/BiocCheck.html) and instructions on [Bioconductor packages: Development, Maintenance, and Peer Review](https://contributions.bioconductor.org/)) + +1. Load the updated package: `devtools::load_all()` + +1. Install the package: `devtools::install()` + +1. Once the package passes all build checks, commit & push to your own + fork, then make a pull request (PR) to the original + repository. This will trigger additional automated checks. If they + fail, inspect the reason, fix accordingly, and update your PR + (simply commit + push new material to your fork; the PR is + already open). + + +Final checks after the R tests; run on command line (replace "R" with +your custom R path if necessary): + +- `R CMD build `miaViz/` +- `R CMD check miaViz_x.y.z.tar.gz` +- `R CMD BiocCheck miaViz_x.y.z.tar.gz` +- `R CMD INSTALL miaViz_x.y.z.tar.gz` + + +After updating the package, remember to do the following (if applicable): + +- update documentation +- update unit tests (in `tests/`) +- update vignettes (in `vignettes`) +- run all checks from above + +After accepted pull request, check if further updates are needed in: + +- OMA +- other related packages (e.g. mia, miaSim, miaTime) + + +## Sync with the official Github development version + +Before new pull request, check that your own fork is in sync with the +main development version. Command line version goes as follows. + +``` +# Add the main repository as your _upstream_ +`git remote add upstream git@github.com:microbiome/miaViz.git` + +# Fetch the changes from upstream (main version) and origin (your fork): +git fetch --all + +# Merge those changes to your fork: +git merge origin/master +git merge upstream/master + +# Finally, commit and push to origin (your version) +# then open a new PR from the Github site if relevant. +git push origin master +``` + +## Syncing with the official Bioconductor version + +This is usually done only by the main developer. In this case, the upstream will be set to Bioconductor: `git remote add upstream git@git.bioconductor.org:packages/miaViz.git` + +- [Workflow to sync between github and bioc versions](https://bioconductor.org/developers/how-to/git/sync-existing-repositories/) + +- [Bioconductor Build reports (devel)](https://bioconductor.org/checkResults/devel/bioc-LATEST/) + +- [Package page in Bioconductor](https://bioconductor.org/packages/miaViz/) + + diff --git a/DESCRIPTION b/DESCRIPTION index 95e74db4..1c3551eb 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: miaViz Title: Microbiome Analysis Plotting and Visualization -Version: 1.9.3 +Version: 1.9.5 Authors@R: c(person(given = "Tuomas", family = "Borman", role = c("aut", "cre"), email = "tuomas.v.borman@utu.fi", @@ -12,6 +12,9 @@ Authors@R: email = "leo.lahti@iki.fi", comment = c(ORCID = "0000-0001-5537-637X")), person(given = "Basil", family = "Courbayre", role = c("ctb")), + person(given = "Giulio", family = "Benedetti", role = c("ctb"), + email = "giulio.benedetti@utu.fi", + comment = c(ORCID = "0000-0002-8732-7692")), person(given = "Muluh", family = "Muluh", role=c("ctb")) ) Description: The miaViz package implements functions to visualize @@ -47,7 +50,9 @@ Imports: tidyr, dplyr, ape, - DirichletMultinomial + DirichletMultinomial, + ggrepel, + SingleCellExperiment Suggests: knitr, rmarkdown, @@ -56,8 +61,7 @@ Suggests: patchwork, vegan, microbiomeDataSets, - bluster, - miaTime + bluster Roxygen: list(markdown = TRUE) RoxygenNote: 7.2.3 VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index c4a449db..387a8e02 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,11 +1,13 @@ # Generated by roxygen2: do not edit by hand export(plotAbundanceDensity) +export(plotCCA) export(plotColGraph) export(plotColTile) export(plotDMNFit) export(plotPrevalence) export(plotPrevalentAbundance) +export(plotRDA) export(plotRowGraph) export(plotRowTile) export(plotSeries) @@ -16,12 +18,14 @@ exportMethods(colTreeData) exportMethods(combineTreeData) exportMethods(plotAbundance) exportMethods(plotAbundanceDensity) +exportMethods(plotCCA) exportMethods(plotColGraph) exportMethods(plotColTile) exportMethods(plotColTree) exportMethods(plotDMNFit) exportMethods(plotPrevalence) exportMethods(plotPrevalentAbundance) +exportMethods(plotRDA) exportMethods(plotRowGraph) exportMethods(plotRowTile) exportMethods(plotRowTree) @@ -45,6 +49,8 @@ importFrom(DelayedArray,rowMeans) importFrom(DelayedArray,rowSums) importFrom(DirichletMultinomial,mixture) importFrom(S4Vectors,metadata) +importFrom(SingleCellExperiment,reducedDim) +importFrom(SingleCellExperiment,reducedDimNames) importFrom(SummarizedExperiment,assay) importFrom(SummarizedExperiment,colData) importFrom(SummarizedExperiment,rowData) @@ -90,6 +96,8 @@ importFrom(ggplot2,scale_y_discrete) importFrom(ggplot2,theme) importFrom(ggplot2,theme_bw) importFrom(ggplot2,theme_classic) +importFrom(ggrepel,geom_label_repel) +importFrom(ggrepel,geom_text_repel) importFrom(ggtree,geom_cladelab) importFrom(ggtree,geom_highlight) importFrom(ggtree,geom_nodelab) @@ -105,6 +113,7 @@ importFrom(rlang,"!!") importFrom(rlang,":=") importFrom(rlang,sym) importFrom(scater,plotExpression) +importFrom(scater,plotReducedDim) importFrom(scater,retrieveCellInfo) importFrom(scater,retrieveFeatureInfo) importFrom(stats,sd) diff --git a/NEWS b/NEWS index f51dd3fa..f30f55fa 100644 --- a/NEWS +++ b/NEWS @@ -20,4 +20,5 @@ Changes in version 1.7.x Changes in version 1.9.x + Updated plotDMN to work with newest mia version ++ Added plotCCA and plotRDA functions + deprecate agglomerateByRank and agglomerateByPrevalence diff --git a/R/plotAbundance.R b/R/plotAbundance.R index 22bdd805..2d485522 100644 --- a/R/plotAbundance.R +++ b/R/plotAbundance.R @@ -110,7 +110,7 @@ #' ## Compositional barplot with top 5 taxa and samples sorted by "Bacteroidetes" #' #' # Getting top taxa on a Phylum level -#' se <- transformCounts(se, method="relabundance") +#' se <- transformAssay(se, method="relabundance") #' se_phylum <- mergeFeaturesByRank(se, rank ="Phylum", onRankOnly=TRUE) #' top_taxa <- getTopTaxa(se_phylum,top = 5, assay.type = "relabundance") #' diff --git a/R/plotAbundanceDensity.R b/R/plotAbundanceDensity.R index b8d3ea8f..b8d314d6 100644 --- a/R/plotAbundanceDensity.R +++ b/R/plotAbundanceDensity.R @@ -94,7 +94,7 @@ #' plotAbundanceDensity(tse, assay.type = "counts") #' #' # Counts relative abundances -#' tse <- transformCounts(tse, method = "relabundance") +#' tse <- transformAssay(tse, method = "relabundance") #' #' # Plots the relative abundance of 10 most abundant taxa. #' # "nationality" information is used to color the points. X-axis is log-scaled. diff --git a/R/plotCCA.R b/R/plotCCA.R new file mode 100644 index 00000000..1ad2eec3 --- /dev/null +++ b/R/plotCCA.R @@ -0,0 +1,582 @@ +#' Plot RDA or CCA object +#' +#' \code{plotRDA} and \code{plotCCA} create an RDA/CCA plot starting from the +#' output of \code{\link[mia:runCCA]{CCA and RDA}} functions, two common methods +#' for supervised ordination of microbiome data. +#' +#' @param object a +#' \code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-constructor]{TreeSummarizedExperiment}} +#' or a matrix of weights. The latter is returned as output from \code{\link[mia:runCCA]{calculateRDA}}. +#' +#' @param dimred A string or integer scalar indicating the reduced dimension to +#' plot. This is the output of \code{\link[mia:runCCA]{runRDA}} and resides in +#' \code{reducedDim(tse, dimred)}. +#' +#' @param add.ellipse One of \code{c(TRUE, FALSE, "fill", "colour", "color")}, +#' indicating whether ellipses should be present, absent, filled or colored. +#' (default: \code{ellipse.fill = TRUE}) +#' +#' @param ellipse.alpha Number between 0 and 1 to adjust the opacity of ellipses. +#' (default: \code{ellipse.alpha = 0.2}) +#' +#' @param ellipse.linewidth Number specifying the size of ellipses. +#' (default: \code{ellipse.linewidth = 0.1}) +#' +#' @param ellipse.linetype Discrete number specifying the style of ellipses. +#' (default: \code{ellipse.linetype = 1}) +#' +#' @param add.vectors TRUE or FALSE, should vectors appear in the plot. +#' (default: \code{add.vectors = TRUE}) +#' +#' @param vec.size Number specifying the size of vectors. +#' (default: \code{vec.size = 0.5}) +#' +#' @param vec.colour String specifying the colour of vectors. +#' (default: \code{vec.color = "black"}) +#' +#' @param vec.color Alias for `vec.colour`. +#' +#' @param vec.linetype Discrete number specifying the style of vector lines. +#' (default: \code{vec.linetype = 1}) +#' +#' @param arrow.size Number specifying the size of arrows. +#' (defaults: \code{arrow.size = 0.25}) +#' +#' @param label.size Number specifying the size of text and labels. +#' (default: \code{label.size = 4}) +#' +#' @param label.colour String specifying the colour of text and labels. +#' (default: \code{label.color = "black"}) +#' +#' @param label.color Alias for `label.colour`. +#' +#' @param sep.group String specifying the separator used in the labels. +#' (default: \code{sep.group = "\U2012"}) +#' +#' @param repl.underscore String used to replace underscores in the labels. +#' (default: \code{repl.underscore = " "}) +#' +#' @param vec.text TRUE or FALSE, should text instead of labels be used to label vectors. +#' (default: \code{vec.text = TRUE}) +#' +#' @param repel.labels TRUE or FALSE, should labels be repelled. +#' (default: \code{repel.labels = TRUE}) +#' +#' @param parse.labels TRUE or FALSE, should labels be parsed. +#' (default: \code{parse.labels = TRUE}) +#' +#' @param add.significance TRUE or FALSE, should explained variance and p-value +#' appear in the labels. (default: \code{add.significance = TRUE}) +#' +#' @param add.expl.var TRUE or FALSE, should explained variance appear on the +#' coordinate axes. (default: \code{add.expl.var = TRUE}) +#' +#' @param ... additional parameters for plotting, inherited from +#' \code{\link[scater:plotReducedDim]{plotReducedDim}}, +#' \code{\link[ggplot2:geom_label]{geom_label}} and +#' \code{\link[ggrepel:geom_label_repel]{geom_label_repel}}. +#' +#' @details +#' \code{plotRDA} and \code{plotCCA} create an RDA/CCA plot starting from the +#' output of \code{\link[mia:runCCA]{CCA and RDA}} functions, two common methods +#' for supervised ordination of microbiome data. Either a +#' \code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-constructor]{TreeSummarizedExperiment}} +#' or a matrix object is supported as input. When the input is a +#' TreeSummarizedExperiment, this should contain the output of runRDA +#' in the reducedDim slot and the argument \code{dimred} needs to be defined. +#' When the input is a matrix, this should be returned as output from +#' calculateRDA. However, the first method is recommended because it provides +#' the option to adjust aesthetics to the colData variables through the +#' arguments inherited from \code{\link[scater:plotReducedDim]{plotReducedDim}}. +#' +#' @return +#' A \code{ggplot2} object +#' +#' @name plotCCA +#' +#' @examples +#' # Load dataset +#' library(miaViz) +#' data("enterotype", package = "mia") +#' tse <- enterotype +#' +#' # Run RDA and store results into TreeSE +#' tse <- runRDA(tse, +#' formula = assay ~ ClinicalStatus + Gender + Age, +#' FUN = vegan::vegdist, +#' distance = "bray", +#' na.action = na.exclude) +#' +#' # Create RDA plot coloured by variable +#' plotRDA(tse, "RDA", +#' colour_by = "ClinicalStatus") +#' +#' # Create RDA plot with empty ellipses +#' plotRDA(tse, "RDA", +#' colour_by = "ClinicalStatus", +#' add.ellipse = "colour") +#' +#' # Create RDA plot with text encased in labels +#' plotRDA(tse, "RDA", +#' colour_by = "ClinicalStatus", +#' vec.text = FALSE) +#' +#' # Create RDA plot without repelling text +#' plotRDA(tse, "RDA", +#' colour_by = "ClinicalStatus", +#' repel.labels = FALSE) +#' +#' # Create RDA plot without vectors +#' plotRDA(tse, "RDA", +#' colour_by = "ClinicalStatus", +#' add.vectors = FALSE) +#' +#' # Calculate RDA as a separate object +#' rda_mat <- calculateRDA(tse, +#' formula = assay ~ ClinicalStatus + Gender + Age, +#' FUN = vegan::vegdist, +#' distance = "bray", +#' na.action = na.exclude) +#' +#' # Create RDA plot from RDA matrix +#' plotRDA(rda_mat) +NULL + +#' @rdname plotCCA +#' @aliases plotRDA +#' @export +setGeneric("plotCCA", signature = c("object"), + function(object, ...) standardGeneric("plotCCA")) + +#' @rdname plotCCA +#' @aliases plotRDA +#' @export +setMethod("plotCCA", signature = c(object = "SingleCellExperiment"), + function(object, dimred, ...){ + # Reproduce plotRDA function + return(plotRDA(object, dimred, ...)) + } +) + +#' @rdname plotCCA +#' @aliases plotRDA +#' @export +setMethod("plotCCA", signature = c(object = "matrix"), + function(object, ...){ + # Reproduce plotRDA function + return(plotRDA(object, ...)) + } +) + +#' @rdname plotCCA +#' @aliases plotCCA +#' @export +setGeneric("plotRDA", signature = c("object"), + function(object, ...) standardGeneric("plotRDA")) + +#' @rdname plotCCA +#' @aliases plotCCA +#' @export +setMethod("plotRDA", signature = c(object = "SingleCellExperiment"), + function(object, dimred, + add.ellipse = TRUE, ellipse.alpha = 0.2, ellipse.linewidth = 0.1, ellipse.linetype = 1, + vec.size = 0.5, vec.color = vec.colour, vec.colour = "black", vec.linetype = 1, + arrow.size = 0.25, label.color = label.colour, label.colour = "black", label.size = 4, + vec.text = TRUE, repel.labels = TRUE, sep.group = "\U2012", repl.underscore = " ", + add.significance = TRUE, add.expl.var = TRUE, add.vectors = TRUE, parse.labels = TRUE, ...){ + ###################### Input check ######################## + if( !(add.ellipse %in% c(TRUE, FALSE, "fill", "color", "colour")) ){ + stop("'add.ellipse' must be one of c(TRUE, FALSE, 'fill', 'color', 'colour').", call. = FALSE) + } + if ( !.is_a_bool(add.vectors) ){ + stop("'add.vectors must be TRUE or FALSE.", call. = FALSE) + } + if ( !add.vectors ){ + warning("'add.vectors' is FALSE, so other arguments for vectors and labels will be disregarded.", call. = FALSE) + } + if( !.is_a_bool(vec.text) ){ + stop("'vec.text' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(repel.labels) ){ + stop("'repel.labels' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(parse.labels) ){ + stop("'parse.labels' must be TRUE or FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.significance) ){ + stop("'add.significance' must be TRUE or FALSE.", call. = FALSE) + } + if( parse.labels && !add.significance ){ + parse.labels <- FALSE + warning("'parse.labels' was turned off because 'add.significance' is FALSE.", call. = FALSE) + } + if( !.is_a_bool(add.expl.var) ){ + stop("'add.expl.var' must be TRUE or FALSE.", call. = FALSE) + } + if( !.are_whole_numbers(ellipse.linetype) ){ + stop("'ellipse.linetype' must be a whole number.", call. = FALSE) + } + if( !.are_whole_numbers(ellipse.linetype) ){ + stop("'vec.linetype' must be a whole number.", call. = FALSE) + } + if( ellipse.alpha < 0 || ellipse.alpha > 1 ){ + stop("'ellipse.alpha' must be a number between 0 and 1.", call. = FALSE) + } + if ( !is.numeric(ellipse.linewidth) && ellipse.linewidth > 0 ) { + stop("'ellipse.linewidth' must be a positive number.", call. = FALSE) + } + if ( !is.numeric(vec.size) && vec.size > 0 ) { + stop("'vec.size' must be a positive number.", call. = FALSE) + } + if ( !is.numeric(arrow.size) && arrow.size > 0 ) { + stop("'arrow.size' must be a positive number.", call. = FALSE) + } + if ( !is.numeric(label.size) && label.size > 0 ) { + stop("'label.size' must be a positive number.", call. = FALSE) + } + if ( !.is_non_empty_string(vec.color) ) { + stop("'vec.color' must be a non-empty string specifying a colour", call. = FALSE) + } + if ( !.is_non_empty_string(label.color) ) { + stop("'label.color' must be a non-empty string specifying a colour", call. = FALSE) + } + if ( !.is_a_string(sep.group) ) { + stop("'sep.group' must be a string specifying a separator.", call. = FALSE) + } + if ( !.is_a_string(repl.underscore) ) { + stop("'repl.underscore' must be a string.", call. = FALSE) + } + ###################### Input check end #################### + # Get data for plotting + plot_data <- .incorporate_rda_vis( + object, dimred, sep.group = sep.group, repl.underscore = repl.underscore, + add.significance = add.significance, add.expl.var = add.expl.var, + add.ellipse = add.ellipse, add.vectors = add.vectors, ... + ) + # Create a plot + plot <- .rda_plotter( + plot_data, ellipse.alpha = ellipse.alpha, ellipse.linewidth = ellipse.linewidth, + ellipse.linetype = ellipse.linetype, vec.size = vec.size, vec.color = vec.color, + vec.colour = vec.colour, vec.linetype = vec.linetype, arrow.size = arrow.size, + label.color = label.color, label.colour = label.colour, label.size = label.size, + vec.text = vec.text, add.ellipse = add.ellipse, repel.labels = repel.labels, + parse.labels = parse.labels, ... + ) + return(plot) + } +) + +#' @rdname plotCCA +#' @aliases plotCCA +#' @export +setMethod("plotRDA", signature = c(object = "matrix"), + function(object, ...){ + # Construct TreeSE from rda/cca object + object <- .rda2tse(object) + # Run plotRDA method for TreeSE + return(plotRDA(object, "RDA", ...)) + } +) + + +################## HELP FUNCTIONS ########################## + +# Construct TreeSE from rda/cca object to pass it to downstream functions +.rda2tse <- function(object) { + # Convert rda/cca object to TreeSE + object <- TreeSummarizedExperiment( + assays = matrix(ncol = nrow(object), dimnames = list(NULL, rownames(object))), + reducedDims = list(RDA = object) + ) + return(object) +} + +# Get data for plotting +#' @importFrom scater plotReducedDim retrieveCellInfo +#' @importFrom SingleCellExperiment reducedDim reducedDimNames +.incorporate_rda_vis <- function( + tse, dimred, ncomponents = 2, colour_by = color_by, color_by = NULL, + add.significance = TRUE, add.expl.var = TRUE, add.ellipse = TRUE, + add.vectors = TRUE, vec.lab = NULL, expl_var = NULL, + sep.group = "\U2012", repl.underscore = " ", ...){ + + # Check dimred + if( !dimred %in% reducedDimNames(tse) ){ + stop("'dimred' must specify reducedDim.", call. = FALSE) + } + # Get reducedDim + reduced_dim <- reducedDim(tse, dimred) + # Check that there are at least 2 coordinates + if( ncol(reduced_dim) < 2 ){ + stop("reducedDim specified by 'dimred' must have at least 2 columns.", call. = FALSE) + } + # Only 2 dimensions are supported currently + ncomponents <- 2 + + # If specified, get explained variance + if( add.expl.var ){ + # Check if data is available + ind <- names(attributes(reduced_dim)) %in% c("rda", "cca") + # If it can be found + if( any(ind) ){ + # Add explained variance + rda <- attributes(reduced_dim)[ind][[1]] + expl_var <- summary(rda)$concont$importance[2, ]*100 + } else{ + # If it cannot be found, give warning + warning(paste("RDA/CCA object was not found. Please compute", + "RDA/CCA by using runCCA or calculateCCA."), + call. = FALSE) + } + } + + # Get scatter plot with plotReducedDim --> keep theme similar between ordination methods + plot <- plotReducedDim( + tse, dimred = dimred, ncomponents = ncomponents, colour_by = colour_by, + percentVar = expl_var, ...) + + # Get data for ellipse + ellipse_data <- NULL + if( add.ellipse != FALSE && !is.null(colour_by) ){ + ellipse_data <- reduced_dim + ellipse_data <- as.data.frame(ellipse_data) + ellipse_data[[colour_by]] <- retrieveCellInfo(tse, colour_by)[["value"]] + attributes(ellipse_data)$colour_by <- colour_by + } + + # Get data for vectors + vector_data <- NULL + if( add.vectors ){ + # Check if data is available + ind <- names(attributes(reduced_dim)) %in% c("rda", "cca") + # If it can be found + if( any(ind) ){ + # Get biplot from cca object + rda <- attributes(reduced_dim)[ind][[1]] + vector_data <- rda$CCA$biplot + vector_data <- as.data.frame(vector_data) + vector_data[["group"]] <- rownames(vector_data) + } else{ + # If it cannot be found, give warning + warning(paste("RDA/CCA object was not found. Please compute RDA/CCA", + "by using runCCA or calculateCCA."), + call. = FALSE) + } + } + + # Get variable names from sample metadata + coldata <- colData(tse) + variable_names <- colnames(coldata) + # Check if variable names can be found metadata + all_var_found <- FALSE + if( !is.null(vector_data) && length(variable_names) > 0 ){ + all_var_found <- all(colSums( + sapply(rownames(vector_data), function(x) + sapply(variable_names, function(y) grepl(y, x)) )) == 1) + } + + # Get vector labels + if( !is.null(vector_data) ){ + # Get those names that are present in data + + # If user has not provided vector labels + if( is.null(vec.lab) ){ + vector_label <- rownames(vector_data) + # Make labels more tidy + if( all_var_found ){ + vector_label <- .tidy_vector_labels( + vector_label, coldata, + sep.group = sep.group, + repl.underscore = repl.underscore + ) + } + # Add to df + vector_data$vector_label <- vector_label + } else{ + # Check that user-provided labels are correct length + if( length(vec.lab) != nrow(vector_data) ){ + stop("Number of labels in 'vec_lab' do not match with number of vectors.", + call. = FALSE) + } + # If they are, add labels to data + vector_data$vector_label <- vec.lab + } + } + + # Get significance data + signif_data <- NULL + if( add.significance && !is.null(vector_data) && all_var_found ){ + # Check if data is available + ind <- names(attributes(reduced_dim)) %in% c("significance") + # If it can be found + if( any(ind) ){ + # Get biplot from cca object + signif_data <- attributes(reduced_dim)[ind][[1]] + signif_data <- signif_data[["permanova"]] + signif_data <- as.data.frame(signif_data) + # Add significance to vector labels + # Get vector labels + vector_label <- vector_data[["vector_label"]] + # Add significance to vector labels + vector_label <- .add_signif_to_vector_labels( + vector_label, variable_names, signif_data, repl.underscore) + vector_data[["vector_label"]] <- vector_label + } else{ + # If it cannot be found, give warning + warning(paste("Significance data was not found. please compute", + "CCA/RDA by using runCCA or calculateCCA."), + call. = FALSE) + } + } + + # Create a list to return + result <- list( + plot = plot, + ellipse_data = ellipse_data, + vector_data = vector_data + ) + return(result) +} + +# Make vector labels more tidy, i.e, separate variable and group names. +# Replace also underscores with space +.tidy_vector_labels <- function( + vector_label, coldata, sep.group = "\U2012", repl.underscore = " ", ...){ + # Get variable names from sample metadata + var_names <- colnames(coldata) + # Loop through vector labels + vector_label <- sapply(vector_label, FUN = function(name){ + # Get the real variable name from sample metadata + var_name <- var_names[ sapply(var_names, function(x) grepl(x, name)) ] + # If the vector label includes also group name + if( !name %in% var_names ){ + # Get the group name + group_name <- unique( coldata[[var_name]] )[ + which( paste0(var_name, unique( coldata[[var_name]] )) == name ) ] + # Modify vector so that group is separated from variable name + new_name <- paste0(var_name, " ", sep.group, " ", group_name) + } else{ + new_name <- name + } + # Replace underscores with space + new_name <- gsub("_", repl.underscore, new_name) + return(new_name) + }) + return(vector_label) +} + +# This function adds significance info to vector labels +.add_signif_to_vector_labels <- function( + vector_label, var_names, signif_data, repl.underscore = " ", ...){ + # Replace underscores from significance data and variable names to match labels + rownames(signif_data) <- sapply(rownames(signif_data), function(x) gsub("_", repl.underscore, x)) + var_names <- sapply(var_names, function(x) gsub("_", repl.underscore, x)) + # Loop through vector labels + vector_label <- sapply(vector_label, FUN = function(name){ + # Get the real variable name from sample metadata + var_name <- var_names[ sapply(var_names, function(x) grepl(x, name)) ] + # Add percentage how much this variable explains, and p-value + new_name <- expr( + paste(!!name, " (", + !!format( + round( signif_data[var_name, "Explained variance"]*100, 1), + nsmall = 1), "%, ", italic("P"), " = ", + !!gsub("0\\.","\\.", format( + round( signif_data[var_name, "Pr(>F)"], 3), + nsmall = 3)), ")")) + + return(new_name) + }) + return(vector_label) +} + +# Plot based on the data +#' @importFrom ggrepel geom_text_repel geom_label_repel +.rda_plotter <- function( + plot_data, ellipse.alpha = 0.2, ellipse.linewidth = 0.1, ellipse.linetype = 1, + vec.size = 0.5, vec.color = vec.colour, vec.colour = "black", + vec.linetype = 1, arrow.size = 0.25, min.segment.length = 5, + label.color = label.colour, label.colour = "black", label.size = 4, + parse.labels = TRUE, vec.text = TRUE, repel.labels = TRUE, add.ellipse = TRUE, + position = NULL, nudge_x = NULL, nudge_y = NULL, direction = "both", + max.overlaps = 10, check_overlap = FALSE, ...){ + + # Get the scatter plot + plot <- plot_data[["plot"]] + # Add ellipse + if( !is.null(plot_data$ellipse_data) ){ + # Get data and variabe names + data <- plot_data$ellipse_data + xvar <- colnames(data)[[1]] + yvar <- colnames(data)[[2]] + colour_var <- attributes(plot_data$ellipse_data)$colour_by + # Add ellipses to plot (fill or colour the edge) + if( add.ellipse %in% c(TRUE, "fill") ){ + plot <- plot + + stat_ellipse(data = data, + aes(x = .data[[xvar]], y = .data[[yvar]], + color = .data[[colour_var]], fill = after_scale(color)), + geom = "polygon", alpha = ellipse.alpha, + size = ellipse.linewidth, linetype = ellipse.linetype) + } else if ( add.ellipse %in% c("color", "colour") ){ + plot <- plot + + stat_ellipse(data = data, + aes(x = .data[[xvar]], y = .data[[yvar]], color = .data[[colour_var]]), + geom = "polygon", alpha = 0, linetype = ellipse.linetype) + } + + } + # Add vectors + if( !is.null(plot_data$vector_data) ){ + # Get data and variabe names + data <- plot_data$vector_data + xvar <- colnames(data)[[1]] + yvar <- colnames(data)[[2]] + # Add vectors + plot <- plot + + geom_segment(data = data, + aes(x = 0, y = 0, xend = .data[[xvar]], yend = .data[[yvar]], + group = .data[["group"]]), + arrow = arrow(length = unit(arrow.size, "cm")), + color = vec.color, linetype = vec.linetype, size = vec.size) + # Add vector labels (text or label) + # Make list of arguments for geom_text/geom_label + label_args <- list( + data = data, + mapping = aes(x = .data[[xvar]], y = .data[[yvar]]), + label = data[["vector_label"]], parse = parse.labels, + color = label.color, size = label.size, stat = "identity", + nudge_x = nudge_x, nudge_y = nudge_y, show.legend = NA, + na.rm = FALSE, inherit.aes = TRUE + ) + # Repel text + if( repel.labels ){ + # Add arguments for geom_text_repel/geom_label_repel to list + label_args <- c( + label_args, min.segment.length = min.segment.length, + box.padding = 0.25, point.padding = 1e-06, force = 1, force_pull = 1, + max.time = 0.5, max.iter = 10000, max.overlaps = max.overlaps, + direction = direction, seed = NA, verbose = FALSE + ) + # repelled text + if( vec.text ){ + plot <- plot + do.call("geom_text_repel", label_args) + # repelled labels + } else{ + plot <- plot + do.call("geom_label_repel", label_args) + } + # Do not repel text + } else{ + # not repelled text + if( vec.text ){ + label_args <- c(label_args, check_overlap = check_overlap) + plot <- plot + do.call("geom_text", label_args) + # not repelled labels + } else{ + plot <- plot + do.call("geom_label", label_args) + } + } + + } + return(plot) +} diff --git a/R/plotSeries.R b/R/plotSeries.R index d61b23d8..2c8adcfd 100644 --- a/R/plotSeries.R +++ b/R/plotSeries.R @@ -67,7 +67,7 @@ #' colour_by = "Family") #' #' # Counts relative abundances -#' object <- transformCounts(object, method = "relabundance") +#' object <- transformAssay(object, method = "relabundance") #' #' # Selects taxa #' taxa <- c("seq_1", "seq_2", "seq_3", "seq_4", "seq_5") diff --git a/man/plotAbundance.Rd b/man/plotAbundance.Rd index 595d87b6..e24a0a1c 100644 --- a/man/plotAbundance.Rd +++ b/man/plotAbundance.Rd @@ -134,7 +134,7 @@ wrap_plots(plot, ncol = 1, heights = c(0.8,0.2)) ## Compositional barplot with top 5 taxa and samples sorted by "Bacteroidetes" # Getting top taxa on a Phylum level -se <- transformCounts(se, method="relabundance") +se <- transformAssay(se, method="relabundance") se_phylum <- mergeFeaturesByRank(se, rank ="Phylum", onRankOnly=TRUE) top_taxa <- getTopTaxa(se_phylum,top = 5, assay.type = "relabundance") diff --git a/man/plotAbundanceDensity.Rd b/man/plotAbundanceDensity.Rd index 27acde66..8b8a7f3b 100644 --- a/man/plotAbundanceDensity.Rd +++ b/man/plotAbundanceDensity.Rd @@ -112,7 +112,7 @@ tse <- microbiomeDataSets::atlas1006() plotAbundanceDensity(tse, assay.type = "counts") # Counts relative abundances -tse <- transformCounts(tse, method = "relabundance") +tse <- transformAssay(tse, method = "relabundance") # Plots the relative abundance of 10 most abundant taxa. # "nationality" information is used to color the points. X-axis is log-scaled. diff --git a/man/plotCCA.Rd b/man/plotCCA.Rd new file mode 100644 index 00000000..714ca498 --- /dev/null +++ b/man/plotCCA.Rd @@ -0,0 +1,188 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotCCA.R +\name{plotCCA} +\alias{plotCCA} +\alias{plotRDA} +\alias{plotCCA,SingleCellExperiment-method} +\alias{plotCCA,matrix-method} +\alias{plotRDA,SingleCellExperiment-method} +\alias{plotRDA,matrix-method} +\title{Plot RDA or CCA object} +\usage{ +plotCCA(object, ...) + +\S4method{plotCCA}{SingleCellExperiment}(object, dimred, ...) + +\S4method{plotCCA}{matrix}(object, ...) + +plotRDA(object, ...) + +\S4method{plotRDA}{SingleCellExperiment}( + object, + dimred, + add.ellipse = TRUE, + ellipse.alpha = 0.2, + ellipse.linewidth = 0.1, + ellipse.linetype = 1, + vec.size = 0.5, + vec.color = vec.colour, + vec.colour = "black", + vec.linetype = 1, + arrow.size = 0.25, + label.color = label.colour, + label.colour = "black", + label.size = 4, + vec.text = TRUE, + repel.labels = TRUE, + sep.group = "‒", + repl.underscore = " ", + add.significance = TRUE, + add.expl.var = TRUE, + add.vectors = TRUE, + parse.labels = TRUE, + ... +) + +\S4method{plotRDA}{matrix}(object, ...) +} +\arguments{ +\item{object}{a +\code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-constructor]{TreeSummarizedExperiment}} +or a matrix of weights. The latter is returned as output from \code{\link[mia:runCCA]{calculateRDA}}.} + +\item{...}{additional parameters for plotting, inherited from +\code{\link[scater:plotReducedDim]{plotReducedDim}}, +\code{\link[ggplot2:geom_label]{geom_label}} and +\code{\link[ggrepel:geom_label_repel]{geom_label_repel}}.} + +\item{dimred}{A string or integer scalar indicating the reduced dimension to +plot. This is the output of \code{\link[mia:runCCA]{runRDA}} and resides in +\code{reducedDim(tse, dimred)}.} + +\item{add.ellipse}{One of \code{c(TRUE, FALSE, "fill", "colour", "color")}, +indicating whether ellipses should be present, absent, filled or colored. +(default: \code{ellipse.fill = TRUE})} + +\item{ellipse.alpha}{Number between 0 and 1 to adjust the opacity of ellipses. +(default: \code{ellipse.alpha = 0.2})} + +\item{ellipse.linewidth}{Number specifying the size of ellipses. +(default: \code{ellipse.linewidth = 0.1})} + +\item{ellipse.linetype}{Discrete number specifying the style of ellipses. +(default: \code{ellipse.linetype = 1})} + +\item{vec.size}{Number specifying the size of vectors. +(default: \code{vec.size = 0.5})} + +\item{vec.color}{Alias for \code{vec.colour}.} + +\item{vec.colour}{String specifying the colour of vectors. +(default: \code{vec.color = "black"})} + +\item{vec.linetype}{Discrete number specifying the style of vector lines. +(default: \code{vec.linetype = 1})} + +\item{arrow.size}{Number specifying the size of arrows. +(defaults: \code{arrow.size = 0.25})} + +\item{label.color}{Alias for \code{label.colour}.} + +\item{label.colour}{String specifying the colour of text and labels. +(default: \code{label.color = "black"})} + +\item{label.size}{Number specifying the size of text and labels. +(default: \code{label.size = 4})} + +\item{vec.text}{TRUE or FALSE, should text instead of labels be used to label vectors. +(default: \code{vec.text = TRUE})} + +\item{repel.labels}{TRUE or FALSE, should labels be repelled. +(default: \code{repel.labels = TRUE})} + +\item{sep.group}{String specifying the separator used in the labels. +(default: \code{sep.group = "\U2012"})} + +\item{repl.underscore}{String used to replace underscores in the labels. +(default: \code{repl.underscore = " "})} + +\item{add.significance}{TRUE or FALSE, should explained variance and p-value +appear in the labels. (default: \code{add.significance = TRUE})} + +\item{add.expl.var}{TRUE or FALSE, should explained variance appear on the +coordinate axes. (default: \code{add.expl.var = TRUE})} + +\item{add.vectors}{TRUE or FALSE, should vectors appear in the plot. +(default: \code{add.vectors = TRUE})} + +\item{parse.labels}{TRUE or FALSE, should labels be parsed. +(default: \code{parse.labels = TRUE})} +} +\value{ +A \code{ggplot2} object +} +\description{ +\code{plotRDA} and \code{plotCCA} create an RDA/CCA plot starting from the +output of \code{\link[mia:runCCA]{CCA and RDA}} functions, two common methods +for supervised ordination of microbiome data. +} +\details{ +\code{plotRDA} and \code{plotCCA} create an RDA/CCA plot starting from the +output of \code{\link[mia:runCCA]{CCA and RDA}} functions, two common methods +for supervised ordination of microbiome data. Either a +\code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-constructor]{TreeSummarizedExperiment}} +or a matrix object is supported as input. When the input is a +TreeSummarizedExperiment, this should contain the output of runRDA +in the reducedDim slot and the argument \code{dimred} needs to be defined. +When the input is a matrix, this should be returned as output from +calculateRDA. However, the first method is recommended because it provides +the option to adjust aesthetics to the colData variables through the +arguments inherited from \code{\link[scater:plotReducedDim]{plotReducedDim}}. +} +\examples{ +# Load dataset +library(miaViz) +data("enterotype", package = "mia") +tse <- enterotype + +# Run RDA and store results into TreeSE +tse <- runRDA(tse, + formula = assay ~ ClinicalStatus + Gender + Age, + FUN = vegan::vegdist, + distance = "bray", + na.action = na.exclude) + +# Create RDA plot coloured by variable +plotRDA(tse, "RDA", + colour_by = "ClinicalStatus") + +# Create RDA plot with empty ellipses +plotRDA(tse, "RDA", + colour_by = "ClinicalStatus", + add.ellipse = "colour") + +# Create RDA plot with text encased in labels +plotRDA(tse, "RDA", + colour_by = "ClinicalStatus", + vec.text = FALSE) + +# Create RDA plot without repelling text +plotRDA(tse, "RDA", + colour_by = "ClinicalStatus", + repel.labels = FALSE) + +# Create RDA plot without vectors +plotRDA(tse, "RDA", + colour_by = "ClinicalStatus", + add.vectors = FALSE) + +# Calculate RDA as a separate object +rda_mat <- calculateRDA(tse, + formula = assay ~ ClinicalStatus + Gender + Age, + FUN = vegan::vegdist, + distance = "bray", + na.action = na.exclude) + +# Create RDA plot from RDA matrix +plotRDA(rda_mat) +} diff --git a/man/plotSeries.Rd b/man/plotSeries.Rd index 75405a2e..e45fa6aa 100644 --- a/man/plotSeries.Rd +++ b/man/plotSeries.Rd @@ -96,7 +96,7 @@ plotSeries(object, colour_by = "Family") # Counts relative abundances -object <- transformCounts(object, method = "relabundance") +object <- transformAssay(object, method = "relabundance") # Selects taxa taxa <- c("seq_1", "seq_2", "seq_3", "seq_4", "seq_5") diff --git a/tests/testthat/test-2plotSeries.R b/tests/testthat/test-2plotSeries.R index 817e3a03..b509ae04 100644 --- a/tests/testthat/test-2plotSeries.R +++ b/tests/testthat/test-2plotSeries.R @@ -4,7 +4,7 @@ test_that("plot series", { # Load data from miaTime package if( !require("miaTime") ){ if( !require("devtools") ){ - install.packages("devtools") + BiocManager::install("devtools") } devtools::install_github("microbiome/miaTime") library("miaTime") diff --git a/tests/testthat/test-plotCCA.R b/tests/testthat/test-plotCCA.R new file mode 100644 index 00000000..2ac58037 --- /dev/null +++ b/tests/testthat/test-plotCCA.R @@ -0,0 +1,80 @@ +test_that("plot RDA/CCA", { + data("Tengeler2020", package = "mia") + tse <- Tengeler2020 + + ### 1). TEST error messages ### + + # Object without reducedDim + expect_error(plotRDA(tse), 'argument "dimred" is missing, with no default') + expect_error(plotRDA(tse, "RDA"), "'dimred' must specify reducedDim.") + + # Run/calculate RDA + tse <- runRDA(tse, assay ~ patient_status + cohort) + rda <- calculateRDA(tse, assay ~ patient_status + cohort) + + # Minimal functionality + expect_no_error(plotRDA(tse, "RDA")) + + # Wrong-entry scenarios + expect_error(plotRDA(tse, "RDA", colour_by = "wrong colname")) + expect_error(plotRDA(tse, "RDA", colour_by = "cohort", shape_by = "wrong colname")) + expect_error(plotRDA(tse, "RDA", add.ellipse = "invalid value"), + "'add.ellipse' must be one of c(TRUE, FALSE, 'fill', 'color', 'colour').", + fixed = TRUE) + expect_error(plotRDA(tse, "RDA", add.significance = "invalid value"), + "'add.significance' must be TRUE or FALSE.") + expect_error(plotRDA(tse, "RDA", add.expl.var = "invalid value"), + "'add.expl.var' must be TRUE or FALSE.") + expect_error(plotRDA(tse, "RDA", repel.labels)) + expect_warning(plotRDA(tse, "RDA", add.significance = FALSE, parse.labels = TRUE), + "'parse.labels' was turned off because 'add.significance' is FALSE.") + expect_warning(plotRDA(tse, "RDA", add.vectors = FALSE), + "'add.vectors' is FALSE, so other arguments for vectors and labels will be disregarded.") + ## add more tests here + + ### 2). TEST plot layers ### + + el_true <- plotRDA(tse, "RDA", colour_by = "patient_status") + el_false <- plotRDA(tse, "RDA", colour_by = "patient_status", add.ellipse = FALSE) + el_col <- plotRDA(tse, "RDA", colour_by = "patient_status", add.ellipse = "colour") + el_fill <- plotRDA(tse, "RDA", colour_by = "patient_status", add.ellipse = "fill") + expect_warning( + vec_false <- plotRDA(tse, "RDA", colour_by = "patient_status", add.vectors = FALSE) + ) + # Filled ellipse has one more layer than no ellipse plot + expect_equal(length(ggplot_build(el_true)[["data"]]), 4) + expect_equal(length(ggplot_build(el_false)[["data"]]), 3) + # No-vector plot has only 2 layers, (points and ellipse) + expect_equal(length(ggplot_build(vec_false)[["data"]]), 2) + # Coloured ellipse but not filled ellipse plot has all 0 alpha values + expect_true(all(ggplot_build(el_col)[["data"]][[2]][["alpha"]] == 0)) + expect_false(all(ggplot_build(el_fill)[["data"]][[2]][["alpha"]] == 0)) + + # Check ggplot aesthetics + p_aes <- plotRDA(tse, "RDA", colour_by = "patient_status", ellipse.alpha = 0.5, + ellipse.linewidth = 0.2, ellipse.linetype = 3, vec.size = 0.6, + vec.colour = "red", vec.linetype = 2, arrow.size = 0.15, + label.colour = "blue", label.size = 5) + # Build plot and get data + p_aes_build <- ggplot_build(p_aes)[["data"]] + # Ellipse aesthetics are correctly defined in ggplot + expect_true(all(p_aes_build[[2]][["alpha"]] == 0.5)) + expect_true(all(ggplot_build(p_aes)[[2]][["linewidth"]] == 0.2)) + expect_true(all(ggplot_build(p_aes)[[2]][["linetype"]] == 3)) + # Vector aesthetics are correctly defined in ggplot + expect_true(all(p_aes_build[[3]][["colour"]] == "red")) + expect_true(all(p_aes_build[[3]][["linewidth"]] == 0.6)) + expect_true(all(p_aes_build[[3]][["linetype"]] == 2)) + # Label aesthetics are correctly defined in ggplot + expect_true(all(p_aes_build[[4]][["colour"]] == "blue")) + expect_true(all(p_aes_build[[4]][["size"]] == 5)) + # Where is arrow size stored? + # expect_true(arrow_size == 0.15) + + # Vector or label text + p_vec <- plotRDA(tse, "RDA", colour_by = "patient_status", vec.text = TRUE) + p_lab <- plotRDA(tse, "RDA", colour_by = "patient_status", vec.text = FALSE) + # Column "fill" is present in p_vec and missing in p_lab, so length differs by 1 + expect_length(ggplot_build(p_vec)[["data"]][[4]], 29) + expect_length(ggplot_build(p_lab)[["data"]][[4]], 28) +}) diff --git a/vignettes/miaViz.Rmd b/vignettes/miaViz.Rmd index f8951e39..e56f29f7 100644 --- a/vignettes/miaViz.Rmd +++ b/vignettes/miaViz.Rmd @@ -63,7 +63,7 @@ However, if the `rank` is set not `NULL` a bar plot is returned. At the same time the `features` argument can be set to `NULL` (default). ```{r} -GlobalPatterns <- transformCounts(GlobalPatterns, method = "relabundance") +GlobalPatterns <- transformAssay(GlobalPatterns, method = "relabundance") ``` ```{r} @@ -256,12 +256,21 @@ plotColGraph(altExp(GlobalPatterns,"Genus"), # Plotting of serial data -```{r} +```{r, eval=FALSE} +if(!requireNamespace("devtools", quietly = TRUE)){ + BiocManager::install("devtools") +} +if(!requireNamespace("miaTime", quietly = TRUE)){ + devtools::install_github("microbiome/miaTime", upgrade = "never") +} +``` + +```{r, eval=FALSE} # Load data from miaTime package library("miaTime") data(SilvermanAGutData, package="miaTime") tse <- SilvermanAGutData -tse <- transformCounts(tse, method = "relabundance") +tse <- transformAssay(tse, method = "relabundance") taxa <- getTopTaxa(tse, 2) ``` @@ -271,7 +280,7 @@ descriptor for ordering the data. The `y` argument selects the feature to show. Since plotting a lot of features is not advised a maximum of 20 features can plotted at the same time. -```{r} +```{r, eval=FALSE} plotSeries(tse, x = "DAY_ORDER", y = taxa, @@ -281,7 +290,7 @@ If replicated data is present, data is automatically used for calculation of the `mean` and `sd` and plotted as a range. Data from different assays can be used for plotting via the `assay.type`. -```{r} +```{r, eval=FALSE} plotSeries(tse[taxa,], x = "DAY_ORDER", colour_by = "Family", @@ -291,7 +300,7 @@ plotSeries(tse[taxa,], Additional variables can be used to modify line type aesthetics. -```{r} +```{r, eval=FALSE} plotSeries(tse, x = "DAY_ORDER", y = getTopTaxa(tse, 5), @@ -329,23 +338,14 @@ plotDMNFit(dmn_se, type = "laplace") # Serial data ordination and trajectories -```{r, error=FALSE, warning=FALSE, results='hide'} -if(!requireNamespace("devtools", quietly = TRUE)){ - BiocManager::install("devtools") -} -if(!requireNamespace("miaTime", quietly = TRUE)){ - devtools::install_github("microbiome/miaTime") -} -``` - Principal Coordinates Analysis using Bray-Curtis dissimilarity on the `hitchip1006` dataset: -```{r MDS, include=FALSE} +```{r MDS, eval=FALSE} library(miaTime) data(hitchip1006, package = "miaTime") tse <- hitchip1006 -tse <- transformCounts(tse, method = "relabundance") +tse <- transformAssay(tse, method = "relabundance") ## Ordination with PCoA with Bray-Curtis dissimilarity tse <- runMDS(tse, FUN = vegan::vegdist, method = "bray", name = "PCoA_BC", assay.type = "relabundance", na.rm = TRUE) @@ -356,7 +356,7 @@ p Retrieving information about all available trajectories: -```{r} +```{r, eval=FALSE} library(dplyr) # List subjects with two time points @@ -369,7 +369,7 @@ table(table(tse$subject)) %>% as.data.frame() %>% Lets look at all trajectories having two time points in the data: -```{r} +```{r, eval=FALSE} # plot p + geom_path(aes(x=X1, y=X2, group=subject), arrow=arrow(length = unit(0.1, "inches")), @@ -383,7 +383,7 @@ p + geom_path(aes(x=X1, y=X2, group=subject), Filtering the two time point trajectories by divergence and displaying top 10%: -```{r} +```{r, eval=FALSE} library(miaTime) # calculating step wise divergence based on the microbial profiles tse <- getStepwiseDivergence(tse, group = "subject", time_field = "time") @@ -411,7 +411,7 @@ p + geom_path(aes(x=X1, y=X2, Plotting an example of the trajectory with the maximum total divergence: -```{r} +```{r, eval=FALSE} # Get subject with the maximum total divergence selected.subject <- data.frame(reducedDim(tse), colData(tse)) %>% group_by(subject) %>%