Sequence clustering in metagenomics

In previous topics, we covered the basics of metagenomic data analysis: from data collection to processing. However, the interpretation of the results can be a challenging task due to potential errors introduced by library preparation and sequencing technology. The goal of metagenomics is to detect true biological variation between sequences, while filtering out spurious sequences. In this topic, we will discuss the most popular methods that address errors in amplicon data, and their application.

OTU and ASV at a glance

Distinguishing biological variation from amplicon sequencing errors presents unique challenges in the field of metagenomics. For accurate metagenomics datasets processing, specific error-correction methods have been developed that infer sample sequences by clustering the dataset or denoising the dataset. Depending on the algorithm, this step of analysis results in either the construction of Operational Taxonomic Units (OTUs) or Amplicon Sequence Variants (ASVs). Both approaches aim to avoid sequencing errors generating spurious taxonomic units, improve quality and reduce dataset size. Let's start with the definition and purpose of OTU:

OTU. Right after read processing, sequences can be grouped into Operational Taxonomic Units (OTUs) — clusters of sequences that share a certain degree of similarity, typically at a threshold of 97%. Clustering sequences on their identity serves several purposes: it reduces the dataset's size, lowers the computational demands of analysis, and helps in addressing sequencing errors by merging erroneous sequences with correct ones. However, in recent years, there has been a noticeable shift in the methodology employed, with a growing preference for denoising methods over traditional clustering methods, which result in Amplicon Sequence Variants (ASV).

Denoising methods aim to achieve a higher level of accuracy compared to 97% identity OTUs. The denoising approaches establish an error model based on the quality of the sequencing data, allowing them to distinguish true biological variation from errors introduced during sequencing. The result of the denoising process is the generation of ASVs, which are unique sequences even if they differ by a single nucleotide.

Denoising algorithms

There are two most popular algorithms for denoising available: UNOISE (UNOISE3), Deblur algorithm or DADA2 algorithm. We will focus on the latter, as Qiime2 wraps DADA2 (or optionally Deblur) to denoise amplicon data.

Divisive Amplicon Denoising Algorithm (DADA2) employs a statistical approach that models the error rates in sequencing reads, allowing it to identify and correct errors. DADA2 outputs ASVs, which are unique sequence variants in the dataset, aiming for a higher level of resolution compared to traditional OTUs. Let's see how this algorithm works!

How are sequencing errors estimated?

The DADA2 algorithm models sequencing errors as occurring independently within a single read and independently between different reads. It calculates the rate at which an amplicon read with sequence $i$ is produced from sample sequence $j$ based on transition probabilities between aligned nucleotides. This value will be further mentioned as error rate.

The transition probability between aligned nucleotides can depend on the original nucleotide, substituting nucleotide, and the associated quality score of the base, e.g. p(A→C, 35), between the $L$ aligned nucleotides. After sequence alignment, the error rate (λji) is computed and saved (1). If sequences $i$ and $j$ have low alignment score, $λ_{ji}$ is set to 0.

$(1)\space\space\spaceλ_{ji} = \prod_{l=0}^{L}p(j(l) \rightarrow i(l), q_i(l))$

To assess whether a specific sequence (denoted as sequence $i$) is too abundant in the dataset to be explained solely by sequencing errors, the abundance p-value ($p_A$) is used (2).

If sequencing errors are assumed to be independent across reads, the number of amplicon reads with sequence $i$ that arises from a sample sequence $j$ follows a Poisson distribution. The expectation (mean) of this distribution is the product of an error rate ($λji$) and the expected number of reads from sample sequence $j$.

$(2)\space\space\space p_A(j \rightarrow i) = \frac {1} {1 - \rho_{pois}(n_jλ_{ji}, 0)} \sum_{a=a_i}^{\infty}\rho_{pois}(n_jλ_{ji}, a)$

An example: consider a unique sequence $i$ with an abundance of $a_{i}$ found in a partition $j$, which contains $n_{j}$ reads. Given that sequence $i$ is observed at least once, the abundance p-value calculates the probability of observing $n_{j}$ or more identical reads. This probability is determined using the Poisson density function ($\rho_{pois}$).

DADA2 algorithm

The algorithm begins by grouping amplicon reads with identical sequences into unique sequences, each having an associated abundance (the number of times it appears in the dataset) and a consensus quality profile. This process is often referred to as "dereplication." Initially, all unique sequences are placed into a single partition, and the most abundant sequence within that partition is designated as the center of that partition.

After that, all unique sequences are compared to the center of their current partition and error rates for each unique sequence are calculated and stored. The sequence abundance p-value is computed for each unique sequence. After applying Bonferroni correction to the p-values, if the smallest p-value falls below a threshold, a new partition is created. The unique sequence with the smallest p-value becomes the center of this new partition. Each unique sequence is then allowed to join the partition where it is most likely to have originated (the partition that would produce the highest expected number of that unique sequence). The division procedure repeats, with each iteration involving the identification of the unique sequence with the smallest p-value, the formation of a new partition with that sequence as its center, and the rearrangement of sequences to their most likely partition.

Which approach is better?

Taken together, there are differences between OTU-based and ASV-based approaches:

• ASVs often exhibit greater sensitivity in detecting bacterial strains, although they may come at the expense of specificity. Notably, when analyzing datasets from soil, rhizosphere, and human microbiomes, both methods tend to produce similar overall biological signals.

• Nonetheless, OTUs frequently tend to overestimate bacterial diversity, especially in evaluating alpha diversity, while beta diversity estimates typically demonstrate greater consistency between the two methods.

Furthermore, the assignment of taxonomic classifications and the estimation of relative abundance for crucial taxa can vary significantly when employing OTU-based versus ASV-based software programs. So, how to be sure that we are estimating relative abundance correctly? First, remember that the reads have to be accurately processed before aggregating. Then, you can use a combination of methods in order to achieve better results. Finally, using control sample abundances, you can estimate how well the OTU or ASV calling performs.

Clustering procedures

Now that we covered how reads are assigned into OTUs or ASVs. Let's now discuss how we can find OTU or ASV groups.

1. Sequence Similarity-Based Clustering: This is one of the most common approaches. In this method, sequences are grouped based on their similarity to each other. Typically, a predefined similarity threshold is set, and sequences that share a certain percentage of similarity (e.g., 97%) are clustered together into an OTU/ASV. Tools like USEARCH tool and VSEARCH tool are commonly used for this purpose. It's important to choose an appropriate threshold as it can significantly impact the results – a higher threshold will yield fewer, more closely related OTUs, while a lower threshold may lead to over-splitting.

2. Reference-Based Clustering: In this approach, sequences are clustered using a reference database. Sequences are compared to a known database (e.g., SILVA database or Greengenes), and they are assigned to existing taxonomic groups based on their similarity. This approach is more conservative and might miss novel microorganisms, but it ensures more accurate taxonomic assignments for known species.

3. De novo Clustering: This approach involves clustering sequences without relying on a reference database. It's particularly useful when dealing with unknown or poorly characterized microorganisms. CD-HIT algorithm and Swarm algorithm are popular choices. De novo clustering can uncover novel diversity, but it might also produce some spurious OTUs, so careful filtering and quality control are essential.

Information about OTU/ASV groups may be helpful to unravel species sample structure.

Conclusion

Clustering of sequences is essential for determining of taxonomic unit, to which the read belongs, and for determination of the structure of the studied community. Modern tools utilize denoising methods over traditional clustering methods, which result in ASV instead of OTU. Denosing methods using sequencing error profile of the read to group them into ASVs.