Skip to content

MIC Module 2 tutorial

Jump to bottom
Diana edited this page Aug 25, 2021 · 34 revisions

This lab is based on the 16S rRNA gene analysis method originally developed by Michael Hall and Robert Beiko, and was updated and modified by Diana Haider for the 2021 Microbiome workshop.

Introduction

The lab in this module is a walkthrough of an end-to-end pipeline using the command line interface for the analysis of high-throughput marker gene data including diversity, taxonomic and statistical explorations. The pipeline described is embedded in the latest version of QIIME2 (version 2021.4, and stands for Quantitative Insights into Microbial Ecology) which is a popular microbiome bioinformatics platform for microbial ecology. It uses software packages called plugins that can be developed by anyone and generates zipped artifacts that can be tracked through a provenance file automatically generated by QIIME2. Two scripts, accessible here, make the whole thing work: “download_Willis_et_al_data.sh” retrieves the sequences from the European Bioinformatics Institute through their relatively awesome API, and “qiime_commands.sh” to execute a chain of QIIME2 commands and generate a number of .qza and .qzv files. These files can be unzipped to reveal the sweet contents inside in case you want to (hypothetically, of course) feed them to a non-QIIME2 application. For this tutorial, we will go through each step individually after the import section, starting with the CBW_Willis_reads.qza file.

Commonly used marker genes for microbiome analysis include the 16S ribosomal RNA (rRNA) for prokaryotes, 18S rRNA for eukaryotes, and the internal transcribed spacer (ITS) for fungi. In this tutorial, we will explore a 16S rRNA dataset from environmental samples collected along three stations from the Halifax Line at four different depths, amplified using two different primers. One interesting question is whether the primer choice under- or overrepresents certain taxa groups, and then impacts diversity measurements. Details on the data collection and processing can be found in the original paper that published the data.

After this lab you will be able to: * * *

Getting started

QIIME2 is recommended to be run in an independent conda environment to ensure consistency between different versions of the packages it uses. Here are the instructions for installation. For this instance, a qiime2 environment was already installed, it just needs to be activated.

source activate qiime2-2021.4

You can visualize all available plugins by typing qiime. Your command line should look like this.

To start your analysis, two files are required:

  • Sequencing data
  • Metadata

The sequencing data will be in multiple FASTQ files. Each FASTQ file contains the raw DNA sequences, and information on the sequences, such as an identifier, the DNA sequence and the quality score of each called bases in text format. These files are usually very large, and are an intermediate for further analyses. The metadata file contains information about the samples allowing us to interpret the data in their biological context.

The structure of your folders should look like this

<ROOT>
|-- sequence data/           # Folder with collected data
    `-- raw reads/           # Folder with the raw reads in FASTQ format
        `-- manifest.txt     # Manifest file
        `-- run*_            # Fastqs
    `-- METADATA.txt         # The metadata file
