Skip to content

Amplicon SOP v2 (qiime2 2019.7)

Jump to bottom
Gavin Douglas edited this page Aug 22, 2019 · 63 revisions

THIS SOP IS CURRENTLY UNDER CONSTRUCTION

This standard operating procedure (SOP) is based on the QIIME 2 tutorials and is meant for users who want to quickly run their amplicon data through the Microbiome Helper virtual box image and for internal use. There are more detailed descriptions and tutorials on the QIIME 2 website, which we recommend you check out!

If you use this workflow make sure to keep track of the commands you use locally as this page will be updated over time (see "revisions" above for earlier versions).

Requirements

This workflow assumes that you have QIIME 2 (version 2019.7) installed in a conda environment and that you also have FASTQC and MultiQC installed in this environment.

This workflow also assumes that the input is raw paired-end MiSeq data in demultiplexed FASTQ format located within a folder called raw_data. The filenames are expected to look like this: 105CHE6WT_S325_L001_R1_001.fastq.gz, where each field separated by an _ character corresponds to:

  1. The sample name (105CHE6WT)
  2. The sample number on the run (S325)
  3. The lane number (L001)
  4. The read number (R1 for forward and R2 for reverse)
  5. The set number (001)

Often you will need to re-name your files to match this syntax. However, QIIME 2 accepts many different formats so if your files are not already in this format (e.g. not demultiplexed) you would just need to use slightly different commands for importing your data.

1. First steps

1.1 Format metadata file

You can format the metadata file (which is in a compatible format to the QIIME 1 mapping files) in a spreadsheet and validate the format using a google spreadsheet plugin as described on the QIIME 2 website. The only required column is the sample id column, which should be first. All other columns should correspond to sample metadata. Below, we assume that your metadata filename has been assigned to a bash variable called "$METADATA".

This can be done like so (for example if your file is called metadata.txt and is found in the folder /home/user):

METADATA="/home/user/metadata.txt"

1.2 Set number of cores

Several commands throughout this workflow can run on multiple cores in parallel. How many cores to use in these cases will be saved to the NCORES variable defined below. We set this variable to 1 below, but you can change this to be however many cores you would like to use.

NCORES=1

1.3 Inspect read quality

Visualize sequence quality across raw reads. This is important as a sanity check that your reads are of reasonable quality and to determine how your reads should be trimmed in downstream steps. QIIME 2 comes with a plugin for visualizing read quality, which we will use at a later step. However, when dealing with raw reads the easiest method to use is a combination of FASTQC and MultiQC. Note that these tools are not packaged with QIIME 2 so you will need to install them separately.

This is an important step for identifying outlier samples with especially low quality, read sizes, read depth, and other metrics.

You can run FASTQC with this command (after creating the output directory).

