fantaxtic

Fantaxtic - Nested Bar Plots for Phyloseq Data

https://github.com/gmteunisse/fantaxtic

Keywords

Repository

Fantaxtic - Nested Bar Plots for Phyloseq Data

README.Rmd

---
output: github_document
---



```{r, echo = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-"
)
```

# fantaxtic
`fantaxtic` contains a set of functions to identify and visualize the most abundant taxa in phyloseq objects. It allows users to identify top taxa using any metric and any grouping, and plot the (relative) abundances of the top taxa using a nested bar plot visualisation. In the nested bar plot, colours or fills signify a top taxonomic rank (e.g. Phylum), and a gradient of shades and tints signifies levels at a nested taxonomic rank (e.g. Species). It is particularly useful to present an overview of microbiome sequencing, amplicon sequencing or metabarcoding data.

Note that `fantaxtic` is essentially a wrapper around [`ggnested`](https://github.com/gmteunisse/ggnested), with some accessory functions to identify top taxa and to ensure that the plot is useful. Thus, the output is `ggplot2` object, and can be manipulated as such.


**Keywords: nested bar plot, phyloseq, taxonomy, most abundant taxa, multiple levels, shades, tints, gradient, 16S, ITS ,18S, microbiome, amplicon sequencing, metabarcoding**

## Installation
```{r, eval = FALSE}
if(!"devtools" %in% installed.packages()){
  install.packages("devtools")
}
devtools::install_github("gmteunisse/fantaxtic")
```

## Basic usage
The \code{fantaxtic} workflow consists of two parts:

1. Identify top taxa using either `top_taxa` or `nested_top_taxa`
2. Visualise the top taxa using `nested_bar_plot`

For basic usage, only a few lines of R code are required. To identify and plot the top 10 most abundant ASVs by their mean relative abundance, using Phylum as the top rank and Species as the nested rank, run:

```{r, results = 'hide', message=F, warning=F, error=F}
require("fantaxtic")
require("phyloseq")
require("tidyverse")
require("magrittr")
require("ggnested")
require("knitr")
require("gridExtra")
```

```{r, fig.width = 6, fig.height = 6}
data(GlobalPatterns)
top_asv <- top_taxa(GlobalPatterns, n_taxa = 10)
plot_nested_bar(ps_obj = top_asv$ps_obj,
                top_level = "Phylum",
                nested_level = "Species")

```

To identify and plot the top 3 most abundant Phyla, and the top 3 most abundant species within those Phyla, run:

```{r, fig.width = 6, fig.height = 6}
top_nested <- nested_top_taxa(GlobalPatterns,
                              top_tax_level = "Phylum",
                              nested_tax_level = "Species",
                              n_top_taxa = 3, 
                              n_nested_taxa = 3)
plot_nested_bar(ps_obj = top_nested$ps_obj,
                top_level = "Phylum",
                nested_level = "Species")

```

### `top_taxa`
This function identifies the top $n$ taxa by some metric (e.g. mean, median, variance, etc.) in a phyloseq object. It outputs a table with the top taxa, as well as a phyloseq object in which all other taxa have been merged into a single taxon. 

#### Taxonomic rank
By default, `top_taxa` runs the analysis at the ASV level; however, if a `tax_level` is specified (e.g. `Species`), it first agglomerates the taxa in the phyloseq object at that rank and then runs the analysis. Note that taxonomic agglomeration makes the assumption that taxa with the same name at all ranks are identical. This also includes taxa with missing annotations (`NA`). By default, `top_taxa` does not considered taxa with an `NA` annotation at `tax_level`, but this can be overcome by setting `include_na_taxa = T`.

```{r}
top_species <- top_taxa(GlobalPatterns,
                        n_taxa = 10, 
                        tax_level = "Species")
top_species$top_taxa %>%
  mutate(abundance = round(abundance, 3)) %>%
  kable(format = "markdown")

```

#### Grouping
Furthermore, if one or more grouping factors are specified in `grouping`, it will calculate the top $n$ taxa using the samples in each group, rather than using all samples in the phyloseq object. This makes it possible to for example identify the top taxa in each sample, or the top taxa in each treatment group.

```{r}
top_grouped <- top_taxa(GlobalPatterns,
                        n_taxa = 1,
                        grouping = "SampleType")
top_grouped$top_taxa %>%
  mutate(abundance = round(abundance, 3)) %>%
  kable(format = "markdown")

```

#### Ranking metric
Lastly, any metric can be used to rank taxa by specifying a function through `FUN`. The mean is used by default, but depending on your analysis, you might want to use the median, variance, maximum or any other function that takes as input a numeric vector and outputs a single number.

```{r}
top_max <- top_taxa(GlobalPatterns,
                        n_taxa = 10,
                        FUN = max)
top_max$top_taxa %>%
  mutate(abundance = round(abundance, 3)) %>%
  kable(format = "markdown")

```

### `nested_top_taxa`
This function identifies the top $n$ taxa at a taxonomic rank (e.g. Phylum) and the top $m$ nested taxa at a lower taxonomic rank (e.g. Species) by some metric (e.g. mean, median, variance, etc.) in a phyloseq object. Internally, it makes use of `top_taxa`, and therefore uses many of the same options. Like `top_taxa`, it agglomerates taxa at the specified taxonomic ranks before identifying the top taxa. It outputs a table with the top taxa, as well as a phyloseq object in which all non-top taxa have been merged, both at the `top_tax_level` and at the `nested_tax_level`. This function is especially nice for providing overviews of your data, as it shows the relative abundance of each select `top_tax_level` taxon.

```{r}
top_nested <- nested_top_taxa(GlobalPatterns,
                              top_tax_level = "Phylum",
                              nested_tax_level = "Species",
                              n_top_taxa = 3, 
                              n_nested_taxa = 3, 
                              nested_merged_label = "NA and other ")
top_nested$top_taxa %>%
  mutate(top_abundance = round(top_abundance, 3),
         nested_abundance = round(nested_abundance, 3)) %>%
  kable(format = "markdown")

```

### `plot_nested_bar`
This function is analogous in use to the phyloseq function `plot_bar`, but plots the abundances of taxa in each sample at two levels: a top level (e.g. Phylum) using colours, and a nested level (e.g. Species) using shades and tints of each colour. It is intended to be used in conjunction with `top_taxa` or `nested_top_taxa`, but works with any phyloseq object. The output is a `ggplot2` object that is generated by `ggnested`, which means it can be further customized by the user, for example with faceting, themes, labels et cetera. [Also see the documentation for `ggnested`.](https://github.com/gmteunisse/ggnested)

```{r, fig.width = 6, fig.height = 6}
plot_nested_bar(top_nested$ps_obj,
                top_level = "Phylum",
                nested_level = "Species",
                nested_merged_label = "NA and other ",
                legend_title = "Phylum and species") +
  facet_wrap(~SampleType,
             scales = "free_x") +
  labs(title = "Relative abundances of the top 3 species for each of the top 3 phyla") +
  theme(plot.title = element_text(hjust = 0.5, 
                                  size = 8, 
                                  face = "bold"),
        legend.key.size = unit(10, 
                               "points"))


```

#### Palette
`plot_nested_bar` automatically generates as many colours as there are `top_level` taxa, starting from a `base_clr`. Furthermore, it colours the merged non-top taxa a different `merged_clr`, by default grey. Sometimes it is necessary to alter the colours of the plot, which is why custom `top_level` palettes can be provided. If the palette is named, the colours will be assigned appropriately. Note that it is not necessary to provide a complete palette; any missing colours will get generated automatically. This can be useful if you are trying to match colours between different plots.

```{r, fig.width = 6, fig.height = 6}
plot_nested_bar(top_nested$ps_obj,
                top_level = "Phylum",
                nested_level = "Species",
                nested_merged_label = "NA and other",
                palette = c(Bacteroidetes = "red", 
                            Proteobacteria = "blue"),
                merged_clr = "black")


```

#### Sample order
By default, `plot_nested_bar` orders samples alphabetically. However, sometimes it may be insightful to alter the sample ordering, for example based on grouping or the abundance of a specific taxon. This can be achieved by supplying a character vector with the sample names in the desired order to `sample_order`.

```{r, fig.width = 6, fig.height = 6}
# Order samples by the total abundance of Proteobacteria
sample_order <- psmelt(top_nested$ps_obj) %>%
  data.frame() %>%
  
  # Calculate relative abundances
  group_by(Sample) %>%
  mutate(Abundance = Abundance / sum(Abundance)) %>%
  
  # Sort by taxon of interest
  filter(Phylum == "Proteobacteria") %>%
  group_by(Sample) %>%
  summarise(Abundance = sum(Abundance)) %>%
  arrange(Abundance) %>% 
  
  # Extract the sample order
  pull(Sample) %>%
  as.character()

# Plot
plot_nested_bar(top_nested$ps_obj,
                "Phylum",
                "Species",
                sample_order = sample_order,
                nested_merged_label = "NA and other")

```

## Advanced usage
The `plot_nested_bar` function will suffice for most purposes. However, sometimes, additional control over the plot is required. While some aspects can be controlled by altering additional arguments of `plot_nested_bar`, it may sometimes be necessary to generate the plot from scratch. `plot_nested_bar` is nothing more than a wrapper around the following functions, and can therefore be recreated manually if required:

* Generate a palette using `taxon_colours`
* Generate names for NA taxa using `name_na_taxa`
* Label identical taxa using `label_duplicate_taxa`
* Convert the phyloseq object to a data frame using `psmelt`
* Relevel the merged taxa using `move_label` and `move_nested_labels`
* Reorder the taxa
* Generate a nested barplot using `ggnested`

Thus, for advanced usage, copy the code chunk below and modify it to your requirements.

```{r, fig.width = 6, fig.height = 6}
# Get the top taxa
top_level <- "Phylum"
nested_level <- "Species"
sample_order <- NULL
top_asv <- top_taxa(GlobalPatterns, n_taxa = 10)

# Create names for NA taxa
ps_tmp <- top_asv$ps_obj %>%
  name_na_taxa()

# Add labels to taxa with the same names
ps_tmp <- ps_tmp %>%
  label_duplicate_taxa(tax_level = nested_level)

# Generate a palette basedon the phyloseq object
pal <- taxon_colours(ps_tmp,
                     tax_level = top_level)

# Convert physeq to df
psdf <- psmelt(ps_tmp)

# Move the merged labels to the appropriate positions in the plot:
# Top merged labels need to be at the top of the plot,
# nested merged labels at the bottom of each group
psdf <- move_label(psdf = psdf,
                   col_name = top_level,
                   label = "Other",
                   pos = 0)
psdf <- move_nested_labels(psdf,
                           top_level = top_level,
                           nested_level = nested_level,
                           top_merged_label = "Other",
                           nested_label = "Other",
                           pos = Inf)

# Reorder samples
if(!is.null(sample_order)){
  if(all(sample_order %in% unique(psdf$Sample))){
    psdf <- psdf %>%
      mutate(Sample = factor(Sample, levels = sample_order))
  } else {
    stop("Error: not all(sample_order %in% sample_names(ps_obj)).")
  }

}

# Generate a base plot
p <- ggnested(psdf,
              aes_string(main_group = top_level,
                         sub_group = nested_level,
                         x = "Sample",
                         y = "Abundance"),
              main_palette = pal) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_nested(theme_light) +
  theme(axis.text.x = element_text(hjust = 1, vjust = 0.5, angle = 90))

# Add relative abundances
p <- p + geom_col(position = position_fill())
p

```


### `name_na_taxa`
More often than not, ASVs will not have a complete taxonomic annotation down to the species level. In phyloseq object, this results in `NA` for any unavailable taxonomic rank. This creates issues when trying to plot at a low taxonomic rank such as Species. `name_na_taxa` resolves this problem by assigning the name of the lowest known taxonomic rank for every `NA` value in each ASV. To make it clear that the newly inferred name is not specific to the rank, it includes the rank from which the name was inferred in the new name. This can be turned off by setting `include_rank = F`.

```{r}
# Fill in names for NA taxa, including their rank
ps_tmp <- name_na_taxa(top_asv$ps_obj)
tax_table(ps_tmp) %>%
  kable(format = "markdown")
```

```{r}
# Leave the rank out and alter the label
ps_tmp <- name_na_taxa(top_asv$ps_obj, 
                       include_rank = F, 
                       na_label = "NA ")
tax_table(ps_tmp) %>%
  kable(format = "markdown")
```

### `taxon_colours`
This function generates a colour for each taxon at a specified rank in a phyloseq object. Custom palettes can also be provided, and if they are named, colours wil be assigned appropriately.

```{r, fig.width = 6, fig.height = 3}
# Function to plot the colours in a palette
plot_colours <- function(pal){
  pal %>%
    data.frame(name = names(.),
               colour = .) %>%
    ggplot(aes(x = 1,
               y = name,
               fill = colour,
               label = paste(name, "-", colour))) +
    geom_tile() +
    geom_text() +
    scale_fill_identity() +
    scale_x_discrete(expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0)) +
    theme(axis.text = element_blank(),
          axis.title = element_blank(),
          axis.ticks = element_blank(),
          plot.title = element_text(hjust = 0.5,
                                    size = 12))
}

# Get top taxa and generate a palette
pal <- taxon_colours(top_asv$ps_obj, tax_level = "Phylum")
p1 <- plot_colours(pal) +
  ggtitle("Default")

# Generate a palette with a different base_clr
pal2 <- taxon_colours(top_asv$ps_obj, tax_level = "Phylum", base_clr = "blue")
p2 <- plot_colours(pal2) +
  ggtitle("Base colour blue")

# Provide a custom incomplete palette
pal3 <- taxon_colours(top_asv$ps_obj, 
                     tax_level = "Phylum", 
                     palette = c(Cyanobacteria = "blue", 
                                 Bacteroidetes = "pink"))
p3 <- plot_colours(pal3) +
  ggtitle("Incomplete custom palette")
grid.arrange(p1, p2, p3, nrow = 1)

```

### `label_duplicate_taxa`
Another issue is that a single taxon may be represented by multiple ASVs, especially when low-rank annotations such as Species are missing. For some studies, it may be important to differentiate between these ASVs. Therefore, ASVs with the same taxonomy need to be assigned unique label to differentiate between them. `label_duplicate_taxa` identifies identical taxa and assigns either a count or the ASV name (taken from `row.names(tax_table(ps_obj))`) to these taxa.

```{r}
# Label the lowest non-NA level
ps_tmp <- label_duplicate_taxa(top_asv$ps_obj, 
                               tax_level = "Genus")
tax_table(ps_tmp)  %>%
  kable(format = "markdown")
```

```{r}
# Use ASVs as ids rather than counts
ps_tmp <- label_duplicate_taxa(top_asv$ps_obj,
                               tax_level = "Genus",
                               asv_as_id = T, 
                               duplicate_label = " ASV ")
tax_table(ps_tmp) %>%
  kable(format = "markdown")

```

### `move_label` & `move_nested_labels`
`move_label` and `move_nested_label` reorder the factors in a specified column of a dataframe, so that in `plot_nested_bar` they appear in the right position. This is mainly used to place the `Other` merged taxa at the top of the plot, and the `Other ` nested merged taxa at the bottom of each group. However, it can be used to move any taxon to any position. For example, if you are ordering your plot by the abundance of a specific Phylum, it may be desirable to make that Phylum appear at the bottom of the plot using `pos = Inf`.

```{r}
# Turn physeq object into a dataframe
ps_tmp <- name_na_taxa(top_asv$ps_obj)
ps_tmp <- label_duplicate_taxa(ps_tmp, tax_level = "Species")
psdf <- psmelt(ps_tmp)
levels(as.factor(psdf$Phylum))

# Move the other label to the start, and Bacteroidetes to the end
psdf <- move_label(psdf, col_name = "Phylum", label =  "Other", pos = 0)
psdf <- move_label(psdf, col_name = "Phylum", label =  "Bacteroidetes", pos = Inf)
levels(psdf$Phylum)

```

Likewise, any nested label can be moved to any desired position. `grep` will be used to find any taxon that contains the `nested_label`, after which it will be moved to the desired nested position. Here we flip the order of the two `Unknown Bacteroides`.
```{r}
levels(as.factor(psdf$Species))
psdf <- move_nested_labels(psdf, 
                           top_level = "Phylum", 
                           nested_level = "Species",
                           top_merged_label = "Other",
                           nested_label = "Unknown Bacteroides 1",
                           pos = Inf)
levels(psdf$Species)


```

### `ggnested`
Finally, the melted phyloseq object can be plotted using `ggnested` and `ggplot`. This gives full control over all aspects of the plot, including details such as the bar width. For more details, see the documentation of [`ggnested`](https://github.com/gmteunisse/ggnested).

```{r, fig.width = 6, fig.height = 6}
# Generate a base plot
p <- ggnested(psdf,
              aes_string(main_group = top_level,
                         sub_group = nested_level,
                         x = "Sample",
                         y = "Abundance"),
              main_palette = pal) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_nested(theme_light) +
  theme(axis.text.x = element_text(hjust = 1, vjust = 0.5, angle = 90))

# Add relative abundances
p <- p + geom_col(position = position_fill(), width = 0.5)
p

```

#### Alternative plot types
By using `ggnested` directly rather than using `plot_nested_bar`, any type of plot can be created that is available through `ggplot2`. For example, you could create a boxplot that is grouped and coloured by phylum, and shaded by species:

```{r, fig.width = 6, fig.height = 6}
# Create a boxplot instead of a barplot
psdf_rel <- psdf %>%
  group_by(Sample) %>%
  mutate(Abundance = Abundance / sum(Abundance)) %>%
  ungroup() %>%
  filter(Phylum != "Other")
p <- ggnested(psdf_rel,
              aes_string(main_group = top_level,
                         sub_group = nested_level,
                         x = "Phylum",
                         y = "Abundance",
                         grouping = "Species"),
              main_palette = pal,
              main_keys = T) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_nested(theme_light) +
  theme(axis.text.x = element_text(hjust = 1, vjust = 0.5, angle = 90))
p + geom_boxplot(alpha = 0.5)

```

## What happened to `fantaxtic_bar`?
All the functions in `fantaxtic` that were required to generate a `fantaxtic_bar`, such as `get_top_taxa`, have been deprecated. However, they are still functional and can be called. Many warnings will be issued, but this does not affect their functionality. See `help("fantaxtic-deprecated")` for deprecated functions.

```{r, fig.width = 6, fig.height = 6, warning = F}
ps_tmp <- get_top_taxa(GlobalPatterns, n = 10)
ps_tmp <- name_taxa(ps_tmp, label = "Unkown", species = T, other_label = "Other")
fantaxtic_bar(ps_tmp, color_by = "Phylum", label_by = "Species", other_label = "Other")

```

Owner

Citation (CITATION.cff)

cff-version: 1.2.0
message: "If you use this software, please cite it as below."
authors:
  - family-names: Teunisse
    given-names: Guus Martijn
    orcid: https://orcid.org/0000-0003-3593-4893
title: "Fantaxtic - Nested Bar Plots for Phyloseq Data"
version: 0.2.1
url: https://github.com/gmteunisse/Fantaxtic
date-released: 2022-09-01

GitHub Events

Issues and Pull Requests

Last synced: 4 months ago