Intro to Metabarcoding using QIIME2

Description

What is the influence of genotype (intrinsic) and environment (extrinsic) on anemone-associated bacterial communities?

Data: Illumina MiSeq v3 paired-end (2 × 300 bp) reads (FASTQ).

Tools: QIIME 2

Pipeline:

• Section 1: Importing, cleaning and quality control of the data
• Section 2: Taxonomic Analysis
• Section 3: Building a phylogenetic tree
• Section 4: Basic visualisations and statistics
• Section 5: Exporting data for further analysis in R
• Section 6: Extra Information

Tutorial Layout

• There is a Table of contents on the top of the page which can be used to easily navigate through the tutorial by clicking the relevant section.

These black coloured boxes are code blocks. The rectangular boxes in the top
right hand corner of this code block/black box can be used to copy the code to
the clipboard.

Requirements and Preparation

Mode of Delivery

This workshop will be run on the SAIAB lab417 server.

Information on how to access the lab417 server is explained in Introduction to the Unix shell and High Performance Computing.

You will need to use a Google Chrome or Mozilla Firefox web browser to view files in QIIME2 View.

Required Data

• No additional data needs to be downloaded for this workshop - it is all located on the lab417 server. FASTQs are located in the directory raw_data and a metadata (metadata.tsv) file has also been provided.
• If you wish to analyse the data independently at a later stage, it can be downloaded from here. This zipped folder contains both the FASTQs and associated metadata file.
• If you are running this tutorial independently, you can also access the classifier that has been trained specifically for this data from here.

Background

What is the influence of genotype (intrinsic) and environment (extrinsic) on anemone-associated bacterial communities?

The Players

• Exaiptasia diaphana - a shallow-water, marine anemone that is often used in research as a model organism for corals. In this experiment, two genotypes (AIMS1 and AIMS4) of E. diaphana were grown in each of two different environments:
  1. sterile seawater OR
  2. unfiltered control seawater
• The anemone-associated bacterial communities or microbiome - these bacteria live on, or within E. diaphana, and likely consist of a combination of commensals, transients, and long-term stable members, and combined with their host, form a mutually beneficial, stable symbiosis.

The Study

The anemone microbiome contributes to the overall health of this complex system and can evolve in tandem with the anemone host. In this data set we are looking at the impact of intrinsic and extrinsic factors on anemone microbiome composition. After three weeks in either sterile or control seawater (environment), anemones were homogenized and DNA was extracted. There are 23 samples in this data set - 5 from each anemone treatment combination (2 genotypes x 2 environments) and 3 DNA extraction blanks as controls. This data is a subset from a larger experiment.

Dungan AM, van Oppen MJH, and Blackall LL (2021) Short-Term Exposure to Sterile Seawater Reduces Bacterial Community Diversity in the Sea Anemone, Exaiptasia diaphana. Front. Mar. Sci. 7:599314. doi:10.3389/fmars.2020.599314 [Full Text].

QIIME 2 Analysis Platform

Quantitative Insights Into Microbial Ecology 2 (QIIME 2™) is a next-generation microbiome bioinformatics platform that is extensible, free, open source, and community developed. It allows researchers to:

• Automatically track analyses with decentralised data provenance
• Interactively explore data with beautiful visualisations
• Easily share results without QIIME 2 installed
• Plugin-based system — researchers can add in tools as they wish

Qiime2 is built around a few key principles:

Core Architecture

1. Artifact System (.qza files)

• Everything in QIIME2 is stored as "artifacts" - special files that contain:
• Your actual data (sequences, tables, trees, etc.)
• Metadata about how the data was created
• Complete provenance tracking (every step that created this data)

2. Plugin Architecture

• QIIME2 is modular - different tools are "plugins"
• Each plugin has specific "actions" (commands you can run)
• Examples: q2-dada2, q2-feature-classifier, q2-diversity

3. Semantic Types

• Data has specific "types" that QIIME2 understands
• Examples: SampleData[PairedEndSequencesWithQuality], FeatureTable[Frequency], Phylogeny[Rooted]
• This prevents you from using the wrong data type with incompatible methods

