Step 1: MacSyFinder Integration

MacSyFinder is used to search for gene clusters representing satellite elements. Genome data in multi-protein FASTA format is scanned against predefined models developed by Institute Pasteur, which users cannot modify. These models specify gene presence, forbidden genes (like structural phage proteins), and the maximum distance between genes to identify satellites.

MacSyFinder uses HMMER to match input sequences to Hidden Markov Model profiles of satellite-associated genes. It then forms gene clusters based on proximity and evaluates them based on quorum requirements to select the most relevant satellite systems (Abby et al., 2014).

Step 2: Post-Processing Script

The post-processing script fixes two issues:

1. Rescuing discarded clusters: Clusters excluded due to the presence of forbidden genes are reevaluated and rescued if they contain only one forbidden protein.
2. Splitting merged satellites: Satellites that are too close together may be mistakenly merged. The script uses integrases as markers to split these clusters correctly.

Prior to the development of SaPhARI, we explored various metagenomic classification tools to identify phages in environmental samples. Our aim was to adapt these methods for satellite phage detection after constructing our novel database of characterized satellites. Early on, we encountered limitations in utilizing metagenomic data due to the dependence of modern classification software on reference databases. Nevertheless, we sought to assess the sensitivity of these tools in detecting viral sequences, in the hopes of eventually leveraging them for satellite identification.

To benchmark different taxonomic classification and profiling software, we used metagenomic data from the soil microcosm experiments conducted by William & Mary iGEM 2023. This dataset was selected due to its well-documented experimental design and expected outcomes, such as increased populations of M. smegmatis and Mycobacterium phage Kampy. By using these tools to detect Mycobacterium phage Kampy, we were able to evaluate their potential utility for satellite detection once a robust satellite database had been established.

Before selecting specific tools, we conducted a review of the latest metagenomic taxonomic identification software. Two particular papers, titled the Critical Assessment of Metagenomic Interpretation, composed of comprehensive benchmark studies from 2017 and 2022 which assessed the precision and recall of various metagenomic classification tools (Meyer et al., 2022; Sczyrba et al., 2017). In addition, we consulted with faculty experts at William & Mary, who recommended the use of k-mer-based alignment software due to its efficiency in reducing computational overhead compared to traditional nucleotide alignment strategies, such as BLAST. K-mer alignment operates by analyzing all possible combinations of nucleotide sequences of a fixed length (k-mers), offering a computationally efficient solution for large-scale metagenomic analyses. 

Tools Evaluated:

Kraken2/Bracken:

Kraken2 assigns taxonomic labels by comparing k-mers from sequences to a reference database, while Bracken refines these classifications by estimating the relative abundance of each taxon. The key difference lies in their focus: Kraken2 is primarily concerned with identifying taxa, while Bracken quantifies their abundance within the dataset. Notably, Kraken2’s confidence thresholds significantly influence its output. Higher precision settings tended to overestimate the percentage of Mycobacterium smegmatis within our iGEM 2023 soil metagenomes, while lower or default confidence levels improved recall but resulted in lower abundance estimates. Despite these adjustments, Kraken2 still struggled to classify the majority of reads, with 70% remaining unclassified at default confidence levels, and up to 99% unclassified at higher confidence settings. Moreover, Kraken2 typically classified sequences only to the lowest taxonomic level it could confidently identify, which was rarely species-specific (Lu et al., 2017; Wood et al., 2019).

Additionally, Kraken2’s memory usage increases substantially with larger databases. Although it is marketed as a compact tool, using the full NCBI Nucleotide database required 1TB of RAM to fit the hash table structure required for classification, which we fortunately had access to via our high-performance computing (HPC) system. We also explored memory-mapping the database to avoid fully loading it onto the computing node, but this dramatically slowed down classification, extending processing times to several days. Ultimately, we had to reduce our database to the RefSeq Archaea, Bacterial, and Viral dataset, which considerably lowered computational demands.

Kraken2 also had difficulty distinguishing phages due to their genetic similarities and the limited representation of phage sequences within our reads. For example, the Kampy phage appeared in very low numbers, while the genus Backyardiganvirus had a 100-fold higher representation. Interestingly, when we created a custom database containing only the Kampy phage, the classification numbers for both Kampy and Backyardiganvirus were almost identical. This suggests Kraken2 struggles to differentiate between species within the same genus, particularly with closely related phages. These limitations indicate that Kraken2 may face significant challenges in distinguishing small, highly similar satellite phage sequences.

Metalign:

Metalign utilizes a two-step algorithm, beginning with a pre-filtering step based on the Jaccard index, which measures the similarity between the k-mers in the query sequences and those in the reference database. This pre-filtering process selects a subset of the database that shares significant k-mer overlap with the sample, thereby reducing the search space and computational load. However, one drawback is that Metalign discarded 80% of the available k-mers in our sequencing data, limiting the amount of genomic information processed. Additionally, Metalign requires that read FASTQs be merged with external software such as Flash or BBmerge before running, which can exacerbate the loss of reads during processing (Bushell et al., 2017; LaPierre et al., 2020; Magoč et al., 2011). 

