2 DNA‐Metabarcoding data analysis

Metabarcoding

Introduction
Metabarcoding approach
Real data analysis

Introduction

The metagenomic approach is based on the analysis of total genetic material (metagenome) directly extracted from an environmental sample (Kunin et al., 2008; Wooley et al., 2010).
Although the metagenome consists of the entire genetic content of species living in a given habitat, this new approach is mainly applied to investigate microbial communities.
Given that approximately 99% of bacterial species cannot be isolated and cultured using standard laboratory techniques (Staley and Konopka, 1985), the metagenomic approach — which requires no cultivation procedure — provides an unprecedented tool for expanding our knowledge of the microbial world. The importance of accessing the "hidden" microbial fraction lies in the fact that approximately two thirds of the biodiversity on Earth is composed of microorganisms (Singh et al., 2009).
Microorganisms, their genomes, and the interactions they establish within every environment they colonize are collectively referred to as the microbiota.
The great interest in microbial communities stems from their ability to colonize virtually any habitat, from soil to ocean waters, including the human body.
Considering the human body specifically, microbiome composition varies depending on the anatomical site or the subject's age (Turnbaugh et al., 2013), as well as lifestyle.
Furthermore, the advent of high-throughput Next-Generation Sequencing (NGS) technologies (Nowrousian, 2010) has supported the study of all possible habitats.
Metagenomics can profile the microbiome through two main approaches:

Taxonomic classification: the primary objective is to identify the living organisms in a given environment (Simon and Daniel, 2011), ideally at the species level. This methodology requires, after sequencing, a classification step to assign reads to their respective taxonomic groups;
Functional characterization: the primary objective is to reveal the repertoire of genetic functions of the community under investigation through analysis of total DNA or RNA, or heterologous expression systems (usually E. coli) (Prakash and Taylor, 2012; Simon and Daniel, 2011).

As for the sequencing method, two main approaches are used:

Analysis of the entire metagenome/meta-transcriptome based on shotgun sequencing of total DNA/RNA extracted from the sample. This approach provides both taxonomic and functional information but is considerably costly in terms of sequencing and computational analysis.
The amplicon-based or metabarcoding strategy, which relies on the selective amplification of specific taxonomic markers using universal primer pairs effective across large taxonomic groups (e.g. bacteria, fungi, etc.) and the subsequent high-throughput sequencing of the resulting amplicon libraries.
This approach can only provide information on the taxonomic composition of the sample, but with superior identification capacity and species-level resolution.
It is also less expensive in terms of sequencing and computational analysis.
The most commonly used marker regions include:

several combinations of the hypervariable regions of the 16S rRNA gene (in bacteria);
the Internal Transcribed Spacers (ITS) of ribosomal RNA gene clusters (in fungi);
specific portions of the mitochondrial Cytochrome c oxidase 1 gene (in Metazoa).

Metabarcoding approach

The aim of today's practical session is the analysis of the human microbiome through the metabarcoding approach, which can ideally be divided into 5 main steps:

Sample collection;
Extraction of the metagenome (the combined genomes of [micro]organisms sharing the same environment);
Amplification of the taxonomic marker;
Sequencing;
Bioinformatic analysis of the data.

Figure 1: schematic representation of the steps characterizing the metabarcoding approach.

An ideal species genetic marker should have the following characteristics:

be ubiquitous, or at least widely shared across the taxonomic range of interest;
be sufficiently variable to allow discrimination at the deepest taxonomic levels;
be flanked by highly conserved regions suitable for universal primer design;
have a size compatible with the read-length limits of current NGS platforms.

Figure 2: representation of the variability of the 16S rRNA gene in prokaryotes.

Real data analysis

For this practical session we will use data from Gaike AH et al. 2020.
According evidence that Diabetes is strongly related to human gut composition, in this study Researcher compared gut microbiome from Type 2 Diabetes and Healthy Control and in particular:

prediabetics (PreDMs) (n = 17)
newly diagnosed diabetics (NewDMs) (n = 11), and
diabetics on antidiabetic treatment (KnownDMs) (n = 39)
compared them with healthy nondiabetics (ND) (n = 35)

