{"type":"doc","content":[{"type":"extension","attrs":{"layout":"default","extensionType":"com.atlassian.confluence.macro.core","extensionKey":"toc","parameters":{"macroParams":{"minLevel":{"value":"1"},"maxLevel":{"value":"6"},"outline":{"value":"false"},"style":{"value":"none"},"type":{"value":"list"},"printable":{"value":"true"}},"macroMetadata":{"macroId":{"value":"a102598c-63b1-4424-bb57-9fc2ed865b05"},"schemaVersion":{"value":"1"},"title":"Table of Contents"}},"localId":"5cafbe00-e748-451c-b94f-fdd9f20ac3e3"}},{"type":"paragraph"},{"type":"heading","attrs":{"level":1},"content":[{"text":"Overview","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Taxonomic diversity can be measured using a variety of statistical models. ","type":"text"},{"text":"Alpha diversity","type":"text","marks":[{"type":"strong"}]},{"text":" represents the diversity (richness, evenness, compositional complexity) ","type":"text"},{"text":"within","type":"text","marks":[{"type":"strong"}]},{"text":" a group (e.g. an experiment group) and beta diversity (next section) examines diversity ","type":"text"},{"text":"between ","type":"text","marks":[{"type":"strong"}]},{"text":"groups.","type":"text"}]},{"type":"paragraph","content":[{"text":"'Diversity' is based on the ","type":"text"},{"text":"amplicon sequence variants (ASV)","type":"text","marks":[{"type":"link","attrs":{"href":"https://en.wikipedia.org/wiki/Amplicon_sequence_variant"}}]},{"text":" detected by DADA2 per sample.","type":"text"}]},{"type":"paragraph","content":[{"text":"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. ","type":"text"}]},{"type":"paragraph","content":[{"text":"A good overview of most of these methods can be seen in this paper:","type":"text"}]},{"type":"paragraph","content":[{"text":"https://www.jmb.or.kr/submission/Journal/027/JMB027-12-02_FDOC_2.pdf","type":"text","marks":[{"type":"link","attrs":{"href":"https://www.jmb.or.kr/submission/Journal/027/JMB027-12-02_FDOC_2.pdf"}}]}]},{"type":"paragraph","content":[{"text":"The plots and tables in each section (where each section = a variable) are as follows:","type":"text"}]},{"type":"orderedList","attrs":{"order":1},"content":[{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Shannon's index","type":"text","marks":[{"type":"strong"},{"type":"link","attrs":{"href":"http://www.tiem.utk.edu/~gross/bioed/bealsmodules/shannonDI.html"}}]},{"text":" 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.","type":"text"}]}]}]},{"type":"paragraph","content":[{"text":"Shannon, C.E. and Weaver, W. (1949). “The mathematical theory of communication”. University of Illonois Press, Champaign, Illonois.","type":"text","marks":[{"type":"code"}]}]},{"type":"orderedList","attrs":{"order":2},"content":[{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Simpson's index","type":"text","marks":[{"type":"strong"},{"type":"link","attrs":{"href":"http://www.tiem.utk.edu/~gross/bioed/bealsmodules/simpsonDI.html"}}]},{"text":" is used to estimate dominance of the species but does not account for species richness.","type":"text"}]}]}]},{"type":"paragraph","content":[{"text":"Simpson, E. H. (1949). Measurement of diversity. nature, 163(4148), 688-688.","type":"text","marks":[{"type":"code"}]}]},{"type":"orderedList","attrs":{"order":3},"content":[{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Chao1 richness","type":"text","marks":[{"type":"strong"},{"type":"link","attrs":{"href":"http://www.pisces-conservation.com/sdrhelp/index.html?chaoquant.htm"}}]},{"text":" is a non-parametric estimator of diversity, based on abundance.","type":"text"}]}]}]},{"type":"paragraph","content":[{"text":"Chao, A. (1984). Nonparametric estimation of the number of classes in a population. Scandinavian Journal of statistics, 265-270.","type":"text","marks":[{"type":"code"}]}]},{"type":"orderedList","attrs":{"order":4},"content":[{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Observed ASVs","type":"text","marks":[{"type":"strong"},{"type":"link","attrs":{"href":"https://www.nature.com/articles/ismej2017119"}}]},{"text":" 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.","type":"text"}]}]}]},{"type":"paragraph","content":[{"text":"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.","type":"text","marks":[{"type":"code"}]}]},{"type":"orderedList","attrs":{"order":5},"content":[{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Combined Shannon and Observed ASV","type":"text","marks":[{"type":"strong"}]},{"text":". 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.","type":"text"}]}]},{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Kruskal-Wallis rank sum test","type":"text","marks":[{"type":"strong"},{"type":"link","attrs":{"href":"https://statistics.laerd.com/stata-tutorials/kruskal-wallis-h-test-using-stata.php"}}]},{"text":" 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 ","type":"text"},{"text":"all groups","type":"text","marks":[{"type":"strong"}]},{"text":".","type":"text"}]}]},{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Pairwise Wilcoxon rank sum test (AKA: Mann-Whitney test","type":"text","marks":[{"type":"strong"},{"type":"link","attrs":{"href":"https://sphweb.bumc.bu.edu/otlt/mph-modules/bs/bs704_nonparametric/bs704_nonparametric4.html"}}]},{"text":" 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).","type":"text"}]}]}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":1},"content":[{"text":"4a. Preparing your data","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"ampvis2 requires 3 components:","type":"text"}]},{"type":"orderedList","attrs":{"order":1},"content":[{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"A table of read counts per sample per ASV (abundance table)","type":"text"}]}]},{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"A table that matches ASVs to taxonomy information","type":"text"}]}]},{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"A metadata table containing your sample IDs and any variables (e.g. experimental groups) associated with each sample.","type":"text"}]}]}]},{"type":"mediaSingle","attrs":{"layout":"center","width":250,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":400,"alt":"image-20240530-200611.png","id":"912d7f60-9409-4827-94ba-b491bc37bb59","collection":"contentId-2593129710","type":"file","height":463}}]},{"type":"paragraph","content":[{"text":"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. ","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":2},"content":[{"text":"Import the metadata table","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"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.","type":"text"}]},{"type":"codeBlock","content":[{"text":"samples_table <- read.table(\"../data/illumina/fastq/metadata.tsv\", sep = \"\\t\", header = T)","type":"text"}]},{"type":"paragraph","content":[{"text":"Have a look at your samples table and variables (metadata):","type":"text"}]},{"type":"codeBlock","content":[{"text":"samples_table","type":"text"}]},{"type":"heading","attrs":{"level":2},"content":[{"text":"Import the ASV abundance table","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"First import the ","type":"text"},{"text":"unfiltered","type":"text","marks":[{"type":"strong"}]},{"text":" ASV table (ampvis2 does internal filtering).","type":"text"}]},{"type":"codeBlock","content":[{"text":"asvtable <- read.table(\"../runs/run2_ampliseq/results/dada2/ASV_table.tsv\", check.names = FALSE, sep = \"\\t\", stringsAsFactors = FALSE, comment.char = \"\", header = T)\ncolnames(asvtable)[1] <- \"ASV\"","type":"text"}]},{"type":"paragraph","content":[{"text":"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).","type":"text"}]},{"type":"codeBlock","content":[{"text":"head(asvtable)","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":2},"content":[{"text":"Import the taxonomy table","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Read in the taxonomy data - i.e. the taxonomic assignments for each ASV","type":"text"}]},{"type":"codeBlock","content":[{"text":"mytax <- read.table(\"../runs/run2_ampliseq/results/dada2/ASV_tax.silva_138.tsv\", sep = \"\\t\", header = T, quote = \"\")","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Have a look again at the top 6 rows","type":"text"}]},{"type":"codeBlock","content":[{"text":"head(mytax)","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Combine the taxonomy and ASV tables","type":"text"}]},{"type":"codeBlock","content":[{"text":"asv_table <- cbind(asvtable, mytax)\n# Also remove the ASV_ID, Confidence and sequence columns, as they are not needed\nasv_table <- subset(asv_table, select=-c(ASV_ID, confidence, sequence))","type":"text"}]},{"type":"codeBlock","content":[{"text":"head(asv_table)","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":2},"content":[{"text":"Convert the imported data to ampvis2 format","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Now combine the samples data with the ASV table using amp_load(). This creates an ampvis2 database that can be used by ampvis2","type":"text"}]},{"type":"codeBlock","content":[{"text":"ampvisdata <- amp_load(otutable = asv_table, metadata = samples_table)","type":"text"}]},{"type":"rule"},{"type":"paragraph","content":[{"type":"emoji","attrs":{"id":"atlassian-red_star","text":":red_star:","shortName":":red_star:"}},{"text":" ","type":"text"},{"text":"Advanced: After running this line of code you might get a Warning Message about ‘Only 37 of 59 unique sample names match between metadata and otutable. The following unmatched samples have been removed:….’ ","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff5630"}}]}]},{"type":"paragraph","content":[{"text":"Q1. Why might that be?","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff5630"}}]}]},{"type":"paragraph","content":[{"text":"Q2. What could have caused it?","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff5630"}}]}]},{"type":"paragraph","content":[{"text":"Q3. How might you investigate to see if you were right about the reason some samples were removed?","type":"text","marks":[{"type":"textColor","attrs":{"color":"#ff5630"}}]}]},{"type":"expand","attrs":{"title":"Answer for Q1:"},"content":[{"type":"paragraph","content":[{"text":"Maybe some samples sequenced poorly and there were not enough reads to pass the default filtering levels of the Nextflow ampliseq pipeline. For more information about these parameters, see here: ","type":"text"},{"type":"inlineCard","attrs":{"url":"https://nf-co.re/ampliseq/2.6.1/parameters/"}},{"text":" ","type":"text"}]}]},{"type":"expand","attrs":{"title":"Answer for Q2:"},"content":[{"type":"paragraph","content":[{"text":"Naturally low diversity levels, poor sample collection, poor sample processing, something going wrong in the lab, etc ","type":"text"}]}]},{"type":"expand","attrs":{"title":"Answer for Q3:"},"content":[{"type":"paragraph","content":[{"text":"You could look at the raw fastq files on the HPC and see how many sequences were in each of them and see if the samples being removed had substantially lower read numbers. This line of code would give you a quick summary of read numbers per fastq file:","type":"text"}]},{"type":"codeBlock","content":[{"text":"for l in `ls *.fastq.gz`; do echo ${l}; zcat ${l} | grep \"^@\" | wc -l; done","type":"text"}]},{"type":"paragraph","content":[{"text":"You could also look to see what was removed at the cutadapt step in the Nextflow pipeline, by looking at the ","type":"text"},{"text":"Illumina11.trimmed.cutadapt.log","type":"text","marks":[{"type":"code"}]},{"text":" and comparing it to the ","type":"text"},{"text":"Illumina52.trimmed.cutadapt.log","type":"text","marks":[{"type":"code"}]},{"text":" in the H:/workshop/2025/S1W1/metagenomics/runs/run2_ampliseq/results/cutadapt folder","type":"text"}]}]},{"type":"rule"},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"You can see information about the ampvis2 object you just created by typing its name","type":"text"}]},{"type":"codeBlock","content":[{"text":"ampvisdata","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":1},"content":[{"text":"4b. Choosing a variable to analyse","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"In your metadata you'll usually have multiple variables. These need to be analysed individually, by selecting the variable in this section, then running the remaining analysis sections on this chosen variable. You can then re-run the analysis on another variable by returning to this section, changing the variable name, then running again the remaining analysis sections.","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"You can view your variables as column names in your samples_table:","type":"text"}]},{"type":"codeBlock","content":[{"text":"colnames(samples_table)","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Enter the column name of the variable you want to analyse (i.e. change ","type":"text"},{"text":"group <- \"Myvariable\"","type":"text","marks":[{"type":"code"}]},{"text":" in the below cell to your chosen variable's column name). This has to be exactly the same as the column name, including capitalisation, characters such as underscores, etc.","type":"text"}]},{"type":"codeBlock","content":[{"text":"group <- \"Nose_size\"","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":2},"content":[{"text":"Ordering your variable","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"The plotting done in ampvis2 is done by the ","type":"text"},{"text":"ggplot2","type":"text","marks":[{"type":"link","attrs":{"href":"https://ggplot2.tidyverse.org/"}}]},{"text":" package. ggplot ","type":"text"},{"text":"factorises","type":"text","marks":[{"type":"link","attrs":{"href":"https://www.datamentor.io/r-programming/factor/"}}]},{"text":" variables and automatically orders them on the plot by alphabetical order.","type":"text"}]},{"type":"paragraph","content":[{"text":"This can cause your groups to be ordered incorrectly on the plot axes (e.g. a time series may not be plotted sequentially).","type":"text"}]},{"type":"paragraph","content":[{"text":"You can manually set the order of your variable here.","type":"text"}]},{"type":"paragraph","content":[{"text":"First have a look at how ggplot will order your variable.","type":"text"}]},{"type":"codeBlock","content":[{"text":"levels(factor(ampvisdata$metadata[[group]]))","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"If these are in the order you want to see them on your plot axes, nothing needs to be done. If they are in the wrong order you need to order them manually by setting the ","type":"text"},{"text":"levels","type":"text","marks":[{"type":"strong"},{"type":"link","attrs":{"href":"https://www.datamentor.io/r-programming/factor/"}}]},{"text":".","type":"text"}]},{"type":"paragraph","content":[{"text":"Choose how you want to order your groups here:","type":"text"}]},{"type":"codeBlock","content":[{"text":"lev <- c(\"Small\", \"Medium\", \"Big\")","type":"text"}]},{"type":"paragraph","content":[{"text":"To order your variable you need to put ","type":"text"},{"text":"all","type":"text","marks":[{"type":"strong"}]},{"text":" the variable levels into the ","type":"text"},{"text":"lev = c(..)","type":"text","marks":[{"type":"code"}]},{"text":". Make sure each level is in double quotes and separated by a comma.","type":"text"}]},{"type":"paragraph","content":[{"text":"Then run the following to apply the levels to your data:","type":"text"}]},{"type":"codeBlock","content":[{"text":"ampvisdata$metadata[[group]] <- factor(ampvisdata$metadata[[group]], levels = lev)","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":1},"content":[{"text":"4c. Rarefaction curve","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"This section is for plotting rarefaction curves for your samples, coloured by your chosen variable (if you want to change variables, go back and re-run section 4b, choosing a different variable).","type":"text"}]},{"type":"paragraph","content":[{"text":"Use the ampvis2 function ","type":"text"},{"text":"amp_rarecurve()","type":"text","marks":[{"type":"code"}]},{"text":" to generate the rarefaction curve plot.","type":"text"}]},{"type":"codeBlock","content":[{"text":"p <- amp_rarecurve(ampvisdata, color_by = group, stepsize = 200)","type":"text"}]},{"type":"paragraph","content":[{"text":"You can print this out as-is simply by:","type":"text"}]},{"type":"codeBlock","content":[{"text":"p","type":"text"}]},{"type":"heading","attrs":{"level":2},"content":[{"text":"Modifying plot attributes","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"You can make additional modifications to the plot colours, axis labels, font size, theme, etc:","type":"text"}]},{"type":"codeBlock","content":[{"text":"p <- p + scale_color_manual(values=c(\"red\", \"green\", \"blue\", \"orange\", \"black\", \"forestgreen\", \"firebrick\", \"goldenrod\", \"darkviolet\", \"olivedrab\", \"purple\", \"tomato\", \"yellow\", \"cornflowerblue\")) +\ntheme_bw() +\nylab(\"Number of observed ASVs\") +\ntheme(text = element_text(size = 14))\np","type":"text"}]},{"type":"paragraph","content":[{"text":"You can edit the above modifications and then re-run the code cell to generate a plot of your liking. The modifications in the above code are:","type":"text"}]},{"type":"paragraph","content":[{"text":"Colours - ","type":"text"},{"text":"scale_color_manual()","type":"text","marks":[{"type":"code"}]},{"text":". See the available R colours here: ","type":"text"},{"text":"http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf","type":"text","marks":[{"type":"link","attrs":{"href":"http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf"}}]},{"text":". If you set colours manually, you need to enter as many colours (or more, never fewer) as there are categories in your variable.","type":"text"}]},{"type":"paragraph","content":[{"text":"Theme - ","type":"text"},{"text":"theme_bw()","type":"text","marks":[{"type":"code"}]},{"text":". There are many themes available, seen here: ","type":"text"},{"type":"inlineCard","attrs":{"url":"https://ggplot2.tidyverse.org/reference/ggtheme.html"}},{"text":" . You can change ","type":"text"},{"text":"theme_bw()","type":"text","marks":[{"type":"code"}]},{"text":" to ","type":"text"},{"text":"theme_minimal()","type":"text","marks":[{"type":"code"}]},{"text":" or ","type":"text"},{"text":"theme_light()","type":"text","marks":[{"type":"code"}]},{"text":" for example, to see how it looks.","type":"text"}]},{"type":"paragraph","content":[{"text":"Y axis label - ","type":"text"},{"text":"ylab()","type":"text","marks":[{"type":"code"}]},{"text":". You can change this to what you like. You can also add an ","type":"text"},{"text":"ylab()","type":"text","marks":[{"type":"code"}]},{"text":" to change the X axis label.","type":"text"}]},{"type":"paragraph","content":[{"text":"Font size - ","type":"text"},{"text":"theme(text = element_text(size = 14))","type":"text","marks":[{"type":"code"}]},{"text":". There are many things you can change in the theme, such as text position, font type, plot margins, etc. See here: ","type":"text"},{"type":"inlineCard","attrs":{"url":"https://ggplot2.tidyverse.org/reference/element.html"}},{"text":" ","type":"text"}]},{"type":"paragraph","content":[{"text":"There are a multitude of other modifications you can make to a ggplot, too many to describe in this workshop. But there are plenty of online guides on how to modify ggplot plots. Here's an example: ","type":"text"},{"type":"inlineCard","attrs":{"url":"http://www.sthda.com/english/wiki/be-awesome-in-ggplot2-a-practical-guide-to-be-highly-effective-r-software-and-data-visualization"}},{"text":" ","type":"text"}]},{"type":"paragraph","content":[{"text":"The above changes and more can be applied to any ampvis2 or ggplot plot in this workshop.","type":"text"}]},{"type":"paragraph","content":[{"text":"Once you have your rarefecation plot looking how you like, you can export it as a 300dpi (i.e. publication quality) tiff or pdf file:","type":"text"}]},{"type":"heading","attrs":{"level":2},"content":[{"text":"Exporting your plot as a file","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"You can save your plot as a 300dpi (i.e. publication quality) tiff or pdf file. ","type":"text"},{"text":"These files can be found in your working directory.","type":"text","marks":[{"type":"strong"}]}]},{"type":"paragraph","content":[{"text":"Tip:","type":"text","marks":[{"type":"strong"}]},{"text":" you can adjust the width and height of the saved images by changing ","type":"text"},{"text":"width =","type":"text","marks":[{"type":"code"}]},{"text":" and ","type":"text"},{"text":"height =","type":"text","marks":[{"type":"code"}]},{"text":" in the code below (and every time ggsave appears in this workflow). ","type":"text"}]},{"type":"paragraph","content":[{"text":"Export as a 300dpi TIFF","type":"text"}]},{"type":"codeBlock","content":[{"text":"tiff_exp <- paste0(\"rarefaction_\", group, \"_samples_.tiff\")\nggsave(file = tiff_exp, dpi = 300, compression = \"lzw\", device = \"tiff\", plot = p, width = 20, height = 20, units = \"cm\")","type":"text"}]},{"type":"paragraph","content":[{"text":"Export as a pdf","type":"text"}]},{"type":"codeBlock","content":[{"text":"pdf_exp <- paste0(\"rarefaction_\", group, \"_samples.pdf\")\nggsave(file = pdf_exp, device = \"pdf\", plot = p, width = 20, height = 20, units = \"cm\")","type":"text"}]},{"type":"paragraph","content":[{"text":"You can now find these files in your working directory (which you originally defined in the 'Setting up your analysis environment' section).","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":1},"content":[{"text":"4d. Diversity index plots and statistics - single categorical variable","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"The overview section outlined (with links and references) the alpha diversity indices that can be examined.","type":"text"}]},{"type":"paragraph","content":[{"text":"Briefly, these are: Shannon’s index, Simpson's index, Chao1 richness and Observed ASVs. Each of these has strengths and weaknesses. It's up to you, the researcher, to explore the literature and decide which is the best index to use for your data.","type":"text"}]},{"type":"paragraph","content":[{"text":"This section will generate the diversity index scores and then plot them.","type":"text"}]},{"type":"paragraph","content":[{"text":"This is based on the variable you chose in the previous sections. If you’ve forgotten what variable group you’ve selected, type:","type":"text"}]},{"type":"codeBlock","content":[{"text":"group","type":"text"}]},{"type":"heading","attrs":{"level":2},"content":[{"text":"Calculate alpha diversity","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"First you need to calculate the alpha diversity index scores using ampvis2 function ","type":"text"},{"text":"amp_alphadiv()","type":"text","marks":[{"type":"code"}]},{"text":". This will calculate all 4 indices.","type":"text"}]},{"type":"codeBlock","content":[{"text":"div_indices <- amp_alphadiv(ampvisdata, richness = T)","type":"text"}]},{"type":"paragraph","content":[{"text":"The raw diversity index scores are now attached to your metadata table. You can view the raw scores by (scroll to last 5 columns for scores):","type":"text"}]},{"type":"codeBlock","content":[{"text":"div_indices","type":"text"}]},{"type":"paragraph","content":[{"text":"Export this table to your working directory as a csv file:","type":"text"}]},{"type":"codeBlock","content":[{"text":"write.csv(div_indices, \"diversity_indices_raw_scores.csv\")","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":2},"content":[{"text":"Choose the index you want to plot","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Choose the diversity index scores you want to plot.","type":"text"}]},{"type":"codeBlock","content":[{"text":"indicname <- \"Shannon\"","type":"text"}]},{"type":"paragraph","content":[{"text":"Shannon's index is used by default (","type":"text"},{"text":"\"Shannon\"","type":"text","marks":[{"type":"code"}]},{"text":") Change this to ","type":"text"},{"text":"\"Simpson\"","type":"text","marks":[{"type":"code"}]},{"text":" to plot Simpson's index scores, ","type":"text"},{"text":"\"Chao1\"","type":"text","marks":[{"type":"code"}]},{"text":" for Chao1 richness or ","type":"text"},{"text":"\"ObservedOTUs\"","type":"text","marks":[{"type":"code"}]},{"text":" for Observed ASVs.","type":"text"}]},{"type":"heading","attrs":{"level":2},"content":[{"text":"Box and whisker plot","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Create a box and whisker plot showing diversity per group. You can view the basic plot like so:","type":"text"}]},{"type":"codeBlock","content":[{"text":"p <- ggplot(div_indices, aes_string(x=group, y=indicname)) + \n geom_boxplot(fill = \"grey\", outlier.size = 2.5, lwd=0.7)\np","type":"text"}]},{"type":"paragraph","content":[{"text":"And you can use the above plot as-is or modify it. For example:","type":"text"}]},{"type":"codeBlock","content":[{"text":"p <- p + theme_bw() + \nylab(indicname) + \nxlab(\"\") + \ncoord_flip() + \ntheme(text = element_text(size = 17)) + \ntheme(legend.position = \"none\") + \ngeom_smooth(se = F, aes(group=1, color = \"red\"))\np","type":"text"}]},{"type":"paragraph","content":[{"text":"Change the above code to modify the plot as you like. See the rarefaction curve section for more details on modifying a plot.","type":"text"}]},{"type":"paragraph","content":[{"text":"Briefly, the above code changes the theme (","type":"text"},{"text":"theme_bw()","type":"text","marks":[{"type":"code"}]},{"text":"), the Y axis label (","type":"text"},{"text":"ylab(indicname)","type":"text","marks":[{"type":"code"}]},{"text":" - note: the current y axis name is inside a variable called ","type":"text"},{"text":"indicname","type":"text","marks":[{"type":"code"}]},{"text":", so you need to replace it with text in quotes, e.g. ","type":"text"},{"text":"\"Chao richness\"","type":"text","marks":[{"type":"code"}]},{"text":"), removes the X axis label (","type":"text"},{"text":"xlab(\"\")","type":"text","marks":[{"type":"code"}]},{"text":"), flips the X and Y axis (","type":"text"},{"text":"coord_flip()","type":"text","marks":[{"type":"code"}]},{"text":"), changes the font size to 17 and removes the legend (both in ","type":"text"},{"text":"theme()","type":"text","marks":[{"type":"code"}]},{"text":"). Finally, a trend line is added (","type":"text"},{"text":"geom_smooth()","type":"text","marks":[{"type":"code"}]},{"text":"). Any of these components can be removed or modified as you like.","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Save your plot as a tiff and pdf file.","type":"text"}]},{"type":"codeBlock","content":[{"text":"tiff_exp <- paste0(\"alpha_div_box_plot_\", group, \"_\", indicname, \"_samples.tiff\")\nggsave(file = tiff_exp, dpi = 300, compression = \"lzw\", device = \"tiff\", plot = p, width = 20, height = 20, units = \"cm\")\npdf_exp <- paste0(\"alpha_div_box_plot_\", group, \"_\", indicname,\"_samples.pdf\")\nggsave(file = pdf_exp, device = \"pdf\", plot = p, width = 20, height = 20, units = \"cm\")","type":"text"}]},{"type":"heading","attrs":{"level":2},"content":[{"text":"Statistical analysis","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"To compare the differences between groups within your variable, a ","type":"text"},{"text":"Kruskal-Wallis test","type":"text","marks":[{"type":"link","attrs":{"href":"https://www.technologynetworks.com/informatics/articles/the-kruskal-wallis-test-370025"}}]},{"text":" (one-way analysis of variance) can be performed to test for overall differences between all groups, and a ","type":"text"},{"text":"Wilcoxon rank sum","type":"text","marks":[{"type":"link","attrs":{"href":"http://www.sthda.com/english/wiki/unpaired-two-samples-wilcoxon-test-in-r"}}]},{"text":" to test pairwise differences between each group.","type":"text"}]},{"type":"paragraph","content":[{"text":"Kruskal-Wallis:","type":"text","marks":[{"type":"strong"}]}]},{"type":"codeBlock","content":[{"text":"kr_all <- kruskal.test(get(indicname) ~ get(group), data = div_indices)","type":"text"}]},{"type":"paragraph","content":[{"text":"This runs the test on the diversity index you chose: ","type":"text"},{"text":"indicname","type":"text","marks":[{"type":"code"}]},{"text":", the variable: ","type":"text"},{"text":"group","type":"text","marks":[{"type":"code"}]},{"text":" and uses the diversity indices table you generated: ","type":"text"},{"text":"div_indices","type":"text","marks":[{"type":"code"}]}]},{"type":"paragraph","content":[{"text":"To remind you what diversity index and variable you're examining:","type":"text"}]},{"type":"codeBlock","content":[{"text":"paste0(\"Diverity index: \", indicname)\npaste0(\"Variable: \", group)","type":"text"}]},{"type":"paragraph","content":[{"text":"To see the results:","type":"text"}]},{"type":"codeBlock","content":[{"text":"kr_all","type":"text"}]},{"type":"paragraph","content":[{"text":"Pairwise Wilcoxon rank sum test:","type":"text","marks":[{"type":"strong"}]}]},{"type":"codeBlock","content":[{"text":"wt_pair <- pairwise.wilcox.test(div_indices[[indicname]], div_indices[[group]], p.adjust.method = \"BH\")\nwt_pair <- data.frame(wt_pair$p.value)\ncolnames(wt_pair) <- gsub(\"X\", \"\", colnames(wt_pair))","type":"text"}]},{"type":"paragraph","content":[{"text":"To see the pairwise results (p values).","type":"text"}]},{"type":"codeBlock","content":[{"text":"wt_pair","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Go to the next section: Beta Diversity","type":"text"}]},{"type":"layoutSection","content":[{"type":"layoutColumn","attrs":{"width":50.0},"content":[{"type":"paragraph","content":[{"text":"Previous section: Setting up your R environment","type":"text","marks":[{"type":"link","attrs":{"href":"https://eresearchqut.atlassian.net/wiki/spaces/EG/pages/2593129682"}}]}]}]},{"type":"layoutColumn","attrs":{"width":50.0},"content":[{"type":"paragraph","content":[{"text":"Next section: Beta diversity","type":"text","marks":[{"type":"link","attrs":{"href":"https://eresearchqut.atlassian.net/wiki/spaces/EG/pages/2593129764"}}]}]}]}]},{"type":"paragraph"},{"type":"paragraph"},{"type":"paragraph"}],"version":1}

Browser not supported