Basic QIIME2 Workflow

                    
                    Raw sequences (.fastq) 
                            ↓ [Import]
                    Sequence artifact (.qza)
                            ↓ [Quality control + Denoising]
                    Feature table + Representative sequences (.qza)
                            ↓ [Taxonomic classification]
                    Feature table + Taxonomy table (.qza)
                            ↓ [Diversity analysis]
                    Alpha/Beta diversity results (.qza/.qzv) 
                    
                

Taxonomic Analysis in Detail

Taxonomic analysis answers the question: "What organisms are in my samples?

The Challenge

You have DNA sequences that look like this:

AGGATTAGATACCCTGGTAGTCCACGCCGTAAACGATGTCGATTTGGA...

You need to know: "Is this from E. coli? Bacteroides? Something unknown?"

How QIIME2 Solves This

Step 1: Prepare a Classifier

A classifier is like a "reference library" that has been trained to recognize sequences:

                    
                qiime feature-classifier classify-sklearn \
                     --i-classifier my-classifier.qza \
                     --i-reads representative-seqs.qza \
                     --o-classification taxonomy.qza
                    
                

• Takes known sequences from databases (SILVA, Greengenes, etc.)
• Uses machine learning (Naive Bayes) to learn patterns
• Creates a "trained model" that can classify new sequences

Viewing QIIME2 Visualisations

As this workshop is being run on a remote lab417 server, you will need to download the visual files (*.qzv) to your local computer and view them in QIIME2 View (q2view).

Section 1: Importing, cleaning and quality control of the data

Initial Set up on lab417 server

Using screen command

To ensure that commands continue to run should you get disconnected from your lab417 server, we'll run a screen session.

Starting a screen session on lab417

On lab417, to start a screen session called workshop, type

screen -S workshop

If you get disconnected from your lab417 Instance, follow the instructions here to resume your session.

Should your SSH session on lab417 terminate, once you log back in to your lab417 instance, list running sessions/screens:

screen -ls

If it says (Detached) next to the workshop session in the list, reattach to workshop by:

screen -r workshop

If it says (Attached) next to the workshop session in the list, you can access workshop which is already attached by:

screen -r -d workshop/

Detaching or Terminating a screen session.

To detach from workshop, type ctrl-a ctrl-d while inside the workshop session. Follow the information on the screen and select 1 for Screen mode.

Activate Qiime2 conda environment

Activate the Qiime2 environment with this command:

 source /usr/local/bin/qiime2-activate
                    

Verify that it is working with this command:

qiime info

Symbolic links to workshop data

Data for this workshop is stored in a central location (/opt/courses/qiime/data_files) on the lab417 file system that we will be using. We will use symbolic links (ln -s) to point to it.

cd
ln -s /opt/courses/qiime/data_files/raw_data raw_data
ln -s /opt/courses/qiime/data_files/metadata.tsv metadata.tsv
ln -s /opt/courses/qiime/data_files/silva_138_16s_v5v6_classifier_2021-4.qza silva_138_16s_v5v6_classifier_2021-4.qza

Import data

These samples were sequenced on a single Illumina MiSeq run using v3 (2 × 300 bp) reagents at the Walter and Eliza Hall Institute (WEHI), Melbourne, Australia.

Data from WEHI came as paired-end, demultiplexed, unzipped *.fastq files with adapters still attached. Following the QIIME2 importing tutorial, this is the Casava One Eight format. The files have been renamed to satisfy the Casava format as SampleID_FWDXX-REVXX_L001_R[1 or 2]_001.fastq e.g. CTRLA_Fwd04-Rev25_L001_1_001.fastq.gz. The files were then zipped (.gzip).

Here, the data files (two per sample i.e. forward and reverse reads 1 and 2 respectively) will be imported and exported as a single QIIME 2 artefact file. These samples are already demultiplexed (i.e. sequences from each sample have been written to separate files), so a metadata file is not initially required.

Start by making a new directory analysis to store all the output files from this tutorial. In addition, we will create a subdirectory called seqs to store the exported sequences.

cd
mkdir -p analysis/seqs

Run the command to import the raw data located in the directory raw_data and export it to a single QIIME 2 artefact file, combined.qza.

qiime tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-path raw_data \
--input-format CasavaOneEightSingleLanePerSampleDirFmt \
--output-path analysis/seqs/combined.qza