The aim of this practical session is to determine whether differences exist between the gut microbiome of healthy and PD subjects.

Data were obtained by amplifying the hypervariable V4 region of the 16S rRNA and subsequently sequencing the amplicons in paired-end (PE) 2x150 mode on an Illumina MiSeq® platform.
All analyses will be performed using QIIME2.
This tool implements both a system for tracking data origin and processing, and one for their visualization, called QIIME 2 View.

Below is a graphical representation of the workflow we will apply.

Workspace setup

As a first step, access the server with the credentials assigned to you;
Then create the PD_microbiome folder and move into it:
mkdir PD_microbiome && cd PD_microbiome
Activate the virtual environment in which QIIME2 is installed:

conda activate qiime2-amplicon-2024.10

Create a folder for temporary data and instruct the operating system to use it:

mkdir tmp
export TMPDIR="$PWD/tmp"

Copy the data into your folder:
cp /data/practices_data/metabarcoding/analysis/PRJNA448494_metadata.tsv .
At this point all preliminary operations required to set up the working environment are complete.
We can proceed with the analyses.

Please Note

In case your ssh connection with virtual machine is timed out, after new ssh connection established, for this tutorial copy and paste the following commands:

cd PD_microbiome
conda activate qiime2-amplicon-2024.10
export TMPDIR="$PWD/tmp"

Data import

Import the data into QIIME2. In order to limit disk usage, you will not import the data into your own folder:
- First, generate a manifest file:

- echo sample-id forward-absolute-filepath reverse-absolute-filepath > manifest_file.tsv

- for sample in (grep -v SampleID metadata.csv|cut -f1);do echo $sample $(pwd)/input_data/
- sample* >> manifest_file.tsv ;done

- sed -i -e 's/ /\t/g' manifest_file.tsv

Then create the QIIME qza artifact:

- qiime tools import \
- --type 'SampleData[PairedEndSequencesWithQuality]' \
- --input-path manifest_file.tsv \
- --output-path pe-demux.qza \
- --input-format PairedEndFastqManifestPhred33V2

Create a visualization file (.qzv extension) to assess data quality:

qiime demux summarize \
--i-data  /data/practices_data/metabarcoding/analysis/pe-demux.qza \
--o-visualization demux-paired_end.qzv

Download the demux-paired_end.qzv file and visualize it with QIIME2 view.
- Looking at the data quality plot, how do you assess the sequencing?

Denoising

Perform Denoising of the data:
Denoising (literally, removal of background noise) is a procedure that allows identification of noise introduced by PCR and sequencing.

PE sequence trimming:

Before deciding whether to trim sequences at the 3' end and by how much, we need to evaluate the overlap between the two reads. This can be estimated using the following formula:

Where:
R is the nominal read length;
L is the average amplicon length.
In our case the amplicon has an average length of ~280 nt, therefore: S = 2*250 - 280 = 500 - 280 = 220.

Given the sequence quality and the estimated overlap region, it is advisable to limit trimming to the last positions of the reads.

- nohup qiime dada2 denoise-paired 
- --i-demultiplexed-seqs pe-demux.qza 
- --p-trunc-len-f 200
- --p-trunc-len-r 200 
- --p-trim-left-f 40 
- --p-trim-left-r 40 
- --p-max-ee-f 2 
- --p-max-ee-r 3 
- --p-n-threads 15 
- --o-table table_16S.qza 
- --o-representative-sequences rep-seqs_16S.qza 
- --o-denoising-stats denoising-stats_16S.qza &

Since this is a very time-consuming operation, denoising has already been performed. We will therefore point to precomputed qza input files.

Create visualization files to verify the outcome of the Denoising procedure:

qiime feature-table summarize \
   --i-table /data/practices_data/metabarcoding/analysis/table_16S.qza \
   --o-visualization table_16S.qzv \
   --m-sample-metadata-file PRJNA448494_metadata.tsv

