Parkinson’s Mouse Tutorial¶

Hypothesis¶

This tutorial will explore the hypothesis that the genetic background of a humanized mouse influences the microbial community. However, we’ll also need to consider other confounders which might drive the shape of the microbiome instead of the mouse genotype.

Set up¶

This tutorial assumes that you have QIIME 2 installed according to the installation instructions.

Before running the tutorial, you will need to make a directory for the tutorial data and navigate into that directory.

mkdir ./mouse_tutorial
cd ./mouse_tutorial

Metadata¶

Before starting any analysis, it’s important to be familiar with the metadata. In this study, the metadata file contains 7 columns.

Even though the mouse ID looks like a number, we will specify that it is categorical using the #q2_type directive.

The metadata is available as a Google Sheet, or you can download it directly and save it as a TSV (tab-separated values) file.

wget \
  -O "metadata.tsv" \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/sample_metadata.tsv"

curl -sL \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/sample_metadata.tsv" > \
  "metadata.tsv"

The sample metadata will be used throughout the tutorial. Let’s run our first QIIME 2 command, to summarize and explore the metadata.

qiime metadata tabulate \
  --m-input-file metadata.tsv \
  --o-visualization metadata.qzv

Importing data into QIIME 2¶

In QIIME 2, all data is structured as an Artifact of a specific semantic type. Artifacts contain the data as well as information about the data, including a record of the original data and the tools used to process it. This allows for better tracking of how you actually got to where you are in your analysis. You can learn more about common QIIME 2 Artifacts and semantic types here.

Our samples were amplified using the EMP 515f-806r primers and sequenced on an Illumina MiSeq with a 2x150bp kit. The hypervariable region covered by the primers we used is 290bp long, so with 150bp reads our sequences will be slightly too short to be able to do paired-end analysis downstream. Therefore, we’re going to work with single-end sequences. We will work with a version of the samples which have already been demultiplexed, for example, by the sequencing center. If you need to demultiplex your sequences, the Moving Pictures tutorial describes how to demultiplex sequences if they were sequenced using the Earth Microbiome Project protocol.

We will import the sequences as SampleData[SequencesWithQuality], which is the demultiplexed single-end sequence format. If we wanted to import paired-sequences, we would specify the semantic type SampleData[PairedEndSequencesWithQuality]. We will import the sequences using the sample manifest format, a versatile way to import demultiplexed data in QIIME 2. We create a tab-separated sample manifest file that maps the sample name we want to use in QIIME 2 to the path of the sequence file. The benefit is that the demultiplexed sequence files can be named anything you want; there are not fixed assumptions about the conventions, and the file names do not dictate the final name. When QIIME 2 reads the file, it ignores any line prefixed with the # symbol. The first line that doesn’t contain a # is the header line and must be sample-id<TAB>absolute-filepath. The sample order after the header line does not matter. Read more about importing data into QIIME 2 artifacts here and more about sample metadata formatting requirements here.

Let’s start by downloading the manifest and corresponding sequences.

wget \
  -O "manifest.tsv" \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/manifest"

curl -sL \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/manifest" > \
  "manifest.tsv"

wget \
  -O "demultiplexed_seqs.zip" \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/demultiplexed_seqs.zip"

curl -sL \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/demultiplexed_seqs.zip" > \
  "demultiplexed_seqs.zip"

You’ll need to unzip sequence archive you just downloaded:

unzip demultiplexed_seqs.zip

You can use the head command to check the first few lines of the sample manifest.

head manifest.tsv

stdout:

