PipeRline

General example

Introduction

This metabarcoding pipeline is based around the targets package, which is a Make-like pipeline tool for R. The benefit of this is that all the code is automatically run, and the pipeline skips costly runtime for tasks that are already up to date.

This page lists the general workflow and all the options that can be adapted for your dataset.

Clone the pipeRline github repository

The first step is to clone this github repository, which contains the required code and directory structure to run the pipeline. To do this, you will need Git installed on your computer. If you are running the pipeline in Rstudio, it is best to create a new project from the github repository. This can be done by going file > new project > version control > git then adding https://github.com/alexpiper/piperline.git then change the name of the project to whatever you are working with.

Alternatively, the repository can be cloned using the command line, change the ‘folder-name’ to the desired name

# Change into the main directory you wish to make the project in
cd metabarcoding

# Clone the repository
git clone https://github.com/alexpiper/piperline.git folder-name

Updating the pipeline

To update to the latest version of the pipeline, run the below code in the terminal.

git pull

Prior to analysis

The PipeRline workflow assumes that your sequencing data meets certain criteria:

1. Samples have been demultiplexed, i.e. split into individual per-sample fastq files. If you want the pipeline to calculate the index-switching rate, the fastq files need to be re-demultiplexed as the miseq does not put indexes in fasta headers by default
2. For paired-end sequencing data, the forward and reverse reads are in separate files (ie. not interleaved) with reads arranged in matched order.

The output directory should be unique for each sequencing run, named as the flowcell id, within a directory called data

For example:

root/
  ├── data/
     ├── CJL7D/
     

You will also need BLAST+ installed on your computer, it can be downloaded from the NCBI website

Sequencing reads

For this workflow to run, we will need some sequencing runs to work with. If you are working with MiSeq data, it is recommended that the data is demultiplexed again using bcl2fastq (see Demultiplex MiSeq run tab below), as the miseq does not put indexes in fasta headers by default which is required for the index swtiching calculation.

If you do not yet have any data, some test sequencing reads have been hosted on Zenodo for testing purposes. The below code in the Download example reads tab below will download these and put them inside the data folder.

Demultiplex MiSeq run

The below code is written for the Agriculture Victoria BASC server, and the locations will be different if you are using a different HPC cluster.

#load module
module load bcl2fastq2/2.20.0-foss-2018b

#raise amount of available file handles
ulimit -n 4000

###Run1

#Set up input and outputs
inputdir=/group/sequencing/210219_M03633_0489_000000000-JDYG3 #CHANGE TO YOUR SEQ RUN
outputdir=/group/pathogens/IAWS/Projects/Metabarcoding/dros_surveillance/data/JDYG3 #CHANGE TO YOUR DATA FOLDER RUN
samplesheet=/group/pathogens/IAWS/Projects/Metabarcoding/dros_surveillance/SampleSheet_JDYG3.csv #CHANGE TO YOUR SAMPLESHEET

# convert samplesheet to unix format
dos2unix $samplesheet

#Demultiplex
bcl2fastq -p 12 --runfolder-dir $inputdir \
--output-dir $outputdir \
--sample-sheet $samplesheet \
--no-lane-splitting --barcode-mismatches 1