|-- qiime_artifacts/         # Folder containing the results
    `-- CBW_reads.qza        # Reads artifact imported in QIIME2
    `-- DADA2_results/       # Results from denoising the sequencing data
        `-- table.qza
        `-- stats.qza
        `-- representative_sequences.qza
|-- qiime_solutions/         # Intermediate files of expected outputs

In the interest of time, we will skip the import step, and the denoising step. The files required will be linked in the tutorial, and available in the qiime_artifacts folder.

16S data processing

Import raw reads to QIIME2

You can do this step yourself, but for this live lab, we will start with the already imported CBW_Willis_reads.qza file. If you want to fetch the dataset from EBI yourself, you can run bash download_Willis_et_al_data.sh bash script, but it is time consuming (~30 mins). The script downloads the FASTQ files from the EBI database, and creates a manifest file which is a text file with three columns. It guides the import in QIIME2 by providing the sample name, the file path and the direction of the read if the data has paired reads (either reverse or forward).

Once the manifest file is ready, you can just use the qiime import function:

qiime tools import 
	--type 'SampleData[PairedEndSequencesWithQuality]' 
	--input-path import_to_qiime 
	--output-path CBW_Willis_reads

Output

Quality control

Let's start here. The sequences undergo quality control, where we decide on a quality trimming threshold based on the quality profiles of the samples. The goal is to retain the most bases so the read is long enough to be informative, while removing bases that have low quality scores since they can increase the risk of false-positives. QIIME2 has a visualizer that helps you decide on a threshold.

qiime demux summarize
	--p-n 10000 
	--i-data CBW_reads.qza
	--o-visualization qual_viz.qzv

Output

This function randomly samples 10 000 of your reads without replacement, and presents the Q scores in a boxplot. Each bar represents the distribution of quality scores from the subsampled reads for one base position. A drop in quality is expected towards the 3' end of reads due to phasing (Schirmer et al. 2016). Since we have paired-end reads, we need to be cautious about trimming too much since we need to have enough overlapping regions between forward and reverse reads for the joining step. Reads that are shorther than the selected trim length will be discarded, and reads that fail to join are also discarded.

Clustering 16S sequences

From the quality score boxplot, we decided to trim our sequences at 240 for both forward and reverse reads. While the original authors clustered their sequences into OTUs, we will make ASVs using DADA2. DADA2 clusters sequences based on quality score, and discards reads that have a high probability to be an error.

qiime dada2 denoise-paired 
  --i-demultiplexed-seqs CBW_reads.qza
  --p-trunc-len-f 240
  --p-trunc-len-r 240 
  --p-n-threads 4
  --o-representative-sequences representative_sequences.qza
  --o-table unfiltered_table.qza
  --o-denoising-stats denoise_stats.qza
  --verbose #verbose can be added to any command line to print the details of what the computer is doing

Output

The outputs of this step are the foundation for the rest of this analysis.

From the feature table

The feature table is a BIOM table of samples x ASVs. We can inspect the summary statistics:

qiime feature-table summarize --i-table unfiltered_table.qza --o-visualization table_summary

Output

Differential abundance

qiime

From the representative sequences

Taxonomy classification

Naive Bayes Classifiers provide highly accurate taxonomic classification along with low run times for microbiome studies. It is recommended to have your classifier trained on your own data, but there are pre-trained classifiers available on QIIME2's resources page, such as the greengenes and SILVA classifiers. Greengenes is a smaller database for Archaea and Bacteria, and was last updated in 2013. SILVA is a is regurlarly updated and includes eukaryotes. If you have limited memory, it's better to use greengenes, but for this tutorial we will use SILVA.

wget https://data.qiime2.org/2021.4/common/silva-138-99-nb-classifier.qza

qiime feature-classifier classify-sklearn
  --i-classifier silva-138-99-nb-classifier.qza
  --i-reads representative_sequences.qza
  --o-classification taxonomy.qza

qiime taxa barplot
  --i-table table.qza
  --i-taxonomy taxonomy.qza
  --m-metadata-file METADATA.txt
  --o-visualization taxa-bar-plots

Diversity analyses (Barplots and heatmaps)

Generate alpha/beta diversity measures at 2000 sequences/sample

qiime diversity core-metrics-phylogenetic
  --i-phylogeny rooted_tree.qza
  --i-table unfiltered_table.qza
  --p-sampling-depth 2000
  --output-dir diversity_2000
  --m-metadata-file METADATA.txt

Visualization outputs

Artifact outputs

  • bray_curtis_distance_matrix.qza : download
  • bray_curtis_pcoa_results.qza : download
  • evenness_vector.qza : download
  • faith_pd_vector.qza :
  • jaccard_distance_matrix.qza
  • jaccard_pcoa_results.qza
  • observed_features_vector.qza
  • rarefied_table.qza
  • shannon_vector.qza
  • unweighted_unifrac_distance_matrix.qza
  • unweighted_unifrac_pcoa_results.qza
  • weighted_unifrac_distance_matrix.qza
  • weighted_unifrac_pcoa_results.qza
qiime diversity alpha-group-significance
  --i-alpha-diversity diversity_2000/faith_pd_vector.qza
  --m-metadata-file METADATA.txt
  --o-visualization diversity_2000/alpha_PD_significance

qiime diversity alpha-group-significance
  --i-alpha-diversity diversity_2000/shannon_vector.qza
  --m-metadata-file METADATA.txt
  --o-visualization diversity_2000/alpha_shannon_significance

qiime diversity alpha-rarefaction
  --i-table unfiltered_table.qza
  --i-phylogeny rooted_tree.qza
  --p-max-depth 2000
  --m-metadata-file METADATA.txt
  --o-visualization diversity_2000/alpha_rarefaction.qzv

Alignment and phylogeny

qiime alignment mafft
  --i-sequences representative_sequences.qza
  --o-alignment aligned_representative_sequences

qiime alignment mask
  --i-alignment aligned_representative_sequences.qza
  --o-masked-alignment masked_aligned_representative_sequences

qiime phylogeny fasttree
  --i-alignment masked_aligned_representative_sequences.qza
  --o-tree unrooted_tree

qiime phylogeny midpoint-root
  --i-tree unrooted_tree.qza
  --o-rooted-tree rooted_tree

Clone this wiki locally