sample-id	absolute-filepath
recip.220.WT.OB1.D7	$PWD/demultiplexed_seqs/10483.recip.220.WT.OB1.D7_30_L001_R1_001.fastq.gz
recip.290.ASO.OB2.D1	$PWD/demultiplexed_seqs/10483.recip.290.ASO.OB2.D1_27_L001_R1_001.fastq.gz
recip.389.WT.HC2.D21	$PWD/demultiplexed_seqs/10483.recip.389.WT.HC2.D21_1_L001_R1_001.fastq.gz
recip.391.ASO.PD2.D14	$PWD/demultiplexed_seqs/10483.recip.391.ASO.PD2.D14_5_L001_R1_001.fastq.gz
recip.391.ASO.PD2.D21	$PWD/demultiplexed_seqs/10483.recip.391.ASO.PD2.D21_1_L001_R1_001.fastq.gz
recip.391.ASO.PD2.D7	$PWD/demultiplexed_seqs/10483.recip.391.ASO.PD2.D7_15_L001_R1_001.fastq.gz
recip.400.ASO.HC2.D14	$PWD/demultiplexed_seqs/10483.recip.400.ASO.HC2.D14_32_L001_R1_001.fastq.gz
recip.401.ASO.HC2.D7	$PWD/demultiplexed_seqs/10483.recip.401.ASO.HC2.D7_22_L001_R1_001.fastq.gz
recip.403.ASO.PD2.D21	$PWD/demultiplexed_seqs/10483.recip.403.ASO.PD2.D21_31_L001_R1_001.fastq.gz

When using this manifest format, a sample name can only appear in one line and can only map to one sequencing file per column (one column for single-end, two columns for paired-end). The absolute-filepath for each sample must be an absolute path, which specifies the “complete” location of the file. We do that here using the $PWD variable, which expands the current directory in absolute terms.

We’ll use the manifest to import our data.

qiime tools import \
  --type "SampleData[SequencesWithQuality]" \
  --input-format SingleEndFastqManifestPhred33V2 \
  --input-path ./manifest.tsv \
  --output-path ./demux_seqs.qza

Let’s check the sequences and the sequencing depth of the samples using the qiime demux summarize command. It provides information about the number of sequences in each sample, as well as the quality of the sequences.

Before running the command, let’s review the help documentation to make sure we understand the arguments for the command.

qiime demux summarize --help

Based on the documentation, we should specify the file (Artifact) with the demultiplexed sequences for the --i-data argument, since this expects data of semantic type SampleData[SequencesWithQuality]. We’ll specify the location we want to save the visualization to by specifying the output path to --o-visualization.

The help documentation is a good reference for any command, and the first place to look if you’re getting errors, especially errors about parameters.

qiime demux summarize \
  --i-data ./demux_seqs.qza \
  --o-visualization ./demux_seqs.qzv

Sequence quality control and feature table¶

There are several ways to construct a feature table in QIIME 2. The first major choice to make is to work with Operational Taxonomic Units (OTUs) or Amplicon Sequence Variants (ASVs). OTUs have been widely used in microbiome research since the mid 2010s, and assign sequences to clusters either based on a reference database or de novo assignment. QIIME 2 offers clustering through q2-vsearch and q2-dbOTU plug-ins, currently.

ASVs are a more recent development and provide better resolution in features than traditional OTU-based methods. ASVs can separate features based on differences of a single nucleotide in sequences of 400 bp or more, a resolution not possibly even with 99% identity OTU clustering. QIIME 2 currently offers denoising via DADA2 (q2-dada2) and Deblur (q2-deblur). The major differences in the algorithms and motivation for denoising are nicely described in Nearing et al, 2018.

It is worth noting in either case that denoising to ASVs and clustering to OTUs are separate, but parallel steps. A choice should be made for a single pathway: either denoising or OTU based clustering; it is not recommended to combine the steps.

In this tutorial, we’ll denoise with DADA2 (using single-end sequences). Please see the Atacama Soil tutorial for an example of using DADA2 on paired-end sequences. For those interested in using Deblur, you can refer to the Moving Pictures tutorial and Alternative methods of read joining tutorial for running Deblur on single- and paired-end sequences, respectively.

The qiime dada2 denoise-single method requires us to set the --p-trunc-len parameter. This controls the length of the sequences and should be selected based on a drop in quality scores. In our dataset, the quality scores are relatively evenly distributed along the sequencing run, so we’ll use the full 150 bp sequences. However, the selection of the trim length is a relatively subjective measurement and relies on the decision making capacity of the analyst.