After establishing the pre-filtered database, Metalign aligns the reads to the filtered reference genomes and estimates the relative abundance of each detected taxon. As a taxonomic profiler rather than a classifier, Metalign focuses on providing an overview of microbial community composition and relative abundances, rather than precisely identifying the organism associated with each read within a dataset. While this broad profiling offered insights into community-level dynamics, Metalign struggled with organisms not represented in its reference database, which could not be updated by the user (Liu et al., 2022).

Phage detection also proved challenging for Metalign, as phages tend to exhibit high genetic variability, making it difficult to accurately profile them using short reads. In our samples, Metalign identified only a single Mycobacterium phage, which was not the expected phage Kampy from our experimental data. Moreover, since Metalign's database cannot be updated, it lacks the flexibility to identify newly discovered satellites, making it less suitable for satellite identification.

Centrifuge:

Centrifuge, like Kraken2, uses a k-mer-based approach for taxonomic classification of metagenomic reads, allowing for rapid assignment of taxa based on reference sequences. Both tools construct an index from a reference database, enabling quick k-mer matching to identify taxa associated with the given reads. These algorithms are designed to achieve high efficiency and speed, outperforming traditional nucleotide alignment tools like BLAST, by leveraging pre-built databases to reduce computational load during classification. (Kim et al., 2016)

While Centrifuge provided results comparable to Kraken2, it lacked the ability to adjust confidence levels, limiting user flexibility. Additionally, Centrifuge required significantly more computational resources and time to index the database compared to Kraken2. Installation through a conda environment took several hours, and indexing of the RefSeq Archaea, Bacterial, and Viral database took several days. Although both tools offered similar outcomes, Kraken2's adjustable confidence intervals and more efficient database construction process give it an advantage as a taxonomic classifier.

mOTUS:

mOTUS is a marker gene analysis profiler that relies on an internal database of species-specific indicative genes. Although tested for its ability to detect Mycobacterium smegmatis, it was ineffective, classifying less than 1% of the sequencing reads. Overall, this classification tool proved inadequate for detecting satellite phages in metagenomic data, and even struggled to classify any bacteria (Ruscheweyh et al., 2021).

Lessons Learned:

The performance of taxonomic classifiers is highly dependent on the confidence thresholds set during analysis. While these tools can help detect specific constructs in metagenomic reads, achieving precise results is challenging—especially when a large portion of reads remain unclassified. Additionally, relying on short genomic reads can be limiting, as they often lack distinct, informative regions. Another key limitation is that these classifiers can only identify satellite phages present in the reference database, restricting their ability to detect novel sequences. Ideally, these tools would allow us to identify satellite constructs within metagenomic data, as we have a clear target in mind for our analysis.

Though taxonomic classifiers like Kraken2 and Centrifuge are useful for estimating the microbial species present in a sample, they often fall short in reliably quantifying species-level data, especially for small or uncharacterized constructs like satellite phages. These tools could be more effective in identifying engineered constructs if we focused on assembling contigs rather than relying on sequencing reads, a more reliable approach for satellite phage discovery through metagenomics. 

These challenges highlight the need for specialized tools like SaPhARI. Unlike traditional phage identification tools that are confined to known reference databases, SaPhARI allows for the discovery of novel satellite phages by going beyond the constructs of existing knowledge. After beginning our development SaPhARI, we shifted our metagenomics strategy to complement this novel tool, better suited for identifying satellites.

Assembly in metagenomics involves using short DNA reads obtained from sequencing to reconstruct longer, contiguous sequences known as contigs. These contigs are formed by aligning overlapping sequences, providing a clearer and more comprehensive representation of the metagenome. A contig is a continuous stretch of DNA pieced together from shorter fragments, which alone are often too brief to represent entire genes or genomic regions. Assembling these reads into contigs enables researchers to analyze the genetic material of microbial communities more effectively (Wang et al., 2019).

One common approach to assembly is through the use of de Bruijn graphs, which are mathematical structures that map shared k-mers—short iterations of nucleotide sequences of a specified length—across different reads. In a de Bruijn graph, each distinct k-mer is represented as a node, and edges connect nodes that share overlapping k-mer sequences. This method allows researchers to trace paths that represent potential contiguous sequences, facilitating genome reconstruction from fragmented reads (Compeau et al., 2011).However, building and analyzing de Bruijn graphs is computationally intensive, as it requires processing vast amounts of sequencing data and managing complex overlaps between k-mers in large datasets. 

As mentioned earlier, existing metagenomic classification tools had significant limitations when it came to satellite phage discovery, particularly due to their reliance on reference-based databases and inability to classify at higher taxonomic orders based off of short read sequencing data, which restricted the identification of novel satellite sequences. These constraints underscored the need for a de novo assembly approach, focuses on reconstructing raw metagenomic data to identify novel satellite phages and other related entities using our custom software, SaPhARI. 

