Creating a relative abundance bar plot in R with qiime2 data

Other Bioinformatics Tools
1

I am trying to use R to create a relative abundance chart using my qiime2 data. However I have having an issue. I want this plot to represent the genus and species levels cause I have 7 bacteria, however two of them are Wolbachia unidentified bacteria (but they are distinct from each other). The bar chart is grouping them together and I do not want this to happen.

Plot:

Rplot11
Rplot111023×686 79.6 KB

Code:
library(phyloseq)
library(qiime2R)
library(tidyverse)
library(permute)
library(lattice)
library(vegan)
library(picante)
library(cowplot)
library(ggplot2)

setwd("C:/Users/lily2/Desktop/Masters")

phy <- qza_to_phyloseq (features = "table_lily-noqueen.qza", "Aphaenogaster_taxonomy_3.qza",
tree = "rooted-tree.qza")
metadata<-import_qiime_sample_data("metadata_lily_forR.txt")
phy_metadata <- merge_phyloseq(phy, metadata)
phy_metadata

shannon <- read_qza("shannon_vector.qza")
rich <- plot_richness(phy)

phy_rel_abund <- transform_sample_counts(phy, function(x) {x / sum(x)})
plot_bar(phy_rel_abund, fill = "Genus") +
geom_bar(aes(color = Genus, fill = Genus), stat = "identity", position = "stack") +
labs(x = "", y="Relative abundance") +
theme_q2r() +
ggtitle("Relative abundance stack bar plot") +
theme(axis.title = element_text(color="black", face="bold", size=10)) +
theme(plot.title = element_text(color="black", face = "bold", size =12, hjust = 0.5))

2 Likes
2

Hi there,

I think maybe the best way to do this is using RESCRIPt before importing the taxonomy QZA to phyloseq. See this RESCRIPt tutorial section explaining how to do it. I think this fits your needs, as the tutorial section says:

If for any reason you still want to do this directly in R, without getting too much into it (because I am not very familiar with phyloseq), AFAIK the phylogeny slot of the Phyloseq object is a character matrix, so some rudimentary base R code like:

mat[mat[, "ID_column"] == "AB000000", "Genus_column"] <- "Wolbachia_1"

mat[mat[, "ID_column"] == "AB000001", "Genus_column"] <- "Wolbachia_2"

should do the trick (I'm not familiar with SILVA either, I assume that, like UNITE, it has an identifier somewhere you can use for searching).

Also, it's quite likely that phyloseq or any other microbiome R package already has a built-in function to do this. But, again, I don't have experience with phyloseq. So I would stick to the RESCRIPt approach.

Best wishes,

Sergio

--

Disclaimer: I'm only another forum user, just like you. Please don't take my answer as a ground truth. A Forum Moderator would probably provide you with a more accurate answer.

2 Likes
3

Hello Lily,

EDIT: Here's how I would do this
EDIT: I no longer use Phloseq. See by solution with Base R + Tidyverse below!

Phyloseq::psmelt() converts a Phyloseq object into a Tidy™ DataFrame :broom:.

Once my data is in this standard R data structure, I can edit the taxonomy and make plots with my favorite R packages.

Here's some real (and very old) example code:

Here's some (newer) code on your data:

phy_rel_abund |>
  psmelt() |>
  mutate(
    Genus = case_when(
      # fix Wolbachia
    )
  ) |>
  ggplot(aes(...)) + # Start ggplot
  geom_bar()         # Make the rest of the plot

Other notes:

Phyloseq has been abandoned for years, so I'm using Tidyverse packages instead.

Please don't take my answers as truth either! There are many ways to make a barplot. :bar_chart:

4 Likes
4

Sorry to double-post!

I've been meaning to post a modern example, so here goes!

Because I usually run qiime taxa barplot anyway, I just import that .qzv file into R.

This uses Base R + Tidyverse (no Phylose, not even qiime2r!)

Full code on GitHub:

image
image1344×960 38.8 KB

3 Likes
5

Hi. To use R microeco package to create bar plot is also a good choice. Some examples can be found here (Chapter 4 Composition-based class | Tutorial for R microeco package (v1.7.1)). The file2meco package provides a function to fast convert qza to microtable object (Chapter 8 file2meco package | Tutorial for R microeco package (v1.7.1)).

1 Like
6

This topic was automatically closed 31 days after the last reply. New replies are no longer allowed.