qiime dada2 denoise-single \
  --i-demultiplexed-seqs ./demux_seqs.qza \
  --p-trunc-len 150 \
  --o-table ./dada2_table.qza \
  --o-representative-sequences ./dada2_rep_set.qza \
  --o-denoising-stats ./dada2_stats.qza

We can also review the denoising statistics using the qiime metadata tabulate command.

qiime metadata tabulate \
  --m-input-file ./dada2_stats.qza  \
  --o-visualization ./dada2_stats.qzv

qiime feature-table summarize \
  --i-table ./dada2_table.qza \
  --m-sample-metadata-file ./metadata.tsv \
  --o-visualization ./dada2_table.qzv

Generating a phylogenetic tree for diversity analysis¶

QIIME 2 analysis allows the use of phylogenetic trees for diversity metrics such as Faith’s Phylogenetic Diversity and UniFrac distance. The tree provides an inherent structure to the data, allowing us to consider an evolutionary relationship between organisms.

QIIME 2 offers several ways to construct a phylogenetic tree. For this tutorial, we’re going to create a fragment insertion tree using the q2-fragment-insertion plugin. The authors of the fragment insertion plugin suggest that it can outperform traditional alignment based methods based on short Illumina reads by alignment against a reference tree built out of larger sequences. Our command, qiime fragment-insertion sepp will use the representative sequences (a FeatureData[Sequence] artifact) we generated during denoising to create a phylogenetic tree where the sequences have been inserted into the greengenes 13_8 99% identity reference tree backbone.

First, we will download the reference database:

wget \
  -O "sepp-refs-gg-13-8.qza" \
  "https://data.qiime2.org/2024.2/common/sepp-refs-gg-13-8.qza"

curl -sL \
  "https://data.qiime2.org/2024.2/common/sepp-refs-gg-13-8.qza" > \
  "sepp-refs-gg-13-8.qza"

qiime fragment-insertion sepp \
  --i-representative-sequences ./dada2_rep_set.qza \
  --i-reference-database sepp-refs-gg-13-8.qza \
  --o-tree ./tree.qza \
  --o-placements ./tree_placements.qza \
  --p-threads 1  # update to a higher number if you can

Alpha Rarefaction and Selecting a Rarefaction Depth¶

We now have a feature table (observation matrix) of ASVs in each sample and a phylogenetic tree representing those ASVs, so we’re almost ready to perform various analyses of microbial diversity. However, first we must normalize our data to account for uneven sequencing depth between samples.

Although sequencing depth in a microbiome sample does not directly relate to the original biomass in a community, the relative sequencing depth has a large impact on observed communities (Weiss et al, 2017). Therefore, for most diversity metrics, a normalization approach is needed.

Rarefaction is a method for normalization via sub-sampling without replacement and is commonly used as a workaround for the issue of uneven sequencing depth. Rarefaction occurs in two steps: first, samples which are below the rarefaction depth are filtered out of the feature table. Then, all remaining samples are subsampled without replacement to get to the specified sequencing depth. It’s both important and sometimes challenging to select a rarefaction depth for diversity analyses. Several strategies exist to figure out an appropriate rarefaction depth - we will primarily consider alpha rarefaction in this tutorial, because it is a data-driven way to approach the problem.

We’ll use qiime diversity alpha-rarefaction to subsample the ASV table at different depths (between --p-min-depth and --p-max-depth) and calculate the alpha diversity using one or more metrics (--p-metrics). When we checked the feature table, we found that the sample with the fewest sequences in the denoised table has 347 features and that the sample with the most has 4996 features. We want to set a maximum depth close to the maximum number of sequences. We also know that if we look at a sequencing depth around 4250 sequences per sample, we’ll be looking at information from 25 samples. So, let’s set this as our maximum sequencing depth.

By default, 10 rarefied tables are calculated at each sampling depth to provide an error estimate. This can be adjusted using the --p-iterations parameter. We can check and see if there is a relationship between the alpha diversity and metadata by specifying the metadata file for the --m-metadata-file parameter.

qiime diversity alpha-rarefaction \
  --i-table ./dada2_table.qza \
  --m-metadata-file ./metadata.tsv \
  --o-visualization ./alpha_rarefaction_curves.qzv \
  --p-min-depth 10 \
  --p-max-depth 4250

