Universal Amplicon Pipeline Tutorial

Created by Jesse McNichol and Nathan LR Williams.

Step 0 • Overview

Welcome

Welcome to the universal amplicon pipeline tutorial and setup workflow. This workflow helps you go from demultiplexed amplicon FASTQ files to a ready-to-run Snakemake configuration package.

The pipeline was designed around the 515F/926R universal primer pair, but it can be used with other amplicon primers provided the primer section is filled out correctly. It supports 16S and 18S reads, internal standards, and the metadata needed to generate the final files.

Use the sections below in order, then download the generated setup package when you are finished. The latest tested workflow is available in the pipeline repository on GitHub.

Why use a universal SSU rRNA primer set?

Metabarcoding using SSU rRNA as a marker gene is a powerful technique for profiling biological communities. Metabarcoding studies can be separated roughly into three groups:

  1. Studies of the “microbiome”—Bacteria and Archaea.
  2. Studies of microbial Eukarya, including phytoplankton and heterotrophic protists.
  3. Studies of macroscopic Eukarya, such as animals and plants.

The first group tends to use 16S SSU rRNA as a marker gene, while the second and third groups use 18S. PCR primer design has followed this division, with many primer sets targeting one group while discriminating against another.

However, 16S and 18S are, in evolutionary terms, the same molecule. SSU rRNA primers can therefore target both 16S and 18S in a single PCR assay. This is not two primer sets: one forward and one reverse primer amplify 16S from Archaea, Bacteria, chloroplasts, and mitochondria alongside eukaryotic nuclear 18S in the same tube.

That this is possible is remarkable. It shows that regions of the SSU rRNA molecule have remained sufficiently conserved across approximately 3.5 billion years of evolutionary history to amplify organisms as different as bacteria, chloroplasts, protists, and jellyfish.

Two example sequence pools

The first dataset comes from DNA extracted from 1 L of whole Pacific Ocean seawater filtered on a 0.2 µm Sterivex filter. It contains material ranging from prokaryotes and microbial eukaryotes to fragments, larvae, or propagules from animals.

LJ-Gradients sequence pool trace with a large 16S peak at 539 base pairs and a smaller 18S peak at 707 base pairs
A large peak at 539 bp represents 16S from Bacteria, Archaea, plastids, and mitochondria. The smaller peak at 707 bp represents 18S from protists, metazoa, picoeukaryotes, and other eukaryotes.

The second dataset comes from the same region of the Pacific Ocean, but larger organisms were concentrated with a net tow. Most of the DNA therefore comes from animals, protists, and other microeukaryotes.

NM-Gradients sequence pool trace with a small 16S peak at 546 base pairs and a large 18S peak at 736 base pairs
The pattern is reversed: the 16S peak at 546 bp is small, while the 18S peak at 736 bp is large. Both datasets were amplified with the same primer set and procedure.

Why this pipeline is needed

One reason this approach has not been used more commonly is the mismatch between sequencing read lengths and amplicon lengths. The primers this pipeline was designed for bind at SSU rRNA coordinates 515 and 926 using E. coli as a reference. The resulting bacterial amplicon is approximately 411 bp, although its exact length varies among organisms. Archaeal amplicons are approximately the same length.

In Eukarya, the corresponding primer coordinates are 562 and 1150 in S. cerevisiae, producing an amplicon of approximately 588 bp. Although 18S length varies across evolutionary and ecological groups, eukaryotic amplicons are about 177 bp longer than prokaryotic amplicons on average.

Modern Illumina sequencing has a maximum paired-end read length of 300 bp, and low-quality bases—especially in the reverse read—must often be trimmed. The remaining forward and reverse reads usually overlap for 16S, but not for 18S. Most conventional amplicon pipelines require that overlap and consequently discard 18S reads from 515Y/926R data.

Based on an initial prototype by Mike Lee, Jesse McNichol refined a method for separating 16S and 18S reads in silico, and the pipeline grew from there. Yi-Chun Yeh, Melody Aleman, and Colette Fletcher-Hoppe have since made important contributions, many inspired by Jed Fuhrman’s ideas.

Step 1 • Getting started

Repository and software setup