mkdir fastqc_out
fastqc -t $NCORES raw_data/*.fastq.gz -o fastqc_out

If you receive the error Value "FASTQ" invalid for option threads (number expected) (where "FASTQ" is an input filename) then make sure you have defined the NCORES variable correctly and re-run the command.

FASTQC generates a report for each individual file. To aggregate the summary files into a single report we can run MultiQC with these commands (including entering and leaving the FASTQC output directory):

cd fastqc_out
multiqc .
cd ..

The full report is found within multiqc_report.html in the FASTQC output directory. You can view this report in a web-browser on your local computer. The most important reason to visualize this report is to ensure that your samples are of high-quality (based largely on whether the per-base quality is >30 across most of the reads) and that there are no outlier samples.

1.4 Activate QIIME 2 conda environment

You should run the rest of the workflow in a conda environment, which makes sure the correct version of the Python packages required by QIIME 2 are being used. You can activate this conda environment with this command (you may need to swap in source for conda if you get an error):

conda activate qiime2-2019.7

1.5 Import FASTQs as QIIME 2 artifact

To standardize QIIME 2 analyses and to keep track of provenance (i.e. a list of what commands were previously run to produce a file) a special format is used for all QIIME 2 input and output files called an "artifact" (with the extension QZA). The first step is to import the raw reads as a QZA file.

mkdir reads_qza
    
qiime tools import \
   --type SampleData[PairedEndSequencesWithQuality] \
   --input-path raw_data/ \
   --output-path reads_qza/reads.qza \
   --input-format CasavaOneEightSingleLanePerSampleDirFmt

All of the FASTQs are now in the single artifact file reads_qza/reads.qza. This file format can be a little confusing at first, but it is actually just a zipped folder. You can manipulate and explore these files better with the qiime tools utilities (e.g. peek and view).

1.6 Trim primers with cutadapt

Screen out reads that do not begin with primer sequence and remove primer sequence from reads using the cutadapt QIIME 2 plugin. The below primers correspond to the 16S V6-V8 region (Bacteria-specific primer set).

qiime cutadapt trim-paired \
   --i-demultiplexed-sequences reads_qza/reads.qza \
   --p-cores $NCORES \
   --p-front-f ACGCGHNRAACCTTACC \
   --p-front-r ACGGGCRGTGWGTRCAA \
   --p-discard-untrimmed \
   --p-no-indels \
   --o-trimmed-sequences reads_qza/reads_trimmed.qza

If you received the error --p-cores option requires an argument make sure that you set the NCORES bash variable as described above.

If using 16S V4/V5 primers, use these options in the above command:

   --p-front-f GTGYCAGCMGCCGCGGTAA \
   --p-front-r CCGYCAATTYMTTTRAGTTT \

If using 16S V6/V8 Archaeal primers, use these options in the above command:

   --p-front-f TYAATYGGANTCAACRCC \
   --p-front-r CRGTGWGTRCAAGGRGCA \

If using 18S V4 primers, use these options in the above command:

   --p-front-f CYGCGGTAATTCCAGCTC \
   --p-front-r AYGGTATCTRATCRTCTTYG \

If using Fungal ITS2 primers, use these options in the above command:

   --p-front-f GTGAATCATCGAATCTTTGAA \
   --p-front-r TCCTCCGCTTATTGATATGC \

1.7 Summarize trimmed FASTQs

You can run the demux summarize command after trimming the reads to get a report of the number of reads per sample and quality distribution across the reads. This generates a more basic output compared to FASTQC/MultiQC, but is sufficient for this step.

qiime demux summarize \
   --i-data reads_qza/reads_trimmed.qza \
   --o-visualization reads_qza/reads_trimmed_summary.qzv

Note that we gave the output file above the extension .qzv since this a special type of artifact file - a visualization. You can look at the visualization by uploading the file to the QIIME2 view website and clicking on the Interactive Quality Plot tab at the top of the page

2. Denoising the reads into amplicon sequence variants

At this stage the main 2 pipelines you can use are based on either deblur or DADA2. Below we will describe the commands for running deblur. See here if you are interested in running DADA2.

2.1 Join paired-end reads

Currently deblur doesn't support unjoined paired-end reads. Forward and reverse reads can be joined with VSEARCH:

qiime vsearch join-pairs \
   --i-demultiplexed-seqs reads_qza/reads_trimmed.qza \
   --o-joined-sequences reads_qza/reads_trimmed_joined.qza

2.2 Filter out low-quality reads

This command will filter out low-quality reads based on the default options.

qiime quality-filter q-score-joined \
   --i-demux reads_qza/reads_trimmed_joined.qza \
   --o-filter-stats filt_stats.qza \
   --o-filtered-sequences reads_qza/reads_trimmed_joined_filt.qza

It is a good idea at this point just to verify that there haven't been any substantial losses of reads, before going through the whole ASV process, at either the joining or quality-filtering steps above:

qiime demux summarize \
   --i-data reads_qza/reads_trimmed_joined_filt.qza \
   --o-visualization reads_qza/reads_trimmed_joined_filt_summary.qzv

2.3 Running deblur

Run the deblur workflow to correct reads and get amplicon sequence variants (ASVs). Note that the below command will retain singletons, which would have been filtered out unless we set --p-min-reads 1, and is for 16S sequences only. For other amplicon regions, you can either use the denoise-other option in the command and specify a reference database of sequences to use for positive filtering (as in the below versions for 18S and ITS) or use DADA2 as described below (the former still being preferred due to speed considerations).

Note that you will need to trim all sequences to the same length with the --p-trim-length option if you get an error. In order to determine the correct length to trim down to, take a look at the summary visualization generated above in order to select a length to trim back to that maintains the largest/acceptable quantity of reads. Finally, input that length in the deblur option --p-trim-length below and rerun the command.

qiime deblur denoise-16S \
   --i-demultiplexed-seqs reads_qza/reads_trimmed_joined_filt.qza \
   --p-trim-length -1 \
   --p-sample-stats \
   --p-jobs-to-start $NCORES \
   --p-min-reads 1 \
   --output-dir deblur_output

When running 18S and ITS data with deblur you will need to specify custom reference files, which can be downloaded here.

For the 18S version (note the first command imports the deblur database and only needs to be run once on your system):

qiime tools import \
   --input-path /home/shared/rRNA_db/gb203_pr2_all_10_28_99p_clean.fasta \
   --output-path /home/shared/rRNA_db/gb203_pr2_all_10_28_99p_clean.qza \
   --type 'FeatureData[Sequence]'

qiime deblur denoise-other \
   --i-demultiplexed-seqs reads_qza/reads_trimmed_joined_filt.qza \
   --i-reference-seqs /home/shared/rRNA_db/gb203_pr2_all_10_28_99p_clean.qza \
    --p-trim-length -1 \
   --p-sample-stats \
   --p-jobs-to-start $NCORES \
   --p-min-reads 1 \
   --output-dir deblur_output

For the ITS version (note the first command imports the deblur database and only needs to be run once on your system):

qiime tools import \
   --input-path /home/shared/rRNA_db/UNITE_sh_refs_qiime_ver8_99_s_all_02.02.2019.fasta \
   --output-path /home/shared/rRNA_db/UNITE_sh_refs_qiime_ver8_99_s_all_02.02.2019.qza \
   --type 'FeatureData[Sequence]'

qiime deblur denoise-other \
   --i-demultiplexed-seqs reads_qza/reads_trimmed_joined_filt.qza \
   --i-reference-seqs /home/shared/rRNA_db/UNITE_sh_refs_qiime_ver8_99_s_all_02.02.2019.qza \
   --p-trim-length -1 \
   --p-sample-stats \
   --p-jobs-to-start $NCORES \
   --p-min-reads 1 \
   --output-dir deblur_output

2.4 Summarizing deblur output

Once a denoising pipeline has been run you can summarize the output table with the below command, which will create a visualization artifact for you to view.

qiime feature-table summarize \
   --i-table deblur_output/table.qza \
   --o-visualization deblur_output/deblur_table_summary.qzv

We will use this visualization later to determine the the cut-offs for filtering the table below, but for now you should mainly take a look at the visualization to ensure that sufficient reads have been retained after running deblur. This denoising tool filters out reads that either do match to known noise or that do not match with low similarity to the expected amplicon region. If your samples have very low depth after running deblur (compared to the input read depth) this could be a red flag that either you ran the tool incorrectly, you have a lot of noise in your data, or that deblur is inappropriate for your dataset.

3. Assign taxonomy to ASVs

You can assign taxonomy to your ASVs using a Naive-Bayes approach implemented in the scikit learn Python library and the SILVA database. Note that we have trained classifiers for a few different amplicon regions already (which are available in the /home/shared/taxa_classifiers folder of the virtual box image). Below we will describe these classifiers, but if your primer set is not there you should check what classifiers are available on the QIIME 2 website or generate your own.

The classifier filename below is for the V6/V8 B969F-BA1406R region. The regions that we have trained classifiers in advance for are:

  • 16S V4/V5 region (classifier_silva_132_99_16S_V4.V5_515F_926R.qza)
  • 16S V6/V8 region (classifier_silva_132_99_16S_V6.V8_B969F_BA1406R.qza)
  • 16S V6/V8 region targeting archaea (classifier_silva_132_99_16S_V6.V8_A956F_A1401R.qza)
  • 16S V3/V4 region targeting cyanobacteria (classifier_silva_132_99_16S_V3.V4_CYA359F_CYA781R.qza)
  • 18S V4 region (classifier_silva_132_99_18S_V4_E572F_E1009R.qza)
  • Full ITS - fungi only (classifier_sh_refs_qiime_ver8_99_s_02.02.2019_ITS.qza)
  • Full ITS - all eukaryotes (classifier_sh_refs_qiime_ver8_99_s_all_02.02.2019_ITS.qza)

You can download these classifiers here, which were created using these commands.

You can run the taxonomic classification with this command, which is one of the longest running and most memory-intensive command of the SOP. If you receive an error related to insufficient memory you will need to increase your memory usage.

qiime feature-classifier classify-sklearn \
   --i-reads deblur_output/representative_sequences.qza \
   --i-classifier /home/shared/taxa_classifiers/qiime2-2019.7_classifiers/classifier_silva_132_99_16S_V6.V8_B969F_BA1406R.qza \
   --p-n-jobs $NCORES \
   --output-dir taxa

As with all QZA files, you can export the output file to take a look at the classifications and confidence scores:

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

4. Filtering resultant table

Filtering the denoised table is an important step of microbiome data analysis. You can see more details on this process in the QIIME 2 filtering tutorial.

4.1 Filter out rare ASVs

Based on the summary visualization created in step 2.4 above you can choose a cut-off for how frequent a variant needs to be (and optionally how many samples need to have the variant) for it to be retained. One possible choice would be to remove all ASVs that have a frequency of less than 0.1% of the mean sample depth. This cut-off excludes ASVs that are likely due to MiSeq bleed-through between runs (reported by Illumina to be 0.1% of reads). To calculate this cut-off you would identify the mean sample depth in the visualization created in step 2.4, multiply it by 0.001, and round to the nearest integer.

Once you've determined how you would like to filter your table you can do so with this command (X is a placeholder for your choice):

qiime feature-table filter-features \
   --i-table deblur_output/table.qza \
   --p-min-frequency X \
   --p-min-samples 1 \
   --o-filtered-table deblur_output/deblur_table_filt.qza

4.2 Filter out contaminant and unclassified ASVs

Now that we have assigned taxonomy to our ASVs we can use that information to remove ASVs which are likely contaminants or noise based on the taxonomic labels. Two common contaminants in 16S sequencing data are mitochondrial and chloroplast 16S sequences, which can be removed by excluding any ASV which contains those terms in its taxonomic label. It can also be useful to exclude any ASV that is unclassified at the phylum level since these sequences are more likely to be noise (e.g. possible chimeric sequences). Note that if your data has not been classified against SILVA you will need to change D_1__ to be a string that enables phylum-level assignments to be identified or simply omit that line. You may also want to omit that line if you are studying a poorly characterized environment where you have a good chance of identifying novel phyla.

qiime taxa filter-table \
   --i-table deblur_output/deblur_table_filt.qza \
   --i-taxonomy taxa/classification.qza \
   --p-include D_1__ \
   --p-exclude mitochondria,chloroplast \
   --o-filtered-table deblur_output/deblur_table_filt_contam.qza

4.3 Exclude low-depth samples

Often certain samples will have quite low depth after these filtering steps, which can be excluded from downstream analyses since they will largely add noise. There is no single cut-off that works best for all datasets, but researchers often use minimum cut-offs within the range of 1000 to 4000 reads. You can also use a cut-off much lower than this if you want to retain all samples except those that failed entirely (e.g. depth < 50 reads). Ideally you would choose this cut-off after visualizing rarefaction curves to determine at what read depth the richness of your samples plateaus and choose a cut-off as close to this plateau as possible while retaining sufficient sample size for your analyses.

To perform this rarefaction curve analysis you would first need to summarize the filtered table we produced in the last step:

qiime feature-table summarize \
   --i-table deblur_output/deblur_table_filt_contam.qza \
   --o-visualization deblur_output/deblur_table_filt_contam_summary.qzv

From this table you need to determine the maximum depth across your samples. You can then generate the rarefaction curves with this command (where X is a placeholder for the max depth across samples).

qiime diversity alpha-rarefaction \
   --i-table deblur_output/deblur_table_filt_contam.qza \
   --p-max-depth X \
   --p-steps 20 \
   --p-metrics 'observed_otus' \
   --o-visualization rarefaction_curves_test.qzv

Take a look at these curves to help decide on a minimum depth cut-off for retaining samples. Once you decide on a hard cut-off you can exclude samples below this cut-off with this command (where SET_CUTOFF is a placeholder for the minimum depth you select):

qiime feature-table filter-samples \
   --i-table deblur_output/deblur_table_filt_contam.qza \
   --p-min-frequency SET_CUTOFF \
   --o-filtered-table deblur_output/deblur_table_final.qza

Alternatively, if you do not wish to exclude any samples then you can simply make a copy of the QZA file with the final table filename (i.e. cp deblur_output/deblur_table_filt_contam.qza deblur_output/deblur_table_final.qza), since this is the filename used for the remaining commands below.

4.4 Subset and summarize filtered table

Once we have our final filtered table we will need to subset the QZA of ASV sequences to the same set. You can exclude the low frequency ASVs from the sequence file with this command:

qiime feature-table filter-seqs \
   --i-data deblur_output/representative_sequences.qza \
   --i-table deblur_output/deblur_table_final.qza \
   --o-filtered-data deblur_output/rep_seqs_final.qza

Finally, you can make a new summary of the final filtered abundance table:

qiime feature-table summarize \
   --i-table deblur_output/deblur_table_final.qza \
   --o-visualization deblur_output/deblur_table_final_summary.qzv

5. Build tree with SEPP QIIME 2 plugin

SEPP is one method for placing short sequences into a reference phylogenetic tree. This is a useful way of determining a phylogenetic tree for your ASVs. For 16S data you can do this with q2-fragment-insertion using the below command:

qiime fragment-insertion sepp \
   --i-representative-sequences deblur_output/rep_seqs_final.qza \
   --o-tree asvs-tree.qza \
   --o-placements insertion-placements.qza \
   --p-threads $NCORES

You can specify custom reference files to place other amplicons, but the easiest approach for 18S and ITS data is to instead create a de novo tree.

6. Generate rarefaction curves

A key quality control step is to plot rarefaction curves for all of your samples to determine if you performed sufficient sequencing. The below command will generate these plots (X is a placeholder for the maximum depth in your dataset, which you can determine by running the summarize command above). Depending on if you decided to exclude any samples you may not want to re-create rarefaction curves. Note however that the below command will output rarefaction curves for a range of alpha-diversity metrics (including phylogenetic metrics), whereas above the curves were based on richness (referred to as observed_otus) only. Also note that the $METADATA bash variable was defined at the first step of this SOP and simply points to the metadata table.

qiime diversity alpha-rarefaction \
   --i-table deblur_output/deblur_table_final.qza \
   --p-max-depth X \
   --p-steps 20 \
   --i-phylogeny asvs-tree.qza \
   --m-metadata-file $METADATA \
   --o-visualization rarefaction_curves.qzv

For some reason, the QIIME 2 default in the above curves with the metadata file (which you can see in the visualization) is to not give you the option of seeing each sample's rarefaction curve individually (even though this is the default later on in stacked barplots!), only the "grouped" curves by each metadata type. As it can be quite important in data QC to see if you have inconsistent samples, we need to rerun the above command, but this time omitting the metadata file (use the same X for the maximum depth, as above).

qiime diversity alpha-rarefaction \
   --i-table deblur_output/deblur_table_final.qza \
   --p-max-depth X \
   --p-steps 20 \
   --i-phylogeny asvs-tree.qza \
   --o-visualization rarefaction_curves_eachsample.qzv

7. Generate stacked barchart of taxa relative abundances

A more useful output is the interactive stacked bar-charts of the taxonomic abundances across samples, which can be output with this command:

qiime taxa barplot \
   --i-table deblur_output/deblur_table_final.qza \
   --i-taxonomy taxa/classification.qza \
   --m-metadata-file $METADATA \
   --o-visualization taxa/taxa_barplot.qzv

Now the QIIME 2 default in the above barplots is to plot each sample individually. However, in the visualizer you can label the samples according to their metadata types, but the plots are not summed together according to their metadata (as we could do in QIIME1 using the summarize_taxa_through_plots.py script), so we have to rerun the above command after running a new command first to group the samples by metadata - this creates a new feature table (you have to do this for each metadata category) which becomes the new input for the above same taxa barplot command (again, run each time for each metadata category). Note that CATEGORY here below is a placeholder for the text label of your category of interest from the metadata file.

qiime feature-table group \
   --i-table deblur_output/deblur_table_final.qza \
   --p-axis sample \
   --p-mode sum \
   --m-metadata-file $METADATA \
   --m-metadata-column CATEGORY \
   --o-grouped-table deblur_output/deblur_table_final_CATEGORY.qza

qiime taxa barplot \
   --i-table deblur_output/deblur_table_final_CATEGORY.qza \
   --i-taxonomy taxa/classification.qza \
   --m-metadata-file $METADATA \
   --o-visualization taxa/taxa_barplot_CATEGORY.qzv

8. Calculating diversity metrics and generating ordination plots

Common alpha and beta-diversity metrics can be calculated with a single command in QIIME 2. In addition, ordination plots (such as PCoA plots for weighted UniFrac distances) will be generated automatically as well. This command will also rarefy all samples to the sample sequencing depth before calculating these metrics (X is a placeholder for the lowest reasonable sample depth; samples with depth below this cut-off will be excluded).

qiime diversity core-metrics-phylogenetic \
   --i-table deblur_output/deblur_table_final.qza \
   --i-phylogeny asvs-tree.qza \
   --p-sampling-depth X \
   --m-metadata-file $METADATA \
   --p-n-jobs $NCORES \
   --output-dir diversity

You can then produce boxplots comparing the different categories in your metadata file. For example, to create boxplots comparing the Shannon alpha-diversity metric you can use this command:

qiime diversity alpha-group-significance \
   --i-alpha-diversity diversity/shannon_vector.qza \
   --m-metadata-file $METADATA \
   --o-visualization diversity/shannon_compare_groups.qzv

Note that you can also export (see below) this or any other diversity metric file (ending in .qza) and analyze them with a different program.

9. Identifying differentially abundant features with ANCOM

ANCOM is one method to test for differences in the relative abundance of features between sample groupings. It is a compositional approach that makes no assumptions about feature distributions. However, it requires that all features have non-zero abundances so a pseudocount first needs to be added (1 is a typical pseudocount choice):

qiime composition add-pseudocount \
   --i-table deblur_output/deblur_table_final.qza \
   --p-pseudocount 1 \
   --o-composition-table deblur_output/deblur_table_final_pseudocount.qza

Then ANCOM can be run with this command; note that CATEGORY is a placeholder for the text label of your category of interest from the metadata file:

qiime composition ancom \
   --i-table deblur_output/deblur_table_final_pseudocount.qza \
   --m-metadata-file $METADATA \
   --m-metadata-column CATEGORY \
   --output-dir ancom_output

10. Exporting the final abundance, profile and sequence files

Lastly, to get the BIOM file (with associated taxonomy) and FASTA file (one per ASV) for your dataset to plug into other programs you can use the commands below.

To export the FASTA:

qiime tools export \
   --input-path deblur_output/rep_seqs_final.qza \
   --output-path deblur_output_exported

To export the BIOM table (with taxonomy added as metadata):

sed -i -e '1 s/Feature/#Feature/' -e '1 s/Taxon/taxonomy/' taxa/taxonomy.tsv

qiime tools export \
   --input-path deblur_output/deblur_table_final.qza \
   --output-path deblur_output_exported

biom add-metadata \
   -i deblur_output_exported/feature-table.biom \
   -o deblur_output_exported/feature-table_w_tax.biom \
   --observation-metadata-fp taxa/taxonomy.tsv \
   --sc-separated taxonomy

biom convert \
   -i deblur_output_exported/feature-table_w_tax.biom \
   -o deblur_output_exported/feature-table_w_tax.txt \
   --to-tsv \
   --header-key taxonomy

Other resources

There are many other possible QIIME 2 analyses that we recommend you look into. You may also find these other resources useful:

Contact

Useful Links

Main SOPs

Amplicon SOP v2 (qiime2-amplicon-2024.5)

PacBio Amplicon SOP v2 (qiime2-2022.2)

Metagenomics SOP v3

Wet-Lab SOPs on Protocols.io

Old SOPs

Tutorials

Microbiome for beginners

Metagenomics Resources

mSystems paper data

Clone this wiki locally