Remove primers

These sequences still have the primers attached - they need to be removed (using cutadapt) before denoising.

qiime cutadapt trim-paired \
--i-demultiplexed-sequences analysis/seqs/combined.qza \
--p-front-f AGGATTAGATACCCTGGTA \
--p-front-r CRRCACGAGCTGACGAC \
--p-error-rate 0.20 \
--output-dir analysis/seqs_trimmed_data \
--verbose

Section 2: Taxonomic Analysis

Assign taxonomy

Here we will classify each identical read or Amplicon Sequence Variant (ASV) to the highest resolution based on a database. Common databases for bacterial datasets are Greengenes, SILVA, Ribosomal Database Project, or Genome Taxonomy Database. See Porter and Hajibabaei, 2020 for a review of different classifiers for metabarcoding research. The classifier chosen is dependent upon:

1. Previously published data in a field
2. The target region of interest
3. The number of reference sequences for your organism in the database and how recently that database was updated.

A classifier in QIIME2 is a pre-trained machine learning model that assigns taxonomic classifications to DNA sequences (specifically ASVs/OTUs) by comparing them against a reference database. It's one of the most crucial components of microbiome analysis. The classifier takes your unknown DNA sequences and answers the question: "What organism does this sequence come from?"

A classifier has already been trained for you for the V5V6 region of the bacterial 16S rRNA gene using the SILVA database. The next step will take a while to run. The output directory cannot previously exist.

n_jobs = 1 This runs the script using all available cores

qiime feature-classifier classify-sklearn \
--i-classifier silva_138_16s_v5v6_classifier_2021-4.qza \
--i-reads analysis/dada2out/representative_sequences.qza \
--p-n-jobs 1 \
--output-dir analysis/taxonomy \
--verbose
        

Generate a viewable summary file of the taxonomic assignments.

The qiime metadata tabulate command creates interactive, searchable tables from QIIME2 metadata and statistical results based on the input files. It's one of the most versatile visualization commands, essential for examining sample metadata, denoising statistics, diversity metrics, and other tabular data. This command is your data exploration hub - essential for quality control, results interpretation, and preparing data for downstream statistical analysis!

Here for example we generate a viewable summary file of the taxonomic assignments

qiime metadata tabulate \
--m-input-file analysis/taxonomy/classification.qza \
--o-visualization analysis/visualisations/taxonomy.qzv \
--verbose 
        

Copy analysis/visualisations/taxonomy.qzv to your local computer and view in QIIME 2 View (q2view).

Filtering

Filter out reads classified as mitochondria and chloroplast. Unassigned ASVs are retained.

According to QIIME developer Nicholas Bokulich, low abundance filtering (i.e. removing ASVs containing very few sequences) is not necessary under the ASV model.

qiime taxa filter-table \
--i-table analysis/dada2out/table.qza \
--i-taxonomy analysis/taxonomy/classification.qza  \
--p-exclude Mitochondria,Chloroplast \
--o-filtered-table analysis/taxonomy/16s_table_filtered.qza \
--verbose
        

Generate a viewable summary file of the new table to see the effect of filtering.

qiime feature-table summarize \
--i-table analysis/taxonomy/16s_table_filtered.qza \
--m-sample-metadata-file metadata.tsv \
--o-visualization analysis/visualisations/16s_table_filtered.qzv \
--verbose
        

Copy analysis/visualisations/16s_table_filtered.qzv to your local computer and view in QIIME 2 View (q2view).

Section 3: Build a phylogenetic tree

The next step does the following:

1. Perform an alignment on the representative sequences.
2. Mask sites in the alignment that are not phylogenetically informative.
3. Generate a phylogenetic tree.
4. Apply mid-point rooting to the tree.

A phylogenetic tree is necessary for any analyses that incorporates information on the relative relatedness of community members, by incorporating phylogenetic distances between observed organisms in the computation. This would include any beta-diversity analyses and visualisations from a weighted or unweighted Unifrac distance matrix.

mkdir analysis/tree
analysis/phylogeny
        

Use one thread only (which is the default action) so that identical results can be produced if rerun.

