qiime2R code example

Example: qiime2R

> # Load packages
> library(tidyverse) 
> library(RColorBrewer)
> library(cowplot)
> library(qiime2R) 
> library(MicrobeR) 
> 
> # Load function
> source("make_taxa_barplot.R")
> 
> # Import data ##################################################################
> # Metadata
> metadata <- read_tsv("metadata.tsv", comment = "#q2") 
> 
> metadata <- metadata %>% 
>   rename_at(vars(contains("-")), ~gsub("-", "_", .x)) %>% # replace dash with underscore 
>   rename(SampleID = "sample_id") # use "SampleID" as sample identifier 
> 
> # Feature table
> table <- read_qza("table.qza")
> 
> count_tab <- table$data %>% as.data.frame() 
> 
> # Taxonomy
> taxonomy <- read_qza("taxonomy.qza")
> 
> tax_tab <- taxonomy$data %>% 
>   as.data.frame() %>%
>   separate(Taxon, sep = ";", c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")) %>% 
>   column_to_rownames("Feature.ID") %>%
>   select(-Confidence)
> 
> # Taxa barplot #################################################################
> # Phylum-level ----------------------------------------------------------------- 
> # Collapse feature table at phylum level
> tab_phy <- Summarize.Taxa(count_tab, tax_tab)$Phylum %>%
>   # the following 4 lines of code prune the taxonomy to contain phylum names only
>   rownames_to_column("tax") %>%
>   filter(!grepl("NA", tax)) %>% # remove taxa without kingdom or phylum level annotation
>   mutate(tax = gsub("\\[|\\]", "", tax), # remove square brackets
>          tax = gsub("k__.*p__", "", tax)) %>%
>   column_to_rownames("tax")
> 
> # Plot taxa on individual basis
> make_taxa_barplot(table = tab_phy, 
>                   metadata = metadata, 
>                   group_by = body_site, 
>                   ntaxa = 10, # number of taxa to display
>                   nrow = 1, # nrow for facetting
>                   plot_mean = FALSE, # plot mean relative abundance within each group
>                   cluster_sample = FALSE, # cluster samples based on similarity
>                   sample_label = SampleID,
>                   italize_taxa_name = TRUE)
> 
> ggsave("taxa_barplot_moving_picture_individual.png")