Heatmap of core microbiome using qiime2R

Other Bioinformatics Tools
,
1

Hello,

I'm trying to build a heatmap using data that was imported using qiime2R -- I'm not quite sure if this is the place to post this question, but I really don't know how to solve this. -- Please see below:

Screen Shot 2022-10-04 at 4.11.57 PM
Screen Shot 2022-10-04 at 4.11.57 PM1378×844 139 KB

The problem is that the label of the y-axis has the "Feature ID" instead of the taxonomy of the organisms. How can I change it?

Thanks!
Giovana

2

Hi @GioFranco,

Let's try something, you'll need to import your taxonomy.qza and feature-table.qza and your sample metadata (in this example called metadata)

library(tidyverse)
library(zoo)
library(qiime2R)

#this chunk only needed if you're going to use the metadata in your analysis
metadata <- read_tsv("sample-metadata.tsv") 
meta <- metadata %>%
  dplyr::arrange(sample_name) #or replace with whatever your first column is called

#import feature-table
asvtab <- qiime2R::read_qza("table.qza")$data %>% # the feature-table
  t() %>% #transpose
  as.data.frame() #make into data.frame

#import taxonomy table
bespoke_taxonomy0 <- qiime2R::read_qza("gtdb89-taxonomy.qza")$data %>% 
  filter(`Feature.ID` %in% colnames(asvtab)) %>% #for sanity just in case
  arrange(Feature.ID)

#rearrange order of the two tables so they are in the same order
asvtab <- asvtab[bespoke_taxonomy0$Feature.ID] #safety step

# Separate the ranks in our taxonomy. The separate below is meant for GTDB. Change accordingly for whatever database you used.  

aa <- bespoke_taxonomy0 %>% 
  tidyr::separate(Taxon, sep = ";", c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species"))

# Not all ASVs have complete taxonomy down to species level. So instead I'm going to grab whatever was the last non-NA value we had between Domain:Species.

aa2 <- aa %>%
  dplyr::select(Feature.ID, Domain:Species) %>% 
  dplyr::mutate(purrr::pmap_df(., ~ zoo::na.locf(c(...)[-1]))) %>% 
  dplyr::left_join(select(aa, !Domain:Species), aa, by="Feature.ID")

# And just in case there are ASVs with identical species name, going to add a suffix to distinguish these from each other

aa3 <- aa2 %>% 
  dplyr::mutate(Species2=make.unique(as.character(aa2$Species), sep = "_ASV"))

Now you have a new object with the highest taxonomic name for each ASV, instead of the hashed ID.
You can use it as you normally would, or if you want to change the names in your feature-table you can do so with:

colnames(asvtab) <- aa3$Species2

Hope that works!

1 Like