The visualization file will display two plots. The upper plot will display the alpha diversity (observed features or shannon) as a function of the sampling depth. This is used to determine whether the richness or evenness has saturated based on the sampling depth. The rarefaction curve should “level out” as you approach the maximum sampling depth. Failure to do so, especially with a diversity-only metric such as observed features or Faith’s PD diversity, may indicate that the richness in the samples has not been fully saturated.

The second plot shows the number of samples in each metadata category group at each sampling depth. This is useful to determine the sampling depth where samples are lost, and whether this may be biased by metadata column group values. Remember that rarefaction is a two-step process and samples that do not meet the rarefaction depth are filtered out of the table. We can use the curves to look at the number of samples by different metadata columns.

If you’re still unsure of the rarefaction depth, you can also use the sample summary to look at which samples are lost by supplying sample metadata to the feature table summary.

After we’ve looked through the data, we need to select a rarefaction depth. In general, selecting a rarefaction depth is a subjective process that requires that analyst’s discretion. Selecting a rarefaction depth is an exercise in minimizing sequence loss while maximizing the sequences retained for diversity analysis. For high-biomass samples (fecal, oral, etc), a general best estimate is a rarefaction depth of no less than 1000 sequences per sample. In low biomass samples where sequencing is shallower, a lower rarefaction depth may be selected although it’s important to keep in mind that the diversity measurements on these samples will be quite noisy and the overall quality will be low.

Diversity analysis¶

The first step in hypothesis testing in microbial ecology is typically to look at within- (alpha) and between-sample (beta) diversity. We can calculate diversity metrics, apply appropriate statistical tests, and visualize the data using the q2-diversity plugin.

We’ll start by using the qiime diversity core-metrics-phylogenetic method, which rarifies the input feature table, calculates several commonly used alpha- and beta-diversity metrics, and produces principal coordinate analysis (PCoA) visualizations in Emperor for the beta diversity metrics. By default, the metrics computed are:

• Alpha Diversity

  • Shannon’s diversity index

  • Observed Features

  • Faith’s phylogenetic diversity

  • Pielou’s evenness

• Beta Diversity

  • Jaccard distance

  • Bray-Curtis distance

  • Unweighted UniFrac distance

  • Weighted UniFrac distance

There is a very good discussion of diversity metrics and their meanings in a forum post by Stephanie Orchanian.

The qiime diversity core-metrics-phylogenetic method wraps several other methods, and it’s worthwhile to note that the steps can also be executed independently.

One important consideration for diversity calculations is the rarefaction depth. Above, we used the alpha rarefaction visualization and the sample summary visualization to pick a rarefaction depth. So, for these analyses, we’ll use a depth of 2000 sequences per sample.

qiime diversity core-metrics-phylogenetic \
  --i-table ./dada2_table.qza \
  --i-phylogeny ./tree.qza \
  --m-metadata-file ./metadata.tsv \
  --p-sampling-depth 2000 \
  --output-dir ./core-metrics-results

qiime diversity alpha-group-significance \
  --i-alpha-diversity ./core-metrics-results/faith_pd_vector.qza \
  --m-metadata-file ./metadata.tsv \
  --o-visualization ./core-metrics-results/faiths_pd_statistics.qzv

qiime diversity alpha-group-significance \
 --i-alpha-diversity ./core-metrics-results/evenness_vector.qza \
 --m-metadata-file ./metadata.tsv \
 --o-visualization ./core-metrics-results/evenness_statistics.qzv

qiime longitudinal anova \
  --m-metadata-file ./core-metrics-results/faith_pd_vector.qza \
  --m-metadata-file ./metadata.tsv \
  --p-formula 'faith_pd ~ genotype * donor_status' \
  --o-visualization ./core-metrics-results/faiths_pd_anova.qzv

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column donor \
  --o-visualization core-metrics-results/unweighted-unifrac-donor-significance.qzv

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column donor \
  --o-visualization core-metrics-results/weighted-unifrac-donor-significance.qzv

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column cage_id \
  --o-visualization core-metrics-results/unweighted-unifrac-cage-significance.qzv \
  --p-pairwise

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column cage_id \
  --o-visualization core-metrics-results/weighted-unifrac-cage-significance.qzv \
  --p-pairwise

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column cage_id \
  --o-visualization core-metrics-results/weighted-unifrac-cage-significance_disp.qzv \
  --p-method permdisp