qiime phylogeny align-to-tree-mafft-fasttree \
--i-sequences analysis/dada2out/representative_sequences.qza \
--o-alignment analysis/tree/aligned_16s_representative_seqs.qza \
--o-masked-alignment analysis/tree/masked_aligned_16s_representative_seqs.qza \
--o-tree analysis/tree/16s_unrooted_tree.qza \
--o-rooted-tree analysis/tree/16s_rooted_tree.qza \
--p-n-threads 1 \
--verbose
        

Section 4: Basic visualisations and statistics

ASV relative abundance bar charts

Create bar charts to compare the relative abundance of ASVs across samples.

qiime taxa barplot \
--i-table analysis/taxonomy/16s_table_filtered.qza \
--i-taxonomy analysis/taxonomy/classification.qza \
--m-metadata-file metadata.tsv \
--o-visualization analysis/visualisations/barchart.qzv \
--verbose
        

Copy analysis/visualisations/barchart.qzv to your local computer and view in QIIME 2 View (q2view). Try selecting different taxonomic levels and metadata-based sample sorting.

Rarefaction curves

Generate rarefaction curves to determine whether the samples have been sequenced deeply enough to capture all the community members.

qiime diversity alpha-rarefaction \
--i-table analysis/taxonomy/16s_table_filtered.qza \
--i-phylogeny analysis/tree/16s_rooted_tree.qza \
--p-max-depth 9062 \
--m-metadata-file metadata.tsv \
--o-visualization analysis/visualisations/16s_alpha_rarefaction.qzv \
--verbose
        

Copy analysis/visualisations/16s_alpha_rarefaction.qzv to your local computer and view in QIIME 2 View (q2view).

Alpha and beta diversity analysis

The following is taken directly from the Moving Pictures tutorial and adapted for this data set. QIIME 2’s diversity analyses are available through the q2-diversity plugin, which supports computing alpha- and beta- diversity metrics, applying related statistical tests, and generating interactive visualisations. We’ll first apply the core-metrics-phylogenetic method, which rarefies a FeatureTable[Frequency] to a user-specified depth, computes several alpha- and beta- diversity metrics, and generates principle coordinates analysis (PCoA) plots using Emperor for each of the beta diversity metrics.

The metrics computed by default are:

• Alpha diversity (operate on a single sample (i.e. within sample diversity)).
• Shannon’s diversity index (a quantitative measure of community richness)
• Observed OTUs (a qualitative measure of community richness)
• Faith’s Phylogenetic Diversity (a qualitative measure of community richness that incorporates phylogenetic relationships between the features)
• Evenness (or Pielou’s Evenness; a measure of community evenness)
• Beta diversity (operate on a pair of samples (i.e. between sample diversity)).
• Jaccard distance (a qualitative measure of community dissimilarity)
• Bray-Curtis distance (a quantitative measure of community dissimilarity)
• unweighted UniFrac distance (a qualitative measure of community dissimilarity that incorporates phylogenetic relationships between the features)
• weighted UniFrac distance (a quantitative measure of community dissimilarity that incorporates phylogenetic relationships between the features)

An important parameter that needs to be provided to this script is --p-sampling-depth, which is the even sampling (i.e. rarefaction) depth that was determined above. As most diversity metrics are sensitive to different sampling depths across different samples, this script will randomly subsample the counts from each sample to the value provided for this parameter. For example, if --p-sampling-depth 500 is provided, this step will subsample the counts in each sample without replacement, so that each sample in the resulting table has a total count of 500. If the total count for any sample(s) are smaller than this value, those samples will be excluded from the diversity analysis. Choosing this value is tricky. We recommend making your choice by reviewing the information presented in the 16s_table_filtered.qzv file that was created above. Choose a value that is as high as possible (so more sequences per sample are retained), while excluding as few samples as possible.

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny analysis/tree/16s_rooted_tree.qza \
  --i-table analysis/taxonomy/16s_table_filtered.qza \
  --p-sampling-depth 5583 \
  --m-metadata-file metadata.tsv \
  --output-dir analysis/diversity_metrics
        

Copy the .qzv files created from the above command into the visualisations subdirectory.

