Overview

Taxonomic diversity can be measured using a variety of statistical models. Alpha diversity represents the diversity (richness, evenness, compositional complexity) within a group (e.g. an experiment group) and beta diversity (next section) examines diversity between groups.

'Diversity' is based on the amplicon sequence variants (ASV) detected by DADA2 per sample.

In this section we can estimate alpha diversity using Shannon’s index, Simpson's index, Chao1 richness and observed ASVs. Different diversity indices were used because each have different strengths and weaknesses and are better contextually - i.e. depending on the structure and nature of your dataset. For example, Shannon’s index is more sensitive to sample size and influenced by rare ASVs. Simpson's index is less sensitive to sample size and or rare ASVs, but is influenced by dominance/abundance of ASVs.

A good overview of most of these methods can be seen in this paper:

https://www.jmb.or.kr/submission/Journal/027/JMB027-12-02_FDOC_2.pdf

The plots and tables in each section (where each section = a variable) are as follows:

1. Shannon's index calculates richness and diversity using a natural logarithm and accounts for both abundance and evenness of the taxa present. Shannon's index is a universally accepted index for diversity as it accounts for entropy in an ecosystem or in representative samples.

Shannon, C.E. and Weaver, W. (1949). “The mathematical theory of communication”. University of Illonois Press, Champaign, Illonois.

2. Simpson's index is used to estimate dominance of the species but does not account for species richness.

Simpson, E. H. (1949). Measurement of diversity. nature, 163(4148), 688-688.

3. Chao1 richness is a non-parametric estimator of diversity, based on abundance.

Chao, A. (1984). Nonparametric estimation of the number of classes in a population. Scandinavian Journal of statistics, 265-270.

4. Observed ASVs calculates diversity based on the number of distinct Amplicon Sequence Variants. Originally this method was developed for observed OTUs, as seen in the below paper.

DeSantis, T.Z., Hugenholtz, P., Larsen, N., Rojas, M., Brodie, E.L., Keller, K. Huber, T., Davis, D., Hu, P., Andersen, G.L. (2006). “Greengenes, a Chimera-Checked 16S rRNA Gene Database and Workbench Compatible with ARB”. Applied and Environmental Microbiology (72): 5069–5072.

5. Combined Shannon and Observed ASV. An additional plot combing both Shannon’s index and Observed ASV indices has been included, to compare similarities and differences between these results. As each index uses different units, results for both have been normalised between 0 and 1.

6. Kruskal-Wallis rank sum test is a rank-based nonparametric test that can be used to determine if there are statistically significant differences between two or more groups. This statistical analysis is provided for each plot, to estimate if there is a significant difference (q.value < 0.05) between all groups.