qiime diversity adonis \
  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --o-visualization core-metrics-results/unweighted_adonis.qzv \
  --p-formula genotype+donor

Taxonomic classification¶

Up until now we have been performing diversity analyses directly on ASVs; in other words, we have assessed the similarity between samples based purely on the unique sequence variants that were observed in each sample. In most experiments we would like to get a sense of what microbial taxa are present — to identify ASVs and give them “names”. To do this, we’ll use the q2-feature-classifier plugin to classify ASVs taxonomically.

For this analysis, we’ll use a pre-trained naive Bayes machine-learning classifier that was trained to differentiate taxa present in the 99% Greengenes 13_8 reference set trimmed to 250 bp of the V4 hypervariable region (corresponding to the 515F-806R primers). This classifier works by identifying k-mers that are diagnostic for particular taxonomic groups, and using that information to predict the taxonomic affiliation of each ASV. We can download the pre-trained classifier here:

wget \
  -O "gg-13-8-99-515-806-nb-classifier.qza" \
  "https://data.qiime2.org/2024.2/common/gg-13-8-99-515-806-nb-classifier.qza"

curl -sL \
  "https://data.qiime2.org/2024.2/common/gg-13-8-99-515-806-nb-classifier.qza" > \
  "gg-13-8-99-515-806-nb-classifier.qza"

It’s worth noting that Naive Bayes classifiers perform best when they’re trained for the specific hypervariable region amplified. You can train a classifier specific for your dataset based on the training classifiers tutorial or download classifiers for other datasets from the QIIME 2 resource page. Classifiers can be re-used for consistent versions of the underlying packages, database, and region of interest.

qiime feature-classifier classify-sklearn \
  --i-reads ./dada2_rep_set.qza \
  --i-classifier ./gg-13-8-99-515-806-nb-classifier.qza \
  --o-classification ./taxonomy.qza

Now, let’s review the taxonomy associated with the sequences using the qiime metadata tabulate method.

qiime metadata tabulate \
  --m-input-file ./taxonomy.qza \
  --o-visualization ./taxonomy.qzv

Let’s also tabulate the representative sequences (FeatureData[Sequence]). Tabulating the representative sequences will allow us to see the sequence assigned to the identifier and interactively blast the sequence against the NCBI database.

qiime feature-table tabulate-seqs \
  --i-data ./dada2_rep_set.qza \
  --o-visualization ./dada2_rep_set.qzv

Taxonomy barchart¶

Since we saw a difference in diversity in this dataset, we may want to look at the taxonomic composition of these samples. To visualize this, we will build a taxonomic barchart of the samples we analyzed in the diversity dataset.

Before doing this, we will first filter out any samples with fewer features than our rarefaction threshold (2000). We can filter samples using the q2-feature-table plugin with the filter-samples method. This lets us filter our table based on a variety of criteria such as the number of counts (frequency, --p-min-frequency and --p-max-frequency), number of features (--p-min-features and --p-max-features), or sample metadata (--p-where). See the filtering tutorial for more details and examples.

For this example, we need to filter out samples with fewer sequences than our rarefaction depth.

qiime feature-table filter-samples \
  --i-table ./dada2_table.qza \
  --p-min-frequency 2000 \
  --o-filtered-table ./table_2k.qza

Now, let’s use the filtered table to build an interactive barplot of the taxonomy in each sample.

qiime taxa barplot \
  --i-table ./table_2k.qza \
  --i-taxonomy ./taxonomy.qza \
  --m-metadata-file ./metadata.tsv \
  --o-visualization ./taxa_barplot.qzv

Differential abundance with ANCOM-BC¶