However, while assembly holds promise, it is a complex and error-prone process. A key challenge in metagenomic assembly is the risk of misassemblies, which can lead to uncertainty about the accuracy of the resulting contigs. Metagenomic samples often contain highly conserved sequences shared across species, making it possible to misclassify sequences from one species as belonging to another. Additionally, repetitive or overlapping sequences may result in the assembly of artificially longer constructs than what is actually reflected in the sample, creating a misleading picture of the genomic content. The complexity of microbial communities further compounds this issue, as low-abundance species may not be fully represented, leading to incomplete assemblies. Moreover, metagenomic assemblies rarely result in full genomes; instead, researchers often obtain fragmented sequences that do not overlap well. As such, the reliability of assembled contigs must be carefully evaluated, as the potential for misassemblies introduces a degree of uncertainty into the analysis.(Wang et al., 2019). Advances in bioinformatic software seek to minimize these potential errors, and the field of metagenomics continues to improve the quality of assemblies with the development of novel algorithms and error-correcting capabilities. While these software have their limitations, they can still provide crucial insights into microbial life in diverse environments and allow for the discovery and documentation of novel species (Meyer et al., 2022; Sczyrba et al., 2017).

Literature Review:

We evaluated several short-read metagenomic assemblers, applying the same benchmarking studies and algorithms previously used for taxonomic classification software. Our goal was to determine the most effective tool for assembling metagenomes to be used with SaPhARI. Based on our analysis, we arrived at the following conclusions:

MEGAHIT:

MEGAHIT is a highly efficient assembler specifically designed for metagenomic data. It utilizes succinct de Bruijn graphs, a memory-optimized variant of traditional de Bruijn graphs that represent overlapping k-mers between reads, enabling the assembly of short reads into longer contiguous sequences. By using this compressed structure, MEGAHIT efficiently handles large datasets with lower memory requirements, making it well-suited for high-throughput metagenomic sequencing. Its key advantages are speed and scalability, making it ideal for large-scale metagenomic projects, although the quality of the resulting assemblies can vary in quality compared to other assemblers (Liu et al., 2015).

metaSPAdes:

metaSPAdes is a more resource-intensive assembler than MEGAHIT, though slower, specifically tailored for metagenomic data and built on the widely-used SPAdes algorithm for single-genome assembly. Like MEGAHIT, metaSPAdes employs de Bruijn graphs, but it enhances the assembly process by incorporating advanced strategies such as read error correction and paired-end read utilization to improve both accuracy and assembly quality. While it requires more computational resources and time, metaSPAdes often delivers comparable or superior assembly quality, particularly for complex or low-abundance microbial communities. Its ability to handle intricate metagenomic datasets with greater precision makes it a strong choice for high-quality assemblies (Bankevich et al., 2012; Meyer et al., 2022; Nurk et al., 2012; Sczyrba et al., 2017).

Velvet:

Velvet, one of the earlier assemblers utilizing de Bruijn graphs, has largely fallen out of favor for metagenomic assembly due to its limitations. While effective for assembling single-genome data, it struggles with the complexity and variability inherent in metagenomic datasets, often resulting in less effective assemblies. Velvet also demands more computational resources and produces lower-quality results compared to modern assemblers like MEGAHIT and metaSPAdes, making it a less viable option for contemporary metagenomic projects (Terrón-Camero et al., 2022).

After reviewing the literature, MEGAHIT and metaSPAdes emerged as the most reliable assemblers for metagenomic data. MEGAHIT consistently outperformed in terms of speed and memory efficiency, while metaSPAdes demonstrated better assembly accuracy and genome fraction coverage, although MEGAHIT occasionally showed a slight advantage in genome fraction during benchmark studies. Ultimately, we chose to use MEGAHIT due to its scalability for large metagenomic projects and comparable quality, aligning with our goal of analyzing as many metagenomes as possible for satellite phage identification (Kumar et al., 2023; Wang et al., 2020).

Evaluating Assemblies:

metaQUAST is a specialized quality assessment tool designed for evaluating the accuracy and completeness of metagenomic assemblies. Building on the functionality of QUAST (Quality Assessment Tool for Genome Assemblies), metaQUAST includes features tailored to address the complexities of metagenomic data. It operates by mapping the original reads back to the assembled sequences, enabling the assessment of key metrics such as genome coverage, the number of misassemblies, contig lengths, and the presence of fragmented or missing regions. (Gurevich, 2013; Mikheenko, 2016).

By using metaQUAST, we were able to map our reads back to the assembled contigs to assess both the coverage and quality of our assemblies. A high proportion of reads mapping back to the contigs indicates a well-supported assembly with minimal gaps or errors. After running metaQUAST on the assemblies produced by MEGAHIT, we observed excellent read mapping coverage, with the vast majority of reads aligning back to the contigs. These results validated the high quality of the MEGAHIT assembly and confirmed that it was an optimal choice for our project.