7. Pairwise Wilcoxon rank sum test (AKA: Mann-Whitney test is the same as the Kruskal-Wallis test, but applied pairwise to each group (technically, The Kruskal-Wallis test is the generalization of the Wilcoxon rank sum test).

4a. Preparing your data

ampvis2 requires 3 components:

1. A table of read counts per sample per ASV (abundance table)

2. A table that matches ASVs to taxonomy information

3. A metadata table containing your sample IDs and any variables (e.g. experimental groups) associated with each sample.

...

nfcore/ampliseq generates the abundance and taxonomy tables, and we created the metadata table in the previous section, so we can import all 3 into R.

Import the metadata table

When you ran your sequences through the ampliseq pipeline, you submitted the samples with a metadata file. This file contains information on your samples and variables. We need to import this metadata file to run our analysis on selected variables.

samples_table <- read.table("../illumina/data/metadata.tsv", sep = "\t", header = T)

Have a look at your samples table and variables (metadata):

samples_table

Import the ASV abundance table

First import the unfiltered ASV table (ampvis2 does internal filtering).

asvtable <- read.table("../illumina/results/dada2/ASV_table.tsv", check.names = FALSE, sep = "\t", skip = 1, stringsAsFactors = FALSE, comment.char = "", header = T)
colnames(asvtable)[1] <- "ASV"

Have a look at the top few rows of the ASV table, to see if it looks right. The sample IDs should be the column names and the ASV IDs in the first column. All the other columns should contain numbers (i.e. the count of the number of times each ASV was found in each sample).

head(asvtable)

Import the taxonomy table

Read in the taxonomy data - i.e. the taxonomic assignments for each ASV

...

Overview

Taxonomic diversity can be measured using a variety of statistical models. Alpha diversity represents the diversity (richness, evenness, compositional complexity) within a group (e.g. an experiment group) and beta diversity (next section) examines diversity between groups.

'Diversity' is based on the amplicon sequence variants (ASV) detected by DADA2 per sample.

In this section we can estimate alpha diversity using Shannon’s index, Simpson's index, Chao1 richness and observed ASVs. Different diversity indices have different strengths and weaknesses and are better contextually - i.e. depending on the structure and nature of your dataset. For example, Shannon’s index is more sensitive to sample size and influenced by rare ASVs. Simpson's index is less sensitive to sample size and or rare ASVs, but is influenced by dominance/abundance of ASVs.

A good overview of most of these methods can be seen in this paper:

https://www.jmb.or.kr/submission/Journal/027/JMB027-12-02_FDOC_2.pdf

The plots and tables in each section (where each section = a variable) are as follows:

1. Shannon's index calculates richness and diversity using a natural logarithm and accounts for both abundance and evenness of the taxa present. Shannon's index is a universally accepted index for diversity as it accounts for entropy in an ecosystem or in representative samples.

Shannon, C.E. and Weaver, W. (1949). “The mathematical theory of communication”. University of Illonois Press, Champaign, Illonois.

2. Simpson's index is used to estimate dominance of the species but does not account for species richness.

Simpson, E. H. (1949). Measurement of diversity. nature, 163(4148), 688-688.

3. Chao1 richness is a non-parametric estimator of diversity, based on abundance.

Chao, A. (1984). Nonparametric estimation of the number of classes in a population. Scandinavian Journal of statistics, 265-270.

4. Observed ASVs calculates diversity based on the number of distinct Amplicon Sequence Variants. Originally this method was developed for observed OTUs, as seen in the below paper.

DeSantis, T.Z., Hugenholtz, P., Larsen, N., Rojas, M., Brodie, E.L., Keller, K. Huber, T., Davis, D., Hu, P., Andersen, G.L. (2006). “Greengenes, a Chimera-Checked 16S rRNA Gene Database and Workbench Compatible with ARB”. Applied and Environmental Microbiology (72): 5069–5072.

5. Combined Shannon and Observed ASV. An additional plot combing both Shannon’s index and Observed ASV indices has been included, to compare similarities and differences between these results. As each index uses different units, results for both have been normalised between 0 and 1.

6. Kruskal-Wallis rank sum test is a rank-based nonparametric test that can be used to determine if there are statistically significant differences between two or more groups. This statistical analysis is provided for each plot, to estimate if there is a significant difference (q.value < 0.05) between all groups.

7. Pairwise Wilcoxon rank sum test (AKA: Mann-Whitney test is the same as the Kruskal-Wallis test, but applied pairwise to each group (technically, The Kruskal-Wallis test is the generalization of the Wilcoxon rank sum test).

4a. Preparing your data

ampvis2 requires 3 components:

1. A table of read counts per sample per ASV (abundance table)

2. A table that matches ASVs to taxonomy information

3. A metadata table containing your sample IDs and any variables (e.g. experimental groups) associated with each sample.

...

nfcore/ampliseq generates the abundance and taxonomy tables, and we created the metadata table in the previous section, so we can import all 3 into R.

Import the metadata table

When you ran your sequences through the ampliseq pipeline, you submitted the samples with a metadata file. This file contains information on your samples and variables. We need to import this metadata file to run our analysis on selected variables.

samples_table <- read.table("../illumina/results/dada2/ASV_tax.silva_138data/metadata.tsv", sep = "\t", header = T, quoteheader = ""T)

Have a look again at the top 6 rowsat your samples table and variables (metadata):

head(mytax)

Combine the taxonomy and ASV tables

asv_table <- cbind(asvtable, mytaxsplit)
# Also remove the ASV_ID, Confidence and sequence columns, as they are not needed
asv_table <- subset(asv_table, select=-c(ASV_ID, Confidence, sequence))

head(asv_table)

...

asv_table$Kingdom <- gsub("D_0", "k", asv_table$Kingdom)
asv_table$Phylum <- gsub("D_1", "p", asv_table$Phylum)
asv_table$Class <- gsub("D_2", "c", asv_table$Class)
asv_table$Order <- gsub("D_3", "o", asv_table$Order)
asv_table$Family <- gsub("D_4", "f", asv_table$Family)
asv_table$Genus <- gsub("D_5", "g", asv_table$Genus)
asv_table$Species <- gsub("D_6", "s", asv_table$Speciessamples_table

Import the ASV abundance table

First import the unfiltered ASV table (ampvis2 does internal filtering).

asvtable <- read.table("../illumina/results/dada2/ASV_table.tsv", check.names = FALSE, sep = "\t", stringsAsFactors = FALSE, comment.char = "", header = T)
colnames(asvtable)[1] <- "ASV"

Have a look at the top few rows of the ASV table, to see if it looks right. The sample IDs should be the column names and the ASV IDs in the first column. All the other columns should contain numbers (i.e. the count of the number of times each ASV was found in each sample).

head(asvtable)

Import the taxonomy table

Read in the taxonomy data - i.e. the taxonomic assignments for each ASV

mytax <- read.table("../illumina/results/dada2/ASV_tax.silva_138.tsv", sep = "\t", header = T, quote = "")

Have a look again at the top 6 rows

head(mytax)

Combine the taxonomy and ASV tables

asv_table <- cbind(asvtable, mytax)
# Also remove the ASV_ID, Confidence and sequence columns, as they are not needed
asv_table <- subset(asv_table, select=-c(ASV_ID, confidence, sequence))

head(asv_table)

Convert the imported data to ampvis2 format

...

pdf_exp <- paste0("rarefaction_", group, "_samples.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")

You can now find these files in your working directory (which you originally defined in the 'Setting up your analysis environment' section).

...

Briefly, the above code changes the theme (theme_bw()), the Y axis label (ylab(indicname) - note: the current y axis name is inside a variable called indicname, so you need to replace it with text in quotes, e.g. "Chao richness"), removes the X axis label (xlab("")), flips the X and Y axis (coord_flip()), changes the font size to 17 and removes the legend (both in theme()). Finally, a trend line is added (geom_smooth()). Any of these components can be removed or modified as you like.

You can save Save your plot as a 300dpi (i.e. publication quality) tiff or pdf file. These files can be found in your working directory.Export as a 300dpi tifftiff and pdf file.

tiff_exp <- paste0("alpha_div_box_plot_", group, "_", indicname, "_samples.tiff")
ggsave(file = tiff_exp, dpi = 300, compression = "lzw", device = "tiff", plot = p, width = 20, height = 20, units = "cm")

Export as a pdf

pdf_exp <- paste0("alpha_div_box_plot_", group, "_", indicname,"_samples.pdf")
ggsave(file = pdf_exp, device = "pdf", plot = p, width = 20, height = 20, units = "cm")

...

To compare the differences between groups within your variable, a Kruskal-Wallis test (one-way analysis of variance) can be performed to test for overall differences between all groups, and a Wilcoxon rank sum to test pairwise differences between each group.

...