Many microbiome investigators are interested in testing whether individual ASVs or taxa are more or less abundant in different sample groups. This is known as differential abundance. Microbiome data present several challenges for performing differential abundance using convential methods. Microbiome abundance data are inherently sparse (have a lot of zeros) and compositional (everything adds up to 1). Because of this, traditional statistical methods that you may be familiar with, such as ANOVA or t-tests, are not appropriate for performing differential abundance tests of microbiome data and lead to a high false-positive rate. ANCOM-BC is a compositionally-aware linear regression-based alternative that allows for testing differentially abundant features across groups that also handles bias correction. Microbiome data are typically subject to two sources of biases: unequal sampling fractions (sample-specific biases) and differential sequencing efficiencies (taxon-specific biases). Methodologies included in the ANCOM-BC package are designed to correct these biases and construct statistically consistent estimators. If you’re unfamiliar with the technique, it’s worthwhile to review the ANCOM-BC paper to better understand the method. We also note that accurately identifying features that are differentially abundant across sample types in microbiome data is a challenging problem and an open area of research. New approaches for differential abundance testing are regularly introduced, and it’s worth assessing the current state of the field when performing differential abundance testing to see if there are new methods that might be useful for your data.

Before we begin, we will filter out low abundance/low prevalence ASVs. Filtering can provide better resolution and limit false discovery rate (FDR) penalty on features that are too far below the noise threshhold to be applicable to a statistical test. A feature that shows up with 10 counts could be: a real feature that is present only in that sample; a feature that’s present in several samples but only got amplified and sequenced in one sample because PCR is a somewhat stochastic process; or it may be noise. It’s not possible to tell, so feature-based analysis may be better after filtering low abundance features. However, filtering also shifts the composition of a sample, further disrupting the relationship. Here, the filtering is performed as a trade off between the model, computational efficiency, and statistical practicality.

qiime feature-table filter-features \
  --i-table ./table_2k.qza \
  --p-min-frequency 50 \
  --p-min-samples 4 \
  --o-filtered-table ./table_2k_abund.qza

ANCOM-BC fundamentally operates on a FeatureTable[Frequency], which contains the frequencies of features in each sample. ANCOM-BC accounts for sampling fraction by introducing a sample-specific offset term in a linear regression framework that is estimated from the observed data. The offset term serves as the bias correction, and the linear regression framework in log scale is analogous to log-ratio transformation to deal with the compositionality of microbiome data.

Let’s use ANCOM-BC to check whether there is a difference in the gut microbiome of the mice based on their donor and their genetic background. We’ll set a significance threshold for each group to narrow down the features based on their q-value (i.e., false discovery rate corrected p-value). Features above that threshold will not be included in each group’s visualization.

qiime composition ancombc \
  --i-table ./table_2k_abund.qza \
  --m-metadata-file ./metadata.tsv \
  --p-formula 'donor' \
  --o-differentials ./ancombc_donor.qza

qiime composition da-barplot \
  --i-data ./ancombc_donor.qza \
  --p-significance-threshold 0.001 \
  --o-visualization da_barplot_donor.qzv

qiime composition ancombc \
  --i-table ./table_2k_abund.qza \
  --m-metadata-file ./metadata.tsv \
  --p-formula 'genotype' \
  --o-differentials ./ancombc_genotype.qza

qiime composition da-barplot \
  --i-data ./ancombc_genotype.qza \
  --p-significance-threshold 0.001 \
  --o-visualization da_barplot_genotype.qzv

qiime composition ancombc \
  --i-table ./table_2k_abund.qza \
  --m-metadata-file ./metadata.tsv \
  --p-formula 'donor * genotype' \
  --o-differentials ./ancombc_donor_genotype.qza

qiime composition da-barplot \
  --i-data ./ancombc_donor_genotype.qza \
  --p-significance-threshold 0.001 \
  --o-visualization da_barplot_donor_genotype.qzv

When you open the differential abundance bar plots generated from your ANCOM-BC results, you’ll see one bar plot per column from the group(s) included in the formula parameter in the ANCOM-BC output, excluding the selected intercept(s) pulled from the reference level parameter. If a particular group is not specified in the reference level parameter, the intercept will default to the highest alphanumeric group (e.g. in alphabetical or numeric order, as applicable) within each formula term.

