...

This will provide you a summary of sequence reads per sample (minimum, maximum, average, etc) and the number and percentage of ASVs resolved to each taxonomic level.

6b. Taxonomic abundance bar plots

This section generates stacked bar plots showing the proportion of taxa per sample, for each taxonomic level.

First, you can output the raw ASV abundance data as a csv file. You can add this as supplemental data to a manuscript or use it in other plotting or statistical programs:

...

tiff_exp <- paste0("community_tax_abund_bar_plot_", taxgroup, "_samples_.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")
pdf_exp <- paste0("community_tax_abund_bar_plot_", taxgroup, "_samples_.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")

6c. Taxonomic abundance box and whisker plots

This section contains box and whisker plots for each taxonomic level (phylum … species).

As with the bar plots in the previous section, these B&W plots show the relative abundance (in %) of taxa. But whereas heatmaps show the mean %, box and whisker plots show the interquartile range between samples (max, min, median, quartiles), thus providing additional statistical information.

Choose the taxonomic level you want to plot (Choose from Phylum, Class, Order, Family, Genus, Species)

taxgroup <- "Order"

Generate the plot (and modify the attributes as you like)

Note the tax_show = 10 attribute. This says to show the top 10 taxa. Change this as desired.

p <- amp_boxplot(ampvisdata, tax_aggregate = taxgroup, tax_empty = "remove", tax_show = 10) +
theme_bw() +
theme(text = element_text(size = 18), axis.text.x = element_text(size=16), axis.text.y = element_text(size=14)) +
labs(x = taxgroup)
p

Export the plots as tiff and pdf:

tiff_exp <- paste0("community_tax_abund_BW_plot_", taxgroup, "_samples_.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")
pdf_exp <- paste0("community_tax_abund_BW_plot_", taxgroup, "_samples_.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")

6d. Choosing a variable to analyse

In the Alpha Diversity section, you already imported your nfcore/ampliseq results and converted them to an ampvis2 object. To have another quick look at your ampvis object:

ampvisdata

 As in the Alpha Diversity section, you need to choose a variable to work with.

You can view your variables as column names in your samples_table:

colnames(samples_table)

 Now enter the variable you want to analyse:

group <- "Nose_size"

And once again you’ll need to order the groups in your variable. See the ‘Ordering your variable’ section in the Alpha diversity section for more details.

Choose how you want to order your groups:

lev <- c("Small", "Medium", "Big")

Then run the following to apply the levels to your data:

ampvisdata$metadata[[group]] <- factor(ampvisdata$metadata[[group]], levels = lev)

6e. Community structure by variable

In this section you can generate heatmaps and box and whisker plots for the variable you selected above, and for each taxonomic level.

If you want to analyse another variable here, go back to the previous section, choose another variable and re-run the code from that point.

Choose the taxonomic level you want to plot (Choose from Phylum, Class, Order, Family, Genus, Species). This applies to both the heatmaps and B&W plots below.

taxgroup <- "Order"

Heatmaps

This section contains heatmaps for each taxonomic level (class - species) for your selected variable. The number in each cell is the relative proportion (in %) of taxa per group.

There are a variety of attributes you can modify:

Note the tax_show = 10 attribute. This says to show the top 10 taxa.

The plot_values_size = defines the size of the text in the heatmap cells.

The tax_add = NULL adds an additional taxonomic group to the taxa names (e.g. changing this to tax_add = "Phylum" and plotting Genus will name the taxa as 'phylum:genus' but leaving it as tax_add = NULL will just give the genus name.

showRemainingTaxa = T will aggregate, in a single row, the remaining taxa not shown on this heatmap. Change to F if you don't want to see this.

color_vector = NULL uses the default colour range (orange -> blue). You can change this by providing your own colour range, e.g. color_vector = c("red", "green"). You can choose from the huge number of R colours here: http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf

Generate the plot:

p <- amp_heatmap(ampvisdata, tax_show = 10, plot_values_size = 4, tax_add = NULL, showRemainingTaxa = T, color_vector = NULL, group_by = group, tax_aggregate = taxgroup, tax_empty = "remove")
p

Add some additional plot modifications. Change or remove these as desired (or add your own modifications)

p <- p +
theme_bw() +
theme(text = element_text(size = 18), axis.text.x = element_text(angle = 90, size=16), axis.text.y = element_text(size=14)) +
labs(y = taxgroup)
p

Export as pdf and tiff

tiff_exp <- paste0("community_tax_abund_heatmap_", group, "_", taxgroup, "_samples_.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")
pdf_exp <- paste0("community_tax_abund_heatmap_", group, "_", taxgroup, "_samples_.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")

Box and whisker plots

This section contains B&W plots for each taxonomic level (class - species) for your selected variable. This is similar to the heatmaps you generated in section 6, but separates the B&W into groups per variable, thus enabling an examination of how the variable groups differ in terms of taxonomic abundance.

Again, there are a variety of attributes you can modify:

Note the tax_show = 7 attribute. This says to show the top 7 taxa (Tip: to look at specific taxa you can provide a vector of taxa names, e.g. tax_show = c("Bacteroidales", "Fibrobacterales")).

The tax_add = NULL adds an additional taxonomic group to the taxa names (e.g. changing this to tax_add = "Phylum" and plotting Genus will name the taxa as 'phylum:genus' but leaving it as tax_add = NULL will just give the genus name.

Generate the plot:

p <- amp_boxplot(ampvisdata, tax_show = 7, tax_add = NULL, tax_aggregate = taxgroup, tax_empty = "remove", group_by = group)
p

You can change additional properties of the plot here:

p$mapping$fill <- as.name(".Group")
p <- p + theme_bw() +
scale_color_manual(values=c("Red", "Green", "Blue")) + scale_fill_manual(values=c("Red", "Green", "Blue")) +
labs(fill = group, x = taxgroup) +
theme(text = element_text(size = 18), axis.text.x = element_text(size=18), axis.text.y = element_text(size=16)) +
guides(color = "none")
p

Export as pdf and tiff

tiff_exp <- paste0("community_tax_abund_BW_", group, "_", taxgroup, "_samples_.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")
pdf_exp <- paste0("community_tax_abund_BW_", group, "_", taxgroup, "_samples_.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")

...