Install the required software, then enter a name for your project directory in the repository section below. The clone command updates automatically and creates a new directory with that project name. Use a different name for each analysis so its configuration, logs, and results stay together.

Install the requirements

Install Git and Conda before continuing. The official Snakemake installation guide recommends a Conda-based installation. Administrator access is not normally required when Conda is installed in your user directory.

After Conda is available, create the small controller environment used to launch the workflow. The pipeline-specific environments are created automatically during the first run.
Install Snakemake 
conda create --name snakemake --channel conda-forge --channel bioconda snakemake
conda activate snakemake
snakemake --version

Clone the pipeline repository

Choose the name of the directory that will contain this analysis. The generated command clones the pipeline from the public pipeline repository into a new directory with that name and then moves into it.

Use letters, numbers, periods, underscores, or hyphens. The command updates automatically.
Copy and run this command from the parent directory where you want the new project directory to be created. HTTPS cloning does not require a GitHub SSH key and is recommended when you only need to download and run the pipeline.
Clone into your project directory 
git clone https://github.com/Nwilliams96/515FY-926R-snakemake-NW-edits.git AMT29-analysis
cd AMT29-analysis
Step 2 • Config

Config.yml

Fill these fields to generate the config/config.yml file used by the Snakemake workflow.

Keep this outside the project clone. If the box above is unchecked, the workflow builds the databases here; if checked, it reuses the completed databases here.
Keep this outside the project clone. Snakemake checks this shared directory and reuses every completed rule environment whose definition has not changed. It creates only missing or updated environments, so future project clones do not reinstall the same software.
Enter the directory containing the FASTQ files. A trailing slash is optional.
Add one row for every internal standard used. These exact names become the <name>_ng sample columns and are printed in the recovery-plot legend and correction-method labels.
Use letters, numbers, periods, underscores, or hyphens. The same ordered list is used in the sample columns, internal-standard table, correction plots, and report.

Internal standard information

Enter the copy number, genome length, and complete 16S sequence for every named standard.

BP, DR, and TT are prefilled from the common preset. If you rename one, its preset definition is cleared so it cannot be mistaken for the newly named standard.

DADA2 read-length controls

These lengths are used for paired 16S denoising and to trim each 18S read before R1 and R2 are concatenated. The presets are the pipeline's established values.

Number of forward-read bases retained after primer removal. Longer values preserve sequence but may retain a poor-quality tail; reads shorter than this are discarded.
Number of reverse-read bases retained. For paired 16S, the retained forward and reverse reads must still overlap sufficiently to merge.

16S DADA2 controls (paired reads)

Forward reads above this expected-error total are discarded. Raising it retains more lower-quality reads but can increase erroneous ASVs.
Reverse-read equivalent of the forward threshold. Reverse reads often decline in quality sooner, so inspect filtering and merging losses when adjusting it.
A read is cut at the first base with a Phred score at or below this value. Raising it trims sooner and may cause more reads to fail the selected lengths.
Minimum matching overlap required to merge denoised forward and reverse reads. Lowering it may rescue pairs but reduces merging confidence.
Independent infers each sample separately. Pseudo-pooling performs a second pass informed by ASVs seen in multiple samples and can improve sensitivity to rare recurring variants, with extra compute.
Consensus detects chimeras per sample and combines the evidence. Pooled uses all samples together. None is useful only as a diagnostic because biological chimeras will remain.
Potential chimera parents must be at least this many times more abundant than the tested sequence. Raising it is less aggressive and can retain more variants, including possible chimeras.
Number of reads used to train DADA2's error model. More reads may stabilize the model for large or diverse datasets but require more time and memory.

18S DADA2 controls (concatenated reads)

The 18S path concatenates trimmed R1 and R2 and then runs single-end DADA2, so it has no paired-read overlap setting.

Concatenated reads above this expected-error total are discarded. Raising it retains more reads but can increase erroneous ASVs.
The current value of 0 disables quality-based truncation inside DADA2 because the component reads were already length-trimmed before concatenation.
Pseudo-pooling can improve sensitivity to rare variants that recur across samples, but runs an additional inference pass.
Controls how concatenated 18S chimeras are identified. Disabling this is intended for diagnostics, not a final biological result.
Higher values make chimera classification less aggressive and may retain both genuine low-abundance variants and unwanted chimeras.
Number of concatenated reads used for error-model training. Larger values increase training work and may help very large datasets.
Step 3 • Sample file