qiime feature-table tabulate-seqs \
    --i-data /data/practices_data/metabarcoding/analysis/rep-seqs_16S.qza \
    --o-visualization rep-seqs.qzv

qiime metadata tabulate \
    --m-input-file /data/practices_data/metabarcoding/analysis/denoising-stats_16S.qza \
    --o-visualization denoising-stats_16S.qzv

Download the files table_16S.qzv, denoising-stats_16S.qzv and rep-seqs.qzv and visualize them with QIIME2 view.

Taxonomic classification

Perform taxonomic classification:

- nohup qiime feature-classifier classify-sklearn \
- --i-classifier /data/practices_data/metabarcoding/analysis/2024.09.backbone.v4.nb.qza \
- --i-reads /data/practices_data/metabarcoding/analysis/rep-seqs_16S.qza \
- --o-classification taxonomy_16S_SKLEARN.qza &

Since this is a very time-consuming operation, taxonomic classification has already been performed. For the following commands we will point to the precomputed qza file.

Generate qzv files to evaluate the results of taxonomic classification:

qiime metadata tabulate \
    --m-input-file /data/practices_data/metabarcoding/analysis/taxonomy_16S_SKLEARN.qza \
    --o-visualization taxonomy_16S_SKLEARN.qzv

qiime taxa barplot \
    --i-table /data/practices_data/metabarcoding/analysis/table_16S.qza \
    --i-taxonomy /data/practices_data/metabarcoding/analysis/taxonomy_16S_SKLEARN.qza \
    --m-metadata-file PRJNA448494_metadata.tsv \
    --o-visualization taxa-bar-plots_16S_SKLEARN.qzv

Download the files taxonomy_16S_SKLEARN.qzv and taxa-bar-plots_16S_SKLEARN.qzv and visualize them with QIIME2 view.

Alpha and Beta Diversity

Diversity is a measure of the complexity of the system under observation, as a function of the number of objects and their relative abundance.

In particular:

alpha: intra-sample diversity. Several metrics can be used:
- Richness: represents the number of different species observed in a biological community;
- Evenness: indicates how uniform the observed community is;
- Shannon Index: a quantitative measure of species richness;
- Faith's Phylogenetic Diversity: a qualitative measure of community richness that accounts for phylogenetic relationships.
beta: inter-sample diversity: we use measures that indicate how different two communities are. These measures compare samples pairwise, ultimately producing a matrix indicating the dissimilarity between each pair of samples. Several types of metrics exist for calculating Beta diversity:
- Jaccard distance: a qualitative measure of dissimilarity;
  
  X and Y represent the samples under analysis, specifically the observed species. This measure therefore only accounts for whether a species is observed or not, not its abundance.
- Bray-Curtis distance: a quantitative measure of dissimilarity;
  
  Where:
i and j are the two samples under analysis;
Si and Sj represent the sum of the species counts observed in sites i and j;
Cij: is the sum of the lower counts of species common to both sites.
- Unweighted UniFrac distance: a qualitative measure of dissimilarity incorporating phylogenetic relationships among ASVs;
- Weighted UniFrac distance: a quantitative measure of dissimilarity incorporating phylogenetic relationships among ASVs;
gamma: total diversity.
Before proceeding with diversity measures, we need to build a phylogenetic tree from the obtained ASVs.

- qiime phylogeny align-to-tree-mafft-fasttree \
-  --i-sequences /data/practices_data/metabarcoding/analysis/rep-seqs_16S.qza \
-  --o-alignment aligned-rep-seqs_16S.qza \
-  --o-masked-alignment masked-aligned-rep-seqs_16S.qza \
-  --o-tree unrooted-tree_16S.qza \
-  --p-n-threads 1 \
-  --o-rooted-tree rooted-tree_16S.qza

Let's copy the output of the previous command in your working directory

cp /data/practices_data/metabarcoding/analysis/aligned-rep-seqs_16S.qza .
cp /data/practices_data/metabarcoding/analysis/masked-aligned-rep-seqs_16S.qza .
cp /data/practices_data/metabarcoding/analysis/unrooted-tree_16S.qza .
cp /data/practices_data/metabarcoding/analysis/rooted-tree_16S.qza .

