0

I am trying to do a relative abundance bar plot, but the bars in the graph have no colors despite the legends having the colors. Below is the script that I used. Kindly assist. I have pasted a photo of how the graph looks like.

library(phyloseq)
library(vegan)
library(ggplot2)
require(vegan)
library(gplots)
library(dplyr)
library(scales)

##create the phyloseq object (I just give the basic command here, but don't forget to rarefy)

table = read.table (file.choose(), header = TRUE, row.names = 1 , sep = ",")
taxonomy = read.table (file.choose(), header = TRUE, row.names = 1 , sep = ",")
metadata = read.table (file.choose(), header = TRUE, row.names = 1 , sep = ",")

physeq_OTU = otu_table(as.matrix(table), taxa_are_rows = T)
physeq_tax = tax_table(as.matrix(taxonomy))
physeq_metadata = sample_data(metadata)

physeq <- phyloseq(physeq_OTU, physeq_tax, physeq_metadata)

physeq_Archaea<- subset_taxa(physeq, Kingdom=="d__Archaea")
physeq_Bacteria<- subset_taxa(physeq, Kingdom=="d__Bacteria")

Archaea=subset_taxa(physeq, Kingdom=="d__Archaea")
Bacteria=subset_taxa(physeq, Kingdom=="d__Bacteria")

#subset to phylum table for the specimens of interest
Phylum_Archaea = tax_glom(Archaea, taxrank = "Phylum")

#create dataframe from phyloseq object
data_glom<- psmelt(Phylum_Archaea)

#convert to character
data_glom$Phylum <- as.character(data_glom$Phylum)

#rename genera with <1% abundance
data_glom$Phylum[data_glom$Abundance < 0.01] <- "<1% abundance"

write.csv(data_glom, file="data_glom_Archaea.csv")

#Count #genera to set colour palette
Count = length(unique(data_glom$Phylum))
Count

#To change order of Archaea, use factor() to change order of levels
#maybe you want to put the most abundant taxa first
#Use unique(data_glom$class) to see exact spelling of each phylum
unique(data_glom$Phylum)
data_glom$Phylum <- factor(data_glom$Phylum, levels = c("p__Halobacterota", 
"p__Thermoplasmatota", "p__Nanoarchaeota", "p__Crenarchaeota", "< 1% abund."))

 Barplot <- ggplot(data=data_glom, aes(x=location, y=Abundance, fill = Phylum)) 
 Barplot2 = Barplot + geom_bar(aes(), stat="identity", position ="fill", color = "black") +
 scale_fill_manual(values = c("p__Halobacterota" = "green", "p__Thermoplasmatota" = 
 "aquamarine", "p__Nanoarchaeota" = "goldenrod1", "p__Crenarchaeota" = "violetred1", "< 1% 
 abund." = "dark cyan")) +
 theme(legend.position="right") + 
 guides(fill=guide_legend(nrow=6)) + 
 theme(text = element_text(size = 10)) + 
 theme(axis.text.x=element_text(angle=90, hjust=1)) +
 scale_y_continuous(labels = scales::percent_format(scale = 100))+ theme_classic()

enter image description here

Quinten
  • 35,235
  • 5
  • 20
  • 53
  • 2
    Your code is not reproducibe since we don't have access to your data. Nevertheless, I believe the `values` in your `scale_fill_manual` should be the colours you want to appear on your plot, not the *labels* you want to associate with each. Check the online help for details. Eg `p + scale_colour_manual(values = c("red", "blue", "green"))` – Limey Jul 22 '22 at 12:05
  • @Limey The approach with using a named vector is fine. That has no effect on the labels (which are set via the `labels` argument). But doing so ensures that colors are assigned to the right categories, i.e. you don't have to think about putting the colors in the right order. Try e.g. `ggplot(mtcars, aes(hp, mpg, color = factor(cyl))) + geom_point() + scale_color_manual(values = c("8" = "blue", "4" = "red", "6" = "yellow"))`. – stefan Jul 22 '22 at 14:18
  • The most likely reason is that something went wrong when converting `Phylum` to a factor. I would guess that by accident you converted your `Phylum`s to NAs and hence when filling by Phylumn your bars are filled according to the `na.value` of scale_fill_manual which is `"grey50"`. – stefan Jul 22 '22 at 14:21
  • Hi @Limey, thank you for your kind response. here is the link to my data to reproduce the code. https://drive.google.com/drive/folders/1kYNSYe0t3-sDfHUMld1qqzgXt9dYkWe- – Nghalipo Elise Jul 22 '22 at 19:58
  • Hi @stefan, thanks for your response, I tried the code you suggested but it also didnt work. Here is the link to my data too, thanks: https://drive.google.com/drive/folders/1kYNSYe0t3-sDfHUMld1qqzgXt9dYkWe- – Nghalipo Elise Jul 22 '22 at 20:01
  • The code was not meant to offer a solution to your issue. Just an example that it is fine to use a named vector (; As I said: My guess is that something went wrong when converting `Phylum` to a `factor`, e.g. you recode values below 0.01 as `"<1% abundance"`. But when converting to a factor you use `"< 1% abund."`. As a consequence at least `"<1% abundance"` gets converted to NAs. – stefan Jul 22 '22 at 20:11
  • 1
    Also. I would suggest to have a look at how to provide [a minimal reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). I and most people don't want to first download three files and running a lot of code to help you fix your issue when only a snippet of your dataset `data_glom` is needed. See the section `Providing a minimal dataset` in the link I referenced on how to do this. – stefan Jul 22 '22 at 20:15
  • Hi @stefan, well noted with thanks. – Nghalipo Elise Jul 22 '22 at 21:00

1 Answers1

0

I was able to fix my issue with the code below..there were a couple of lines in my previous code that were converting my phylum column to NAs.

Phylum_Archaea = tax_glom(Archaea, taxrank = "Phylum")

#Transform sample counts
glom_Phylum_Archaea = transform_sample_counts(Phylum_Archaea, function(x) x / 
sum(x))

#create dataframe from phyloseq object
data_glom<- psmelt(glom_Phylum_Archaea)

Barplot <- ggplot(data=data_glom, aes(x=location, y=Abundance, fill = Phylum)) 
Barplot2 = Barplot + geom_bar(aes(), stat="identity", position ="fill") +
scale_fill_manual(values = c("green", "aquamarine", "goldenrod1", "violetred1")) 
+
theme(legend.position="right") + 
guides(fill=guide_legend(nrow=6)) + 
theme(text = element_text(size = 10)) + 
theme(axis.text.x=element_text(angle=90, hjust=1)) +
scale_y_continuous(labels = scales::percent_format(scale = 100))+ theme_classic ()