Each plot visualizes features in each group compared to the intercept as LFC (log-fold change), sorted by the most relatively enriched feature to the most relatively depleted feature. As mentioned above, this visualization can be filtered by q-value (i.e., false discovery rate corrected p-value).

Taxonomic classification again¶

It is possible to increase taxonomic classification accuracy by showing the taxonomic classifier what a typical animal stool sample looks like before attempting classification. To do that we will have to retrain the naive Bayes classifier. Fortunately, a representation of a typical stool sample that is derived from Qiita data is available from the readytowear collection.

If you feel that these samples are not typical stool samples, it is possible to, for instance, assemble data on just mouse or just human (or just human and mouse) stool samples using q2-clawback. We will not attempt that here because it takes a while to run, but details are available in the tutorial.

Start by downloading the stool data, along with the 99% Greengene 13_8 reference data.

wget \
  -O "ref_seqs_v4.qza" \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/ref_seqs_v4.qza"

curl -sL \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/ref_seqs_v4.qza" > \
  "ref_seqs_v4.qza"

wget \
  -O "ref_tax.qza" \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/ref_tax.qza"

curl -sL \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/ref_tax.qza" > \
  "ref_tax.qza"

wget \
  -O "animal_distal_gut.qza" \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/animal_distal_gut.qza"

curl -sL \
  "https://data.qiime2.org/2024.2/tutorials/pd-mice/animal_distal_gut.qza" > \
  "animal_distal_gut.qza"

Next retrain the classifier.

qiime feature-classifier fit-classifier-naive-bayes \
  --i-reference-reads ./ref_seqs_v4.qza \
  --i-reference-taxonomy ./ref_tax.qza \
  --i-class-weight ./animal_distal_gut.qza \
  --o-classifier ./bespoke.qza

We can use the new classifier in exactly the same way as the standard classifier that we downloaded above.

qiime feature-classifier classify-sklearn \
  --i-reads ./dada2_rep_set.qza \
  --i-classifier ./bespoke.qza \
  --o-classification ./bespoke_taxonomy.qza

We’ll use this new taxonomy along with the original taxonomy to examine the taxonomic resolution for some of our most abundant ASVs from da_barplot_donor.qzv. Since we’ll be using multiple taxonomies, we’ll need to first put them in a directory that we’ll pass in to tabulate-seqs below as a Collection.

mkdir ./multi-taxonomy
mv ./taxonomy.qza ./multi-taxonomy
mv ./bespoke_taxonomy.qza ./multi-taxonomy

qiime feature-table tabulate-seqs \
  --i-data dada2_rep_set.qza \
  --i-taxonomy multi-taxonomy/ \
  --o-visualization dada2_rep_set_multi_taxonomy.qzv

Longitudinal analysis¶

This study includes a longitudinal component; samples from each mouse were collected 7, 14, 21, and 49 days post fecal transplant. We can use the q2-longitudinal plug-in to explore the hypothesis that a mouse’s genetic background affected the change in the microbial community of each mouse. For this longitudinal analysis, we’re going to focus on beta diversity. Alpha diversity changes wildly in infants, but it’s often stable in adults over short time periods. We’re dealing with an adult fecal community over a relatively short time period, and there is no difference in alpha diversity with time. The longitudinal analysis tutorial is an excellent resource for exploring longitudinal analyses of microbiome samples.

qiime longitudinal volatility \
  --m-metadata-file ./metadata.tsv \
  --m-metadata-file ./core-metrics-results/unweighted_unifrac_pcoa_results.qza \
  --p-state-column days_post_transplant \
  --p-individual-id-column mouse_id \
  --p-default-group-column 'donor_status' \
  --p-default-metric 'Axis 2' \
  --o-visualization ./pc_vol.qzv

qiime longitudinal first-distances \
  --i-distance-matrix ./core-metrics-results/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file ./metadata.tsv \
  --p-state-column days_post_transplant \
  --p-individual-id-column mouse_id \
  --p-baseline 7 \
  --o-first-distances ./from_first_unifrac.qza