In order to compare the obtained data, we need to normalize them. The normalization procedure used for microbiome data is Rarefaction.
This procedure is based on choosing a rarefaction value (sampling depth) to which all samples will be subsampled. In practice, sampling without replacement is performed until the chosen rarefaction value is reached.
All samples with counts below the sampling depth will be discarded.

qiime diversity alpha-rarefaction \
    --i-table /data/practices_data/metabarcoding/analysis/table_16S.qza \
    --i-phylogeny /data/practices_data/metabarcoding/analysis/rooted-tree_16S.qza \
    --p-max-depth 100000 \
    --p-min-depth 1000 \
    --m-metadata-file PRJNA448494_metadata.tsv \
    --o-visualization alpha-rarefaction.qzv

Proceed with the estimation of diversity indices:

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny /data/practices_data/metabarcoding/analysis/rooted-tree_16S.qza \
  --i-table /data/practices_data/metabarcoding/analysis/table_16S.qza \
  --p-sampling-depth 50000 \
  --m-metadata-file PRJNA448494_metadata.tsv \
  --p-n-jobs-or-threads 1 \
  --output-dir core-metrics-results_16S

Check whether there are significant differences between alpha diversity indices:

qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results_16S/shannon_vector.qza \
  --m-metadata-file PRJNA448494_metadata.tsv \
  --o-visualization core-metrics-results_16S/shannon-Condition-significance_16S.qzv

qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results_16S/faith_pd_vector.qza \
  --m-metadata-file PRJNA448494_metadata.tsv \
  --o-visualization core-metrics-results_16S/faith_pd-Condition-significance.qzv

Check whether there are significant differences between beta diversity indices:

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

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

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results_16S/bray_curtis_distance_matrix.qza \
  --m-metadata-file PRJNA448494_metadata.tsv \
  --m-metadata-column status \
  --o-visualization core-metrics-results_16S/bray_curtis-Condition-significance.qzv \
  --p-pairwise

Comparative analysis between conditions

We will use ANCOM-BC (analysis of composition of microbiomes with bias-correction) to perform a comparative analysis between the two conditions.

Perform the comparative analysis:

qiime composition ancombc \
--i-table /data/practices_data/metabarcoding/analysis/table_16S.qza \
--m-metadata-file PRJNA448494_metadata.tsv \
--p-formula 'Diet + sex + status' \
--p-reference-levels status::Control \
--o-differentials dataloaf.qza

Create tables to visualize the results:

qiime composition tabulate \
--i-data dataloaf.qza \
--o-visualization T2D_ASV_condition.qzv

Create barplots for statistically relevant taxa:

qiime composition da-barplot \
--i-data dataloaf.qza \
--p-significance-threshold 0.05 \
--o-visualization T2D_ASV_status_barplot.qzv

The comparison was performed at the ASV level, but we may be interested in comparing at higher taxonomic levels.
4. Group the data at the taxonomic level of interest. In this case, level 6 of the GreenGenes2 taxonomy (genus):

qiime taxa collapse \
--i-table /data/practices_data/metabarcoding/analysis/table_16S.qza \
--i-taxonomy /data/practices_data/metabarcoding/analysis/taxonomy_16S_SKLEARN.qza \
--p-level 6 \
--o-collapsed-table table-l6.qza

Perform the comparative analysis:

qiime composition ancombc \
--i-table table-l6.qza \
--m-metadata-file PRJNA448494_metadata.tsv \
--p-formula 'Diet + sex + status' \
--p-reference-levels status::Control \
--o-differentials data_genus.qza

Visualize the results:

qiime composition tabulate \
--i-data data_genus.qza \
--o-visualization T2D_Genus_status.qzv

qiime composition da-barplot \
--i-data data_genus.qza \
--p-significance-threshold 0.05 \
--o-visualization T2D_Genus_status_barplot.qzv

If we wanted to perform the same analysis at the family level, what would we need to do?

Navigation Menu