{"type":"doc","content":[{"type":"paragraph","content":[{"text":"Overview","type":"text"}]},{"type":"extension","attrs":{"layout":"default","extensionType":"com.atlassian.confluence.macro.core","extensionKey":"toc","parameters":{"macroParams":{"style":{"value":"none"}},"macroMetadata":{"macroId":{"value":"323e193e-8832-4822-83cb-986c20a34802"},"schemaVersion":{"value":"1"},"title":"Table of Contents"}},"localId":"83d3b029-25aa-469e-a421-6c5981425095"}},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"In this section we’re going to:","type":"text"}]},{"type":"bulletList","content":[{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Install and load the R packages we need to run the analysis","type":"text"}]}]},{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"Import our taxonomic abundance table into R","type":"text"}]}]},{"type":"listItem","content":[{"type":"paragraph","content":[{"text":"View a summary of our abundance table","type":"text"}]}]}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"We’re going to be running various commands in R. To do this, copy and paste the code into the R script you created, highlight all the code you want to run, then press the run button:","type":"text"}]},{"type":"mediaSingle","attrs":{"layout":"align-start","width":122,"widthType":"pixel"},"content":[{"type":"media","attrs":{"width":71,"id":"1b1f54f3-01b5-45db-a1fd-78bb143e5993","collection":"contentId-2236678156","type":"file","height":26}}]},{"type":"heading","attrs":{"level":1},"content":[{"text":"Install required R packages","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Copy and paste the following code into the R script you just created, then run the code. This will install all the required packages and dependencies and may take 45 minutes or more to complete. It may prompt you occasionally to update packages - select 'a' for all if/when this occurs.","type":"text"}]},{"type":"paragraph","content":[{"text":"NOTE: you only need to run this section once if you’re running this analysis on your own laptop/PC, and you don’t need to run it if you’re using an rVDI machine as all the packages are already installed.","type":"text","marks":[{"type":"underline"},{"type":"textColor","attrs":{"color":"#ff5630"}},{"type":"strong"}]}]},{"type":"paragraph"},{"type":"codeBlock","content":[{"text":"#### Metagenomics analysis ####\n\n# When you see '## USER INPUT', this means you have to modify the code for your computer or dataset. All other code can be run as-is (i.e. you don't need to understand the code, just run it)\n\n#### 1. Installing required packages ####\n\n# **NOTE: this section only needs to be run once (or occasionally to update the packages)\n# Install devtools\ninstall.packages(\"devtools\", repos = \"http://cran.us.r-project.org\")\n# Install R packages. This only needs to be run once.\n\n# Make a vector of CRAN and Bioconductor packages\nbioconductor_packages <- c(\"VariantAnnotation\", \"biomaRt\", \"clusterProfiler\", \"org.Hs.eg.db\")\ncran_packages <- c(\"devtools\", \"tidyverse\", \"DT\", \"gt\", \"openxlsx\", \"dplyr\", \"scales\", \"ggplot2\", \"plotly\", \"tidyr\", \"ggsci\", \"viridis\", \"vcfR\", \"data.table\", \"remotes\", \"vegan\")\n\n# Compares installed packages to above packages and returns a vector of missing packages\nnew_packages <- bioconductor_packages[!(bioconductor_packages %in% installed.packages()[,\"Package\"])]\nnew_cran_packages <- cran_packages[!(cran_packages %in% installed.packages()[,\"Package\"])]\n\n# Install missing bioconductor packages\nif (!requireNamespace(\"BiocManager\", quietly = TRUE))\n install.packages(\"BiocManager\")\nBiocManager::install(new_packages)\n\n# Install missing cran packages\nif (length(new_cran_packages)) install.packages(new_cran_packages, repos = \"http://cran.us.r-project.org\")\n\n# Update all installed packages to the latest version\nupdate.packages(bioconductor_packages, ask = FALSE)\nupdate.packages(cran_packages, ask = FALSE, repos = \"http://cran.us.r-project.org\")\n\n# Install ampvis2 (needs to be installed from Github)\nremotes::install_github(\"kasperskytte/ampvis2\")\n","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":1},"content":[{"text":"Load required R packages","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"This section loads the packages you’ve installed in the previous section. Unlike installing packages, this needs to be run every time and should only take a few seconds to run.","type":"text"}]},{"type":"paragraph"},{"type":"codeBlock","content":[{"text":"#### 2. Loading required packages ####\n\n# This section needs to be run every time\n# Load packages\nbioconductor_packages <- c(\"VariantAnnotation\", \"biomaRt\", \"clusterProfiler\", \"org.Hs.eg.db\")\ncran_packages <- c(\"devtools\", \"tidyverse\", \"DT\", \"gt\", \"openxlsx\", \"dplyr\", \"scales\", \"ggplot2\", \"plotly\", \"tidyr\", \"ggsci\", \"viridis\", \"vcfR\", \"data.table\", \"remotes\", \"vegan\")\nlapply(cran_packages, require, character.only = TRUE)\nlapply(bioconductor_packages, require, character.only = TRUE)\nlibrary(ampvis2)\n","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":1},"content":[{"text":"Set your working directory","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"‘Working directory’ is an important concept in R. It defines where R automatically looks for data files and where it outputs results (tables, figures, etc).","type":"text"}]},{"type":"paragraph","content":[{"text":"To set your working directory, click ‘","type":"text"},{"text":"Session","type":"text","marks":[{"type":"strong"}]},{"text":"’ → ‘","type":"text"},{"text":"Set working directory","type":"text","marks":[{"type":"strong"}]},{"text":"’ → ‘","type":"text"},{"text":"Choose working directory","type":"text","marks":[{"type":"strong"}]},{"text":"’ and then choose the ","type":"text"},{"text":"H:/meta_workshop/R_analysis","type":"text","marks":[{"type":"textColor","attrs":{"color":"#36b37e"}},{"type":"strong"}]},{"text":" directory.","type":"text"}]},{"type":"paragraph"},{"type":"heading","attrs":{"level":1},"content":[{"text":"Summary of DADA2 quality filtration","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"DADA2 ","type":"text","marks":[{"type":"link","attrs":{"href":"https://benjjneb.github.io/dada2/"}}]},{"text":"is the main tool used by the nfcore/ampliseq workflow to examine taxonomy in amplicon datasets. It initially completes a ","type":"text"},{"text":"series of quality filtration steps","type":"text","marks":[{"type":"link","attrs":{"href":"https://benjjneb.github.io/dada2/tutorial.html"}}]},{"text":", including quality filtration, denoising and chimeric read removal. ","type":"text"}]},{"type":"paragraph","content":[{"text":"As an initial task, we’re going to generate a table and figure that shows the percentage of reads were kept and removed in each filtration step. Percentage columns are the filtered, merged, and non-chimeric reads compared to the number of original, unfiltered reads. Thus the number of non-chimeric reads and percentage is the number and % of reads remaining after all filtration steps.","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Copy, paste and run the following into your RStudio script:","type":"text"}]},{"type":"codeBlock","content":[{"text":"#### 2. Summary of DADA2 quality filtration ####\n\n## Make a summary table\n\n# Import the DADA2 filtration summary information\ndada <- read.table(\"../illumina/results/dada2/dada2_stats.tsv\", sep = \"\\t\", header = T)\n\n# Calculate the filtration percentages and add as new columns\ndada$percentage.of.input.passed.filter <- paste0(round((dada$filtered/dada$DADA2_input)*100,2), \"%\")\ndada$percentage.of.input.denoised <- paste0(round((dada$denoised/dada$DADA2_input)*100,2), \"%\")\ndada$percentage.of.input.non.chimeric <- paste0(round((dada$nonchim/dada$DADA2_input)*100,2), \"%\")\n\n# Rearrange the columns and change the column names\ndada <- dada[c(1,2,3,6,4,7,5,8)]\ncolnames(dada) <- c(\"Sample\", \"Unfiltered_reads\", \"Filtered\", \"%\", \"Denoised\", \"%\", \"Non-chimeric\", \"%\")\n\n# Export as a csv file\nwrite.csv(dada, \"DADA2_filtration.csv\")\n","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"The final line in the above code box will write out the results as a table called ","type":"text"},{"text":"DADA2_filtration.csv","type":"text","marks":[{"type":"code"}]},{"text":". It will write this out in your working directory (","type":"text"},{"text":"H:/meta_workshop/R_analysis","type":"text","marks":[{"type":"textColor","attrs":{"color":"#36b37e"}},{"type":"strong"}]},{"text":"). You can go there in Windows File Explorer and open the table in Excel.","type":"text"}]},{"type":"heading","attrs":{"level":2},"content":[{"text":"Stacked bar plot of DADA2 filtration","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Now we'll plot the DADA2 filtration results (filtered vs unfiltered reads) as a stacked barplot.","type":"text"}]},{"type":"paragraph","content":[{"text":"We'll do this using the ","type":"text"},{"text":"ggplot ","type":"text","marks":[{"type":"link","attrs":{"href":"https://ggplot2.tidyverse.org/"}}]},{"text":"package.","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Run the following in your R script:","type":"text"}]},{"type":"codeBlock","content":[{"text":"## Plot the summary results\n\n# Pull out the required columns\ndada_bp <- data.frame(dada$Sample, dada$Unfiltered_reads - dada$`Non-chimeric`, dada$`Non-chimeric`)\ncolnames(dada_bp) <- c(\"Sample\", \"Removed reads\", \"Remaining reads\")\ndf <- tidyr::pivot_longer(dada_bp, cols =c(\"Removed reads\", \"Remaining reads\"))\ncolnames(df) <- c(\"Sample\", \"Reads\", \"Read_count\")\n\n# Plot the DADA2 filtration results\np <- ggplot(df, aes(Sample, Read_count, fill=Reads)) + geom_bar(stat = \"identity\", position = 'stack')\np <- p + scale_fill_brewer(palette=\"Dark2\") + ylab(\"Read count\") + theme_bw() +\n theme(text = element_text(size = 14), axis.text.x = element_text(angle = 90, size = 10))\np\n","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"The above will display the plot in RStudio.","type":"text"}]},{"type":"paragraph","content":[{"text":"You can also 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. ","type":"text"}]},{"type":"paragraph"},{"type":"codeBlock","content":[{"text":"# Export as a 300dpi tiff\ntiff_exp <- \"DADA2_summary.tiff\"\nggsave(file = tiff_exp, dpi = 300, compression = \"lzw\", device = \"tiff\", plot = p, width = 20, height = 20, units = \"cm\")\n# Export as a pdf\npdf_exp <- \"DADA2_summary.pdf\"\nggsave(file = pdf_exp, device = \"pdf\", plot = p, width = 20, height = 20, units = \"cm\")\n","type":"text"}]},{"type":"paragraph"},{"type":"paragraph","content":[{"text":"Go to the next section: Alpha diversity","type":"text"}]}],"version":1}

Browser not supported