qiime longitudinal volatility \
  --m-metadata-file ./metadata.tsv \
  --m-metadata-file ./from_first_unifrac.qza \
  --p-state-column days_post_transplant \
  --p-individual-id-column mouse_id \
  --p-default-metric Distance \
  --p-default-group-column 'donor_status' \
  --o-visualization ./from_first_unifrac_vol.qzv

qiime longitudinal linear-mixed-effects \
  --m-metadata-file ./metadata.tsv \
  --m-metadata-file ./from_first_unifrac.qza \
  --p-metric Distance \
  --p-state-column days_post_transplant \
  --p-individual-id-column mouse_id \
  --p-group-columns genotype,donor \
  --o-visualization ./from_first_unifrac_lme.qzv

Machine-learning classifiers for predicting sample characteristics¶

As an alternative (or complementary) approach to the methods we have used in this tutorial for testing if and how samples are different from one another, we can utilize machine-learning methods to determine how predictive microbiome composition is of other characteristics about a sample. For example, we may use machine-learning classifiers to predict a patient’s susceptibility to disease, or predict the treatment group that a sample belongs to. Additionally, many machine-learning methods report which features are most important for predicting sample characteristics, making this a useful approach for determining which features (ASVs, species, etc) are associated with a particular treatment, disease state, or other category of interest. All of this and much more can be found in the q2-sample-classifier plugin. Here we will use this plugin to predict each mouse’s genotype and donor status based on their ASV composition using a Random Forest classifier (this pipeline can access many different machine-learning methods via the estimator parameter, but Random Forest classifiers are used by default).

qiime sample-classifier classify-samples \
  --i-table ./dada2_table.qza \
  --m-metadata-file ./metadata.tsv \
  --m-metadata-column genotype_and_donor_status \
  --p-random-state 666 \
  --p-n-jobs 1 \
  --output-dir ./sample-classifier-results/

This pipeline generates a number of output artifacts and visualizations. You can read more about these in the sample classifier tutorial but right now let’s just focus on ./sample-classifier-results/accuracy_results.qzv. This visualization tells you how well your sample classifier performed via a confusion matrix and accompanying table of accuracy scores. This tells you how frequently each sample type is classified to each sample class, including the correct class label. Overall error rates are also reported in the table below.

Looks like we did pretty well! So we can see what features are most predictive of each sample class (donor and genotype groups). The importance scores are stored in the ./sample-classifier-results/feature_importance.qza artifact (pro tip: this can be viewed with the qiime metadata tabulate command we covered earlier). Here we will generate a heatmap showing the mean abundance of the 100 most important ASVs in each genotype and donor group.

qiime sample-classifier heatmap \
  --i-table ./dada2_table.qza \
  --i-importance ./sample-classifier-results/feature_importance.qza \
  --m-sample-metadata-file ./metadata.tsv \
  --m-sample-metadata-column genotype_and_donor_status \
  --p-group-samples \
  --p-feature-count 100 \
  --o-heatmap ./sample-classifier-results/heatmap_100-features.qzv \
  --o-filtered-table ./sample-classifier-results/filtered-table_100-features.qza

Synthesis¶

Based on the results of the analysis, we can say that there is a difference in the microbial communities of these mice based on their donor and genetic background. (This recapitulates the results of the original analysis.)

We found that the donor is the primary driver of alpha diversity.

But, we saw differences by donor and genotype based on beta diversity. Using the PCoA emperor plots, we can see clear separation between the mice from the two donors. After adjusting for the donor, we saw a significant difference between the genotypes.

Although there wasn’t a clear pattern in the barchart at the phylum level between donors or genotypes, we were still able to find ASVs which differentiated the genotypes using ANCOM-BC and Random Forest classification. There was no overlap between these ASVs in the donor and genetic background, supporting the hypothesis that the difference due to genotype is separate from the difference due to donor.

The volatility plots and temporal analysis showed that the microbiome in different genetic backgrounds changed differently over time.

This suggests that there is a genotype-specific effect on the microbiome of mice receiving fecal transplants.

💩🐁