# Copy other necessary files and move fastqs
cd $outputdir
cp -r $inputdir/InterOp $outputdir
cp $inputdir/RunInfo.xml $outputdir
cp $inputdir/[Rr]unParameters.xml $outputdir
cp $samplesheet $outputdir
mv **/*.fastq.gz $outputdir

# Append fcid to start of sample names if missing
fcid=$(echo $inputdir | sed 's/^.*-//')
for i in *.fastq.gz; do
  if ! [[ $i == $fcid* ]]; then
  new=$(echo ${fcid} ${i}) #append together
  new=$(echo ${new// /_}) #remove any white space
  mv -v "$i" "$new"
  fi
done

Optional: Run R on BASC

You may wish to run this workflow through the command line on Agriculture Victorias BASC server in order to take advantage of more processing power. To do this, you can start a new SLURM interactive session. Press the CODE button to the lower right to display the code for this optional step.

# Create new interactive SLURM session
sinteractive --ntasks=1 --cpus-per-task=10 --mem-per-cpu=10GB --time=72:00:00

module load R/4.2.0-foss-2021b
module load pkgconfig/1.5.1-GCCcore-9.3.0-Python-3.8.2
module load GDAL/3.3.0-foss-2021a
module load BLAST+/2.11.0-gompi-2020a
module load Pandoc/2.5
module load ZeroMQ/4.3.2-GCCcore-9.3.0

# Load R
R

Download example reads

Example tephritid metabarcoding reads have been hosted on Zenodo, either download them manually or with R below:

# Source the functions file which contains code to download from zenodo
source("R/functions.R")

# Create directory for data
if(!dir.exists("data/K77JP")) {dir.create("data/K77JP", recursive = TRUE)}
if(!dir.exists("data/K77JP/InterOp")) {dir.create("data/K77JP/InterOp", recursive = TRUE)}

# Download files from zenodo
download_zenodo(
  doi = "10.5281/zenodo.7112162",
  path = "data/K77JP"
)

# Move the interop files to the interop folder
fs::dir_ls(path="data/K77JP", glob="*.bin") %>%
  purrr::map(function(x){
    fs::file_copy(path = x, new_path = x %>% str_replace("data/K77JP", "data/K77JP/InterOp"))
    file.remove(x)
  })

Install and load R packages and setup directories

This pipeline depends on a number of other R packages which this step will install. The versions of these packages are managed using renv to ensure they match the versions the pipeline was developed on. The first time installing packages may take some time, but they should be quick to load for any future runs.

# Load the targets and renv packages
library(targets)
library(tarchetypes)

# Load all packages using renv::restore() and source package list
renv::restore()
source("_targets_packages.R")

# Source ancillary functions
source("R/functions.R")
source("R/themes.R")

Reference databases

The pipeRline workflow uses the IDTAXA and BLAST software to assign taxonomy to the sequence reads. These require a pre-trained IDTAXA model, and a fasta file with heirarchial taxonomy. For insects, these

Download reference database

Reference databases for insects and arachnids have been hosted on Zenodo, either download them manually, or with R below:

# Source the functions file which contains code to download from zenodo
source("R/functions.R")

# Download files from zenodo
download_zenodo(
  doi = "10.5281/zenodo.7655352",
  path = "reference"
)

Create your own reference database

If you need a new database, these can be created from genbank and blast using taxreturn or other simialr database curation packages.

>Accession|Root;Kingdom;Phylum;Class;Order;Family;Genus;Species
ACCTAGAAAGTCGTAGATCGAAGTTGAAGCATCGCCCGATGATCGTCTGAAGCTGTAGCATGAGTCGATTTTCACATTCAGGGATACCATAGGATAC
>Root;Kingdom;Phylum;Class;Order;Family;Genus;Species
CGCTAGAAAGTCGTAGAAGGCTCGGAGGTTTGAAGCATCGCCCGATGGGATCTCGTTGCTGTAGCATGAGTACGGACATTCAGGGATCATAGGATAC

The IDTAXA model can then be trained on this reference fasta - See the IDTAXA manual

Create sample tracking sheet

In order to track samples and relevant QC statistics throughout the metabarcoding pipeline, we will first create a new sample tracking sheet from our input illumina samplesheets. This function requires both the SampleSheet.csv used for the sequencing run, and the runParameters.xml, both of which should have been automatically obtained from the demultiplexed sequencing run folder in the bash step above

runs <- dir("data/") #Find all directories within data
SampleSheet <- list.files(paste0("data/", runs), pattern= "SampleSheet", full.names = TRUE)
runParameters <- list.files(paste0("data/", runs), pattern= "[Rr]unParameters.xml", full.names = TRUE)

# Create samplesheet containing samples and run parameters for all runs
samdf <- create_samplesheet(SampleSheet = SampleSheet, runParameters = runParameters, template = "V4") %>%
  distinct()

# Check that sample_ids contain fcid, if not; attatch
samdf <- samdf %>%
  mutate(sample_id = case_when(
    !str_detect(sample_id, fcid) ~ paste0(fcid,"_",sample_id),
    TRUE ~ sample_id
  ))

# Check that samples match samplesheet
fastqFs <- purrr::map(list.dirs("data", recursive=FALSE),
                      list.files, pattern="_R1_", full.names = TRUE) %>%
  unlist() %>%
  str_remove(pattern = "^(.*)\\/") %>%
  str_remove(pattern = "(?:.(?!_S))+$")

# Filter undetermined reads from sample sheet
fastqFs <- fastqFs[!str_detect(fastqFs, "Undetermined")]

# Check for fastq files that are missing from samplesheet
if (length(setdiff(fastqFs, samdf$sample_id)) > 0) {warning("The fastq file/s: ", setdiff(fastqFs, samdf$sample_id), " are not in the sample sheet") }

# Check for sample_ids that dont have a corresponding fastq file
if (length(setdiff(samdf$sample_id, fastqFs)) > 0) {
  warning(paste0("The fastq file: ",
                 setdiff(samdf$sample_id, fastqFs),
                 " is missing, dropping from samplesheet \n")) 
  samdf <- samdf %>%
    filter(!sample_id %in% setdiff(samdf$sample_id, fastqFs))
}

# Write out sample tracking sheet
write_csv(samdf, "sample_data/Sample_info.csv")

Add PCR primers to sample sheet

This can either be done manually by editing the sample_data/Sample_info.csv file, or it can be done in R as below.

If a single primer set was used across all samples, these can simply be added using a mutate call (tab 1).

If different primers were used for different samples, these can be set using pattern matching on the sample names with case_when (tab 2).

If multiple primer sets are used per sample, before indexing, the pipeline will conduct an extra round of demultiplexing to split each sample by primer. This option can be set by splitting each primer set with a semicolon (tab 3)

Single primer set

# Add primers to sample sheet
samdf <- samdf %>%
  mutate(pcr_primers = "fwhF2-fwhR2n",
  for_primer_seq = "GGDACWGGWTGAACWGTWTAYCCHCC",
  rev_primer_seq = "GTRATWGCHCCDGCTARWACWGG"
  )

write_csv(samdf, "sample_data/Sample_info.csv")

Different primer sets per sample

# Add PCR primers to sample sheet
samdf <- samdf %>%
  mutate(pcr_primers = case_when(
    str_detect(sample_name, "primer1") ~ "fwhF2-fwhR2n",
    str_detect(sample_name, "primer2") ~ "EIF3LminiF4-EIF3lminiR4"
  ),
  for_primer_seq = case_when(
    str_detect(sample_name, "primer1") ~ "GGDACWGGWTGAACWGTWTAYCCHCC",
    str_detect(sample_name, "primer2") ~ "GATGCGYCGTTATGCYGATGC"
  ),
  rev_primer_seq = case_when(
    str_detect(sample_name, "primer1") ~ "GTRATWGCHCCDGCTARWACWGG",
    str_detect(sample_name, "primer2") ~ "TTRAAYACTTCYARATCRCC"
  ))

write_csv(samdf, "sample_data/Sample_info.csv")

Multiple primer sets per sample

# Add PCR primers to sample sheet
samdf <- samdf %>%
  mutate(
    pcr_primers = "fwhF2-fwhR2nDac;EIF3LminiF4-EIF3lminiR4",
    for_primer_seq = "GGDACWGGWTGAACWGTWTAYCCHCC;GATGCGYCGTTATGCYGATGC",
    rev_primer_seq = "GTRATWGCHCCIGCTAADACHGG;TTRAAYACTTCYARATCRCC"
    )

write_csv(samdf, "sample_data/Sample_info.csv")

Create parameters file

The parameters file table the respective target gene, reference databases, and filtering parameters for each primer set used to amplify the samples. As the pcr_primers column is used in the pipeline to match the parameters to the respective sample, it is critical that the primer names match those set in the previous step.

If a single primer set was used across all samples, these can simply be added to the table as below (tab 1).

If different primers were used for different samples, these can be set as different rows in the table using c() when creatign the tibble (tab 2).

If multiple primer reference databases are to be used for each sample, these can be set by splitting each a semicolon. The taxonomic assignment will be conducted sequentially through the databases from left to right, with the second reference database only being used for those ASVs that couldnt be assigned to species level using the first databases (tab 3)

Parameter options

• Primer parameters:
  • pcr_primers - Name of PCR primers - must match samdf file
  • target_gene - Name of target gene
  • max_primer_mismatch - How much mismatch to allow when detecting primer sequences
• Read filtering:
  • read_min_length - Minimum length of primer trimmed reads
  • read_max_length - Maximum length of primer trimmed reads
  • read_max_ee - Maximum expected errors of primer trimmed reads
  • read_trunc_length - Length to cut all longer reads to
  • read_trim_left - Remove this many bp from left side of primer trimmed reads
  • read_trim_right - Remove this many bp from right side of primer trimmed reads
• ASV filtering
  • asv_min_length - Minimum length of amplicon
  • asv_max_length - Maximum length of amplicon
  • high_sensitivity - Option to turn on pseudo pooling, setting this to FALSE greatly speeds up ASV inference at the expense of detecting fewer rare (<5 reads) taxa.
  • concat_unmerged - Retain any unmerged read pairs by concatenating them together separated by a string of 10 N bases
  • genetic_code - Genetic code for amplicon - see Biostrings::GENETIC_CODE_TABLE
  • coding - Is the amplicon from a protein coding gene
  • phmm - Path to profile hidden markov model (Optional)
• Taxonomic assignment:
  • idtaxa_db - Path to trained IDTAXA model
  • ref_fasta - Path to fasta file of reference database
  • idtaxa_confidence - Minimum bootstrap confidence for IDTAXA
  • run_blast - Whether a blast top hit search should be conducted in addition to IDTAXA
  • blast_min_identity - Minimum nucleotide identity for BLAST
  • blast_min_coverage - Minimum query coverage for BLAST
  • target_kingdom - Subset to target kingdom
  • target_phylum - Subset to target phylum
  • target_class - Subset to target class
  • target_order - Subset to target order
  • target_family - Subset to target family
  • target_genus - Subset to target genus
  • target_species - Subset to target species
• Sample & Taxon filtering:
  • min_sample_reads - Minimum reads per sample after filtering
  • min_taxa_reads - Minimum reads per ASV to retain
  • min_taxa_ra - Minimum relative abundance per ASV to retain. 1e-4 is 0.01%
• General pipeline parameters
  • threads - Number of CPU threads to run the pipeline across

Single primer set

# Params to add in step_add_parameters
params <- tibble(
  # Primer parameters
  pcr_primers = "fwhF2-fwhR2n",
  target_gene="COI",
  max_primer_mismatch=0,

  # Read filtering
  read_min_length = 20,
  read_max_length = Inf,
  read_max_ee = 1,
  read_trunc_length = 150,
  read_trim_left = 0, 
  read_trim_right = 0,
  
  # ASV filtering
  asv_min_length = 195, 
  asv_max_length = 215,
  high_sensitivity = TRUE,
  concat_unmerged = FALSE,
  genetic_code = "SGC4",
  coding = TRUE,
  phmm = "reference/folmer_fullength_model.rds",
  
  # Taxonomic assignment
  idtaxa_db = "reference/idtaxa_bftrimmed.rds",
  ref_fasta = "reference/insecta_hierarchial_bftrimmed.fa.gz",
  idtaxa_confidence = 60,
  run_blast=TRUE,
  blast_min_identity = 97,
  blast_min_coverage = 90,
  target_kingdom = "Metazoa",
  target_phylum = "Arthropoda",
  target_class = NA,
  target_order = NA,
  target_family = NA,
  target_genus = NA,
  target_species= NA,
  
  # Sample & Taxon filtering
  min_sample_reads = 1000,
  min_taxa_reads= NA,
  min_taxa_ra = 1e-4, #1e-4 is 0.01%
    
  # General pipeline parameters
  threads = 1
)

write_csv(params, "sample_data/loci_params.csv")

Different or multiple primer sets per sample

# Params to add in step_add_parameters
params <- tibble(
  # Primer parameters
  pcr_primers = c("fwhF2-fwhR2nDac", "EIF3LminiF4-EIF3lminiR4"),
  target_gene=c("COI", "EIF3L"),
  max_primer_mismatch=1,

  # Read filtering
  read_min_length = 20,
  read_max_length = Inf,
  read_max_ee = 1,
  read_trunc_length = 150,
  read_trim_left = 0,
  read_trim_right = 0,
  
  # ASV filtering
  asv_min_length = c(195, 207),
  asv_max_length = c(215, 227),
  high_sensitivity = TRUE,
  concat_unmerged = FALSE,
  genetic_code = c("SGC4", "SGC0"),
  coding = c(TRUE, TRUE),
  phmm = c("reference/phmm/Bactrocera_COI.rds", "reference/phmm/Bactrocera_EIF3L.rds"),
  
  # Taxonomic assignment
  idtaxa_db = c("reference/COI_idtaxa.rds","reference/EIF3L_idtaxa.rds"),
  ref_fasta = c("reference/COI_hierarchial.fa.gz", "reference/EIF3L_hierarchial.fa.gz"),
  idtaxa_confidence = 60,
  run_blast=TRUE,
  blast_min_identity = 97,
  blast_min_coverage = 90,
  target_kingdom = c("Metazoa", "Metazoa"),
  target_phylum = c("Arthropoda","Arthropoda"),
  target_class = NA,
  target_order = NA,
  target_family = NA, 
  target_genus = NA,  
  target_species= NA,  
  
  # Sample & Taxon filtering
  min_sample_reads = c(1000, 1000),
  min_taxa_reads= NA, 
  min_taxa_ra = c(1e-4, 1e-4),
)

write_csv(params, "sample_data/loci_params.csv")

Multiple reference databases

params <- tibble(
  # Primer parameters
  pcr_primers = c("fwhF2-fwhR2nDac", "EIF3LminiF4-EIF3lminiR4"),
  target_gene=c("COI", "EIF3L"),
  max_primer_mismatch=0,

  # Read filtering
  read_min_length = 20,
  read_max_length = Inf,
  read_max_ee = 1,
  read_trunc_length = 150,
  read_trim_left = 0,
  read_trim_right = 0,
  
  # ASV filtering
  asv_min_length = c(195, 207),
  asv_max_length = c(215, 227),
  high_sensitivity = TRUE,
  concat_unmerged = FALSE,
  genetic_code = c("SGC4", "SGC0"),
  coding = c(TRUE, TRUE),
  phmm = c("reference/phmm/Bactrocera_COI.rds", "reference/phmm/Bactrocera_EIF3L.rds"),
  
  # Taxonomic assignment
  idtaxa_db = c("reference/COI_internal_idtaxa.rds;reference/COI_idtaxa.rds",
             "reference/EIF3L_internal_idtaxa.rds;reference/EIF3L_idtaxa.rds"),
  ref_fasta = c("reference/COI_internal.fa.gz;reference/COI_hierarchial.fa.gz",
               "reference/EIF3L_internal.fa.gz;reference/EIF3L_hierarchial.fa.gz"),
  idtaxa_confidence = 60,
  run_blast=TRUE,
  blast_min_identity = 97,
  blast_min_coverage = 90,
  target_kingdom = c("Metazoa", "Metazoa"),
  target_phylum = c("Arthropoda","Arthropoda"),
  target_class = NA,
  target_order = NA,
  target_family = NA, 
  target_genus = NA,  
  target_species= NA,  
  
  # Sample & Taxon filtering
  min_sample_reads = c(1000, 1000),
  min_taxa_reads= NA,
  min_taxa_ra = c(1e-4, 1e-4),
)

write_csv(params, "sample_data/loci_params.csv")

Run pipeline

Now that the sample data sheet and parameters are defined, the pipeline steps can now be run automatically using the tar_make() command

tar_make()

Results

After the run, two directories should have been made within the output folder, output/results/unfiltered which contains the results before any taxonomic subsetting, sample minimum abundance, and ASV minimum abundance filtering, and output/results/filtered containing the results after these steps. Within each of these folders there should be a series of outputs:

• summary.csv - This is the main summary file which lists the abundance of each ASV, the taxonomy assigned to it, and the read counts for each sample it was detected in.

• ps.rds - This phyloseq object contains the same information as the summary file, but in a format that can be analysed further using the phyloseq R package, a useful tool for plotting and further analysis of metabarcoding and microbiome datasets.

• asvs.fasta - This contains the amplification sequence variants inferred from the data set

• seqtab.csv - This lists the abundance of each ASV and the sample it was detected in. This file is generally for development purposes only, and all information contained here is present in the summmary csv.

• taxtab.csv - This lists the heirarchial taxonomy assigned to each ASV. This file is generally for development purposes only, and all information contained here is present in the summmary csv.

• samdf.csv - This is the sample tracking sheet that lists the sequencing details for each sample. This file is generally for development purposes only, and all information contained here is present in the summmary csv.

• raw.csv - This is a large file containing all the above results in a single csv in a ‘long’ format appropriate for further statistical analysis.

Quality control

The pipeline also outputs some quality control plots that should be checked to ensure the pipeline has run correctly and the outputs are as desired

• Per-flowcell quality checks

  • output/logs/FCID/FCID_flowcell_qc.pdf - This set of plots summarises the raw outputs from the sequencing run. The first page shows the average imaging intensity for the 90% percentile of the data for each tile on the sequencing flow cell. The second page shows the number of clusters (reads) imaged by the seqeuncer, and the number of clusters passing the machines filtering. The third page shows the Qscore by cycle for all reads.

  • output/logs/FCID/FCID_index_switching.pdf - This is a heatmap of the amount of index switching per index set and flowcell used. On the diagonal of the plot is the indexes applied to the samples, while the other squares on the grid display instances where either the i5 or i7 index has switched causing this read to end up in the undetermined reads file. The text above the plot lists the total reads demultiplex from that flowcell, the switch rate (only one of the indexes switching causing the read to end up in the undetermined reads file) as a percentage, the cross-contamination rate (both indexes switching the read ending up in another sample) as a percentage, and the number of “other reads” which were those where no applied indexes were detected. The cross-contamination rate is calculated as the square of the index switching rate as per (costello et al), and should be a rough guide to determine an appropriate ASV minimum abundance threshold for filtering. NOTE: This plot is only created when unique dual indexes were used on the sample.

  • output/logs/FCID/FCID_prefilt_qualplots.pdf - This multiplot summarises the read quality scores for the forward and reverse reads, pre filtering. These plots can serve as an indicator for adjusting read filtering parameters of the pipeline. The top two plots summarise the mean quality score by position for the forward and reverse reads respectively. For illumina data, these will start to decrease towards the end of the reads, particularly for the reverse reads. The two bottom plots show the cumulative expected errors along the reads, with different coloured dots for different quantiles of the read data (i.e. lower 10% to upper 90%). The red lines show the potential max expecter error cutoffs, with all reads above the set maximum expected error filter beign removed (param read_max_ee). For example, if the lower 10% quantile of reads goes above the red line for maxEE=1, you would expect 10% of the reads to be removed if read_max_ee was set to 1. As more errors commonly occur in the ends of the read, if there is sufficient overlap between the reads to cover the length of the amplicon, the reads could instead be truncated to smaller lengths to reduce the amount of reads that cross this minimum expected error filter which would retain more data. Note that when truncating reads you should aim to have at least 20bp of overlap to facilitate forward and reverse read merging.

  • output/logs/FCID/FCID_postfilt_qualplots.pdf - This multiplot summarises the read quality scores for the forward and reverse reads, pre filtering. These plots list the same information as the pre-filtering quality plots, and should be checked to ensure that the filtering has been correctly applied. i.e., if reads were truncated or trimmed on the left or right (params read_trunc_length, read_trim_left, read_trim_right) they should be shorter then the pre-filter plots, and if a maximum expected error filter (param read_max_ee) was applied, none of the EE quantiles should be above the corresponding red line

  • output/logs/FCID/FCID_pcrprimers_errormodel.pdf - There is a separate PDF for the forward and reverse error models for each primer set and flow cell. Each plot summarises the error model used for dada2 that was inferred from the data for the respective flow cell. The error rates for each possible transition (A→C, A→G, etc) are shown in each subplot. The dot-points are the error rates observed in the data for each quality score. The black line shows the estimated error rates from the machine-learning algorithm. The red line shows the error rates that would be expected if the Q-scores assigned by the sequencer were completely accurate. This plot can be more subjective to interpret, but generally you want to see the estimated error rates (black line) decrease as the quality score (x axis) increases, in a relatively linear fashion. Minor humps in the black line are no problem, but a very jagged line may indicate that the error model did not fit correctly. Note: for NextSeq and NovaSeq data, the binned quality scores provided by the sequencer can cause large humps or dips in the error model fit. This is generally not a big problem, see this github issue for a deeper discussion on fitting error models for binned quality scores.

• Whole dataset quality checks:

  • output/logs/ASV_cleanup_summary.pdf - This plot lists the abundance, and number of unique sequence variance inferred from dada2, their lengths, and whether they were retained or filtered out during the ASV filtering step. Generally you should see a lot of unique sequences filtered out in this step, but the majority of the abundance should remain.

  • output/logs/ASV_cleanup_summary.csv - This table lists the individual sequence variants, their abudnance, length and whether they were retained or filtered out during the ASV filtering step. This table is useful for deeper investigation to ensure that real sequences were not removed in any of the ASV filtering steps.

  • output/logs/taxonomic_assignment_summary.pdf - This plot summarises the genetic distance of the ASVs to the reference database (from BLAST top hit identity), and what taxonomy was assigned to each. This plot gives a good overview of how well the analysed sequences are represented in the reference database, and how successful the taxonomic assignment was.

  • output/logs/taxonomic_assignment_summary.csv - This table lists the individual sequence variants and their assigned taxonomy and confidence scores with IDTAXA. The columns to the right also indicate the BLAST top hit and their distance to the closest blast hit. This table is useful for deeper investigation of how different sequences were assigned.

  • output/logs/read_tracker.pdf This plot summarises the number of reads retained through all steps of the pipeline, as well as how many reads were classified to different taxonomic ranks. This plot provides a useful overview of all pipeline steps to determine if too many reads are being lost at a certain step

  • output/logs/read_tracker.csv - Similar to the above plot, this table summarises the number of reads retained through all steps of the pipeline, as well as how many reads were classified to different taxonomic ranks.

  • output/logs/accumulation_curve.pdf - This plot shows the accumulation of ASVs within each sample as sequencing depth increases. The curves should start to plateau if adequate sequencing depth has been achieved to characterise the community. If the curves have not reached a plateau it may indicate that these samples should be sequenced again so that all taxa are captured in the data.