Sample File

Enter the number of samples to create the table, then type directly or paste a block of cells copied from Excel, Google Sheets, or equivalent software.

The table updates automatically. A spreadsheet range can be pasted starting in any cell.
Optional. This list defines every metadata column besides sample and the internal-standard fields. Add, rename, or remove columns—including depth, latitude, and longitude—then apply the list. The chosen columns are preserved downstream and become filters in the HTML taxonomy report.

No additional columns configured.

Upload requires only sample and the configured internal-standard amount columns. Optional uploaded columns are retained, and omitted optional columns are left blank.

Step 4 • Amplicon molarities

Amplicon Molarities

The 515F/926R primer pair targets both 16S and nuclear 18S rRNA. Unfortunately, there is a bias against 18S sequences when 16S and 18S sequences are mixed. This is because the 16S and 18S amplicons have different lengths (see the figures in the preamble) and are therefore sequenced at different rates, with a bias toward the shorter amplicons.

To correct for this, fill out this section of the configuration with the amounts of 16S and 18S molecules in the starting pool. These values can be obtained by running the pool on a TapeStation or Bioanalyzer, which most sequencing cores perform as part of quality control. The pipeline treats these values as the starting ratio, then calculates the ratio observed on the sequencing platform from the final 16S and 18S read counts. It uses the difference to increase the estimated 18S abundance and decrease the estimated 16S abundance accordingly.

Step 5 • Export

Export

Once everything is filled in, export the package you can upload into the Snakemake directory.

Suggested folder structure 
config/
  README.md
  config.yml
  samples.tsv
  prok_and_euk_SSU_amplicon_molarities.tsv
  internal_stds.tsv  # included when internal standards are enabled
  schemas/
    config.schema.yml
    samples.schema.yml
  setup-scripts/
    setup-analysis-dir.sh
Final validation before download

    Your package will contain one complete config/ folder, including the README, schemas, and setup script.

    Enter the path to the downloaded configuration ZIP on your computer.
    For a remote server, use username@hostname:/full/path/. You can also enter a local directory.
    The command updates from the input and output paths above. Run it in Terminal on the computer that contains the downloaded ZIP.
    Copy the configuration package 
    scp ~/Downloads/AMT29-analysis-config.zip <OUTPUT_DESTINATION>

    If you prefer a graphical file-transfer program, upload the ZIP to the same project directory using Cyberduck or FileZilla. Connect with SFTP using the hostname and login details supplied by your institution or server administrator.

    After the ZIP reaches its destination, run these commands from inside the cloned project folder. The first creates the config/ directory and the second refreshes every extracted config-file timestamp to prevent clock-skew errors on the compute system.
    Extract the configuration package 
    unzip AMT29-analysis-config.zip
find config -type f -exec touch {} +

    Run the pipeline

    Once you have downloaded the configuration package, transferred it to the cloned project, and extracted it so that the config/ directory is in the main project folder, activate the Snakemake environment and run the included run_snakemake.sh script.

    Run these commands from the main directory of your cloned project.
    Run run_snakemake.sh 
    conda activate snakemake
bash run_snakemake.sh

    The included script is intentionally minimal. Depending on your HPC system, you may need to add the appropriate executor and resource allocation settings, such as jobs, CPUs, memory, runtime, account, or partition. The script uses the shared conda_envs_dir from config/config.yml; existing matching environments are reused automatically, and only missing or changed environments are created.

    After Snakemake finishes, download the <project-name>-Results-Export/ folder. Its project-specific name prevents downloads from different runs being confused. It contains the formatted data tables, phylum and order summaries, and the self-contained HTML report.
    Final pipeline report 
    AMT29-analysis-Results-Export/AMT29.pipeline-report.html

    Explore your ASV in GRUMP

    If you used this pipeline with the same 515FY–926R primers and the same forward and reverse trim lengths used for GRUMP, you can copy an ASV hash from your results and search for it in the GRUMP Explorer. The explorer will show the ASV’s distribution across GRUMP samples.

    Search an ASV hash in GRUMP