cp analysis/diversity_metrics/*.qzv analysis/visualisations
        

To view the differences between sample composition using unweighted UniFrac in ordination space, copy analysis/visualisations/unweighted_unifrac_emperor.qzv to your local computer and view in QIIME 2 View (q2view).

Other metrics calculated above can also be viewed by copying the files in analysis/visualisations/ to your local computer and view in QIIME 2 View (q2view).

Next, we’ll test for associations between categorical metadata columns and alpha diversity data. We’ll do that here for the Faith Phylogenetic Diversity (a measure of community richness) and Evenness metrics.

qiime diversity alpha-group-significance \
  --i-alpha-diversity analysis/diversity_metrics/faith_pd_vector.qza \
  --m-metadata-file metadata.tsv \
  --o-visualization analysis/visualisations/faith-pd-group-significance.qzv
        

Copy analysis/visualisations/faith-pd-group-significance.qzv to your local computer and view in QIIME 2 View (q2view).

qiime diversity alpha-group-significance \
  --i-alpha-diversity analysis/diversity_metrics/evenness_vector.qza \
  --m-metadata-file metadata.tsv \
  --o-visualization analysis/visualisations/evenness-group-significance.qzv
        

Copy analysis/visualisations/evenness-group-significance.qzv to your local computer and view in QIIME 2 View (q2view).

Finally, we’ll analyse sample composition in the context of categorical metadata using a permutational multivariate analysis of variance (PERMANOVA, first described in Anderson (2001)) test using the beta-group-significance command. The following commands will test whether distances between samples within a group, such as samples from the same genotype, are more similar to each other then they are to samples from the other groups. If you call this command with the --p-pairwise parameter, as we’ll do here, it will also perform pairwise tests that will allow you to determine which specific pairs of groups (e.g., AIMS1 and AIMS4) differ from one another, if any. This command can be slow to run, especially when passing --p-pairwise, since it is based on permutation tests. So, unlike the previous commands, we’ll run beta-group-significance on specific columns of metadata that we’re interested in exploring, rather than all metadata columns to which it is applicable. Here we’ll apply this to our unweighted UniFrac distances, using two sample metadata columns, as follows.

qiime diversity beta-group-significance \
  --i-distance-matrix analysis/diversity_metrics/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column Genotype \
  --o-visualization analysis/visualisations/unweighted-unifrac-genotype-significance.qzv \
  --p-pairwise
        

Copy analysis/visualisations/unweighted-unifrac-genotype-significance.qzv to your local computer and view in QIIME 2 View (q2view).

qiime diversity beta-group-significance \
  --i-distance-matrix analysis/diversity_metrics/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column Environment \
  --o-visualization analysis/visualisations/unweighted-unifrac-environment-significance.qzv \
  --p-pairwise
        

Copy analysis/visualisations/unweighted-unifrac-environment-significance.qzv to your local computer and view in QIIME 2 View (q2view).

Section 5: Exporting data for further analysis in R

You need to export your ASV table, taxonomy table, and tree file for analyses in R. Many file formats can be accepted.

Export unrooted tree as .nwk format as required for the R package phyloseq.

qiime tools export \
--input-path analysis/tree/16s_unrooted_tree.qza \
--output-path analysis/export

Create a BIOM table with taxonomy annotations. A FeatureTable[Frequency] artefact will be exported as a BIOM v2.1.0 formatted file.

qiime tools export \
--input-path analysis/taxonomy/16s_table_filtered.qza \
--output-path analysis/export
            

Then export BIOM to TSV

biom convert \
-i analysis/export/feature-table.biom \
-o analysis/export/feature-table.tsv \
--to-tsv
                

Export Taxonomy as TSV

qiime tools export \
--input-path analysis/taxonomy/classification.qza \
--output-path analysis/export
                

Delete the header lines of the .tsv files

sed '1d' analysis/export/taxonomy.tsv > analysis/export/taxonomy_noHeader.tsv
sed '1d' analysis/export/feature-table.tsv > analysis/export/feature-table_noHeader.tsv
                

Some packages require your data to be in a consistent order, i.e. the order of your ASVs in the taxonomy table rows to be the same order of ASVs in the columns of your ASV table. It’s recommended to clean up your taxonomy file. You can have blank spots where the level of classification was not completely resolved.

Author Information

This tutorial was modifed from the original developed by Ashley Dungan and Gayle Philip, School of Biosciences, University of Melbourne.