MycorrhizaTracer

MycorrhizaTracer is a bioinformatic pipeline to automate the processing, quality control and taxonomic classification of large amounts of Sanger sequence data, optimised for fungal and plant DNA barcode sequences from the ITS region of the nuclear ribosomal genes and/or the rbcL gene. This pipeline is structured in five main steps following pre-pipeline preparations: (1) consensus calling; (2) classification; (3) clustering; (4) salvaging and (5) functional annotation (Fig. 1). Version 2.3.1 is highlighted here and all defaults and parameters mentioned here are accurate to that version, but MycorrhizaTracer is, as of publication, a living codebase undergoing active development and strict version control and, as such, the specific parameters, default values and outputs are subject to change. For the current iteration, use, parameter default values etc., please see the GitHub repository and README pages available there.

Pre-pipeline

The pipeline by default expects that each sample has been sequenced in the forward and reverse directions for some or all of the following genomic regions: full ITS (fungi), ITS2 (plants) and rbcL (plants). We have tested specific primer combinations and databases (Table 1), but the pipeline can be adapted to any taxa and genomic region by using an appropriate database (i.e. specific to that DNA region/marker) and PCR primers that target the correct genomic region in the correct taxa. For further details on how to adapt the pipeline to other taxa, please see the README on the GitHub repository which has specific instructions.

Consensus calling

Forward and reverse chromatograms are loaded and quality trimmed using PHRED quality scores (Ewing and Green 1998, Ewing et al. 1998) and chromatogram peak heights. A sliding window approach is used to remove low-quality regions. High-quality reads are aligned using Biopython’s PairwiseAligner() to produce consensus sequences. If only one read passes quality filters, that alone is used; if none of them passes the quality filters, then the sample is discarded. This supports high confidence in downstream classification. There may be times when either chromatograms are not available or when a researcher may want to manually call a consensus sequence; for these instances, the pipeline can accept as input user-defined consensus sequences in fastq format that will be used instead of the chromatogram-derived consensuses. This is implemented on a sample-by-sample basis so the pipeline can accept a mix of chromatograms and user-defined consensuses if needed.

Classification

Consensus or individual reads are classified taxonomically using BLASTn (Altschul et al. 1990). BLAST hits are filtered based on length, percent identity and coverage thresholds. A taxonomic group is assigned by finding agreement across top hits within a differential of 0.5% (by default) of the highest percentage identity. The output taxonomic rank is truncated at the last common level shared by all top hits and further constrained by minimum percentage identity cutoffs (e.g. 97% for species-level). Plant and fungal ITS sequences classified with the UNITE database (Abarenkov et al. 2023) are assigned an SH number. When databases other than UNITE are used, species-level assignments are reported as genus–species binomials. Fungal samples that do not achieve the 97% similarity species-level classification or which have not been assigned a species-level taxonomy are passed to the clustering phase for further resolution. Thus, all samples fall into one of two groups: those assigned to the species level by BLAST and those clustered into OTUs. These two types of samples are clearly demarcated in the output (the former with a genus-species binomial and the latter with a "denovo_X" designation) and can be analysed together or separately as the end user sees fit.

Clustering

Samples passed to the clustering workflow are quality-filtered and then processed through ITSx (Bengtsson-Palme et al. 2013) to retain only those with full ITS regions. These are then clustered into operational taxonomic units (OTUs) using vsearch (Rognes et al. 2016) at a default identity threshold of 97%. Each cluster is assigned a name in the format "denovo_X", representing a sequence cluster not present in the reference database. Together, these de novo clusters and the database-verified species assignments constitute the final profile of taxonomic units for each project.

Salvaging

The salvaging workflow attempts to recover sequences that failed initial classification or clustering due to sequencing errors and low sequence quality. This is achieved by ‘blasting’ the as-yet unclassified samples against the reference sequences that have already been positively identified in the study and tolerating low percentage identity of matches (assuming many sequencing errors occur in these samples). In order to help prevent artificially inflated diversity estimates, the user can define “salvage groups” that restrict the possible classifications that could be assigned to one of those already present amongst samples within the same group (often a soil core, sampling locality or geographic area). This minimises misclassification due to low-quality data and prevents inflation of local diversity estimates.

Functional Annotation

Fungal sequences with taxonomic assignments at specific taxonomic levels, such as species or genus, can be functionally annotated using FUNGuild (Nguyen et al. 2016). This optional step provides ecological context, such as trophic mode or lifestyle, aiding downstream biodiversity analysis.

Output

The primary output is a classification table detailing consensus results, taxonomic assignments, clusters and functional annotations. A log file summarises metrics describing how many samples were successfully classified, clustered, salvaged and functionally annotated. The output package also includes FASTA/FASTQ files, BLAST results, clustering output and FUNGuild results (for fungi) for transparency and manual validation.

Other

MIT licence

MycorrhizaTracer is implemented as an executable Python script (“Classify_Sanger_OTUs.py”) with a command-line interface. The pipeline is distributed with a conda environment containing all necessary packages and executables, with the one exception of FUNGuild (Nguyen et al. 2016), ensuring reproducible setup and deployment. Each major component of the workflow - input/output handling, consensus generation, taxonomic classification, clustering, salvaging and functional annotation (Fig. 1) is configurable via command-line arguments, allowing the user to fine-tune pipeline behaviour to suit specific datasets or experimental goals. The up-to-date default values of all parameters can be accessed by using the --help (or -h) flag when calling the programme.

Metadata and databases

A metadata file is required to specify the location of the chromatogram files and define salvaging groups. These groups are used during the optional salvaging step to limit taxonomic inference from low-quality sequences to predefined ecological or geographic contexts. The expected metadata format is described in Table 2.

While any .fasta file with sequences which are targeted by the primers chosen will work as a reference database, the header lines of each sequence have a specific required format. The databases for fungal ITS (--ITS_db), plant ITS2 (--ITS2_db) and plant rbcL (--RBCL_db) must have headers formatted as found in the UNITE general release .fasta (Abarenkov et al. 2023), though some of the fields are never used in the pipeline here: “Genus_species|{unused_field}|{SH_number_or_similar}|{unused_field}|k__kingdom;p__phylum;c__class;o__order;f__family;g__Genus;s__Genus_species”.

Consensus generation

The pipeline begins by loading forward and reverse chromatograms and applying quality trimming. Trimming is based on two metrics: PHRED quality scores and chromatogram peak height. These are evaluated using a sliding window approach, with user-definable parameters, such as --qual (minimum PHRED score), --window_size (window length) and --stddev_cutoff (chromatogram peak height threshold, in standard deviations below the mean). The window is moved base by base across the sequence ends, trimming low-quality regions until both metrics are satisfied. The peak height trimming is implemented specifically to deal with chromatograms that do not include quality scores and (when left at default) has minimal impact on chromatograms with PHRED scores.

After trimming, reads are assessed against a minimum length requirement, set via the --min_read_length flag. If both forward and reverse reads exceed this threshold, they are aligned using a local pairwise aligner from Biopython's Bio.Align module. A consensus sequence is then generated: matching bases are retained with the higher quality score, while mismatches are replaced with an “N”. If only one read passes filtering, it is used directly as the representative sequence for the sample and will, hereafter, be referred to as a “consensus” for simplicity and, to avoid confusion, the number of raw sequences that underlie each samples' "consensus" is reported in the output of the pipeline. Samples with no passing reads are omitted from downstream analyses.

Taxonomic classification

Taxonomic classification is performed using BLASTn (Camacho et al. 2009), with the consensus sequence or high-quality individual reads queried against a user-specified database. To refine the classification to regionally relevant taxa, the user may optionally supply species lists (via --species_list_ITS, --species_list_ITS2 or --species_list_rbcL). These lists must contain genus and species columns and only database entries with exact matches are retained during classification.

BLAST hits are initially filtered by alignment length (--minMatchLength, default 50 bp) and coverage thresholds for both query and subject sequences (--minQCOV and --minSCOV). A match passes if either coverage exceeds 0.95 or both exceed the defined thresholds. Additional filtering steps are available: the --no_Incertae_Sedis flag removes BLAST hits with ambiguous taxonomy from consideration (unless every hit is thus removed, in which case, the original list of hits is used), while --adjustPidents adjusts identity scores to refine the percentage identity of partial alignments.

The highest percentage identity from the remaining hits is saved and any hits within a user-defined identity differential (--pidentDiff, default 0.5%) of this top percentage identity are selected as candidate matches. These are compared taxonomically from kingdom to species. Where disagreement is found, the taxonomic ranking is truncated to the last level at which all hits agree. Final classification is further constrained by a series of minimum identity cutoffs for each rank (e.g. --min_pident_species, default 97.0). The final classification is stored for each sample and the top BLAST hit’s full taxonomy is also recorded separately for users requiring higher resolution who can tolerate lower accuracy, though this should be interpreted with caution.

Plant and fungal sequences that receive a species-level match against the UNITE database are assigned a UNITE SH number. For plant sequences classified using the rbcL region, a genus-species label is used instead. Fungal samples failing to achieve species-level resolution are passed to the clustering phase.

Clustering

Sequences not classified to the species level at the designated cut off (the default of 97.0 was used in this case study) enter the clustering workflow. These are subjected to stringent quality and length filtering using thresholds set by --minSeqLengthForCluster (default 300 bp) and --minAverageQualityForCluster (default PHRED 20). High-quality reads are then processed with ITSx to extract the ITS region (ITS1, 5.8S, ITS2); sequences lacking a full ITS region are saved for later (optional) salvaging.

Valid sequences are clustered using VSEARCH --cluster_fast with an identity threshold defined by --percentIdentityForCluster (default 97.0). The resulting OTUs are given names of the form denovo_X, representing species not found in the original reference database. These clusters are best understood as study-specific sequences that are not assigned to a named taxa in the reference database (which could occur for a variety of reasons including database incompleteness, unresolved taxonomy, sequence variation and/or sequencing error) and may be analysed alongside UNITE’s SH numbers for diversity analyses.

Salvaging

The optional salvaging workflow, triggered via the --salvage flag, attempts to recover information from sequences excluded during consensus or clustering. These sequences are aligned against a custom BLAST database composed of SH reference sequences and cluster centroids identified earlier in the pipeline. Each salvaged sample is assigned the taxonomy of its best hit by bitscore, provided that a sample from the same salvage group has already been assigned that taxonomy. This requirement is intended to prevent poor-quality reads from introducing artificial taxa to the dataset and helps maintain local consistency in diversity metrics.

Functional annotation

The final step, enabled with --run_funguild, performs functional annotation of fungal taxa using FUNGuild (Nguyen et al. 2016). This tool assigns trophic mode, lifestyle and ecological function, based on fungal taxonomic assignments. The parsed guild assignments and confidence scores are summarised in the output, with the full FUNGuild report also included for reference.

Benchmarking bespoke algorithms

MycorrhizaTracer uses established programmes to do many tasks (BLAST, VSEARCH, FUNGuild) and also incorporates bespoke implementations for other tasks (quality trimming, alignment and consensus generation, BLAST-ouput parsing and salvaging). The algorithms that were written specifically for the pipeline are benchmarked and validated using a dataset of 144 chromatograms obtained from 72 morphologically identified fungal fruiting bodies. These analyses can be found in Suppl. material 1 and the ground-truthed data used can be found in Suppl. material 2.

The dataset used for development and testing of this pipeline contains 2,357 individual ectomycorrhizas collected from 63 soil plots across England, for fungal and plant identification following previously published methods (Suz et al. 2014, van der Linde et al. 2018). Soil samples were collected through Defra's Natural Capital and Ecosystem Assessment (NCEA) programme. The PCR amplification used the fungal specific primer ITS1F (Gardes and Bruns 1993), the eukaryotic primer ITS4 (White et al. 1990) and two sets of primers for plants: ITS2F and ITS3R (Chen et al. 2010) targeting ITS2 and rbcLa-F (Levin et al. 2003) and rbcLa590 (Kress et al. 2009) which targets rbcL. The PCR products were visualised on an agarose gel and submitted for bidirectional Sanger sequencing to Azenta UK Ltd (GENEWIZ Sanger Sequencing Laboratory, Takeley, United Kingdom). The resulting 14,134 (a few samples were not sequenced at one of the redundant plant targets) chromatograms contain a range of sequence lengths and qualities as may be expected from environmental samples (Fig. 2). We have included poor-quality sequences here to showcase the pipeline's ability to find and appropriately handle such occurrences. The anonymised data, two metadata sheets and the commands to recreate the data below are provided in the supplement: one with site-specific salvage groups (63 groups) and one with country-specific salvage groups (one group). For the statistics reported below, version 2.3.1 was run twice, once with each metadata sheet using the flags --metadata, --prefix, --dirDepth, --salvage, --run_funguild and --sequence_stats. The databases used, defined with --ITS_db, --ITS2_db and --RBCL_db, are those referenced above.

The sample dataset and exact parameters used are available in Suppl. material 3.

Of the 2,357 samples sequenced at the three target regions, consensus sequences were successfully generated for 88% of targets on average (86% for fungal full ITS, 94% for plant ITS2 and 83% for rbcL). The majority (85%) of successful consensus sequences were constructed from both forward and reverse DNA sequences, while the remaining 15% were based on a single sequence (either forward or reverse). These numbers are primarily dependent on the sample and sequencing quality and successful consensus building is unsurprisingly far more likely with high-quality sequences (Table 3).

To be classified, a sample had to exceed the user-defined match length (default = 50 bp) and coverage thresholds (default of 0.3 for scov and qcov), meet or exceed the user-specified percentage identity threshold (default = 97% for species and 94% for genus), and show taxonomic consensus among all BLAST hits within a percetage identity differential (default = 0.5%) of the top match. For each target region, ca. 17% of samples (19.6% ITS, 15.4% ITS2, 17.7% rbcL) were not classified to any taxonomic rank due to poor sequencing quality. Of those that had some classification with the initial BLAST, 62% (fungal ITS), 54% (plant ITS2) and 47% (plant rbcL) achieved classification to the species level and 79% (ITS), 65% (ITS2) and 73% (rbcL) were classified at the genus level (Fig. 3). These taxonomic classifications are strongly dependent on the database completeness (i.e. the proportion of the species present in the study which have suitable reference sequences in the database).

Operational taxonomic units were defined as the combination of: (1) the species names of samples with species-level classifications (or SH numbers) and (2) the cluster centroids. For fungi, OTUs are a combination of species-level classifications from the BLASTn (i.e. SH numbers) and the results of clustering the remaining high-quality sequences at 97% identity. The clustering step always occurs for fungal ITS, as many fungal species are yet undescribed and do not appear in the databases, but it is optional for plants and left to the users’ discretion to include or not. The pipeline identified 677 fungal OTUs (samples = 1,796, mean size = 2.6, max. size = 117, singletons = 467), 19 plant OTUs, based on ITS2 (samples = 1,077, mean size = 56.7, max. size = 753, singletons = 8) and 15 plant OTUs based on rbcL (samples = 919, mean size = 61.2, max. size = 637, singletons = 5).

Our salvaging workflow attempts to rescue low-quality sequences by using a BLASTn search against existing OTUs with a lenient use of the highest bitscore to decide the best hit and then a stringent filter for whether another sample in the same salvage group has already been assigned to that OTU. The highest bitscore is used to allow tolerance for poor-quality sequences that may contain sequencing errors or ambiguities. The strict already-present filter is then used to keep low-quality sequences from increasing the alpha diversity of the site; low-quality sequences should not be the only instance of an OTU at a site. The number of successful salvages is, thus, mainly constrained by two factors: (1) the sequence quality impacting BLAST’s ability to identify a potential hit and (2) the size of the dataset at each site (i.e. the number and size of salvaging groups specified in the metadata). Having a finer-scale grouping structure (i.e. using study sites as groups) may result in fewer successful salvages than a coarser-scale grouping structure (i.e. using countries as groups). When using individual sites for salvage groups (n = 63 sites), 58 of the attempted 257 fungal samples were salvaged, whereas when using the country as a group (n = 1 country; all sites are in the UK), then 122 of 257 samples were successfully salvaged. The decision how fine to set the salvage grouping is an important consideration for those using the pipeline, as it will influence the balance between treating rare (or singleton) sequences as novel taxa or discarding them as low-quality sequencing attempts. Indeed, low-frequency and singleton sequences in any ecological dataset are always worth special consideration, as they are likely to be enriched for sequencing errors, but rare and undescribed species (particularly for microbes) exist and removing all singletons may result in underestimated diversity. In our study dataset, each sample is an individual ectomycorrhiza whose morphology was confirmed by microscopy, so high-quality singletons assigned as ectomycorrhizal fungi are not removed from the results.

Sanger sequencing has been used for decades and pipelines are designed with a focus on Sanger data. Table 4 presents a comparison between MycorrhizaTracer and other pipelines with similar functionality. Most similar is SangerFlow (Prodhan et al. 2025), a nextflow-based programme designed and optimised for high-throughput pest and pathogen screening which differs from MycorrhizaTracer in that it takes in FASTA files (necessitating chromatogram processing elsewhere), does not cluster the samples and does no functional annotation of fungal samples. Another similar to MycorrhizaTracer is SCATA (Durling et al. 2011), which is also optimised for fungal ITS barcoding. SCATA is only available as a web-interface and so may be more accessible in the absence of bioinformatic expertise, but it does not accept chromatograms necessitating the major bioinformatic hurdle of sequence editing and consensus building. Neither does SCATA interface with FUNGuild for functional annotation. PipeBar (Oliveira et al. 2018) and SangerAnalyzeR (Chao et al. 2021) are both designed for processing Sanger chromatograms at scale, but both stop at the consensus generation stage and, thus, classification and clustering are left to the user. Another alternative is Geneious (www.geneious.com) which has much of the same functionality (large-batch chromatogram processing, consensus generation, sequence classification with BLAST or similar, but not native clustering nor functional annotation), but is expensive and requires the user to build the pipeline steps themselves. To our knowledge, MycorrhizaTracer is the only pipeline that takes raw chromatograms through to the functional annotation (of fungal taxonomic assignments) stage. In summary, these tools fill distinct niche roles across the Sanger processing workflow rather than serving as direct functional equivalents.

Despite its strengths, MycorrhizaTracer has some limitations that warrant consideration.

Notably, many of the failures to assemble a consensus or achieve species-level classification stem from poor-quality input sequences. This underscores the importance of generating high-quality chromatograms at the sequencing stage, as downstream performance is tightly coupled to input data quality. Occasionally, using a different base-caller may improve sequence quality and could be considered in the case of low-quality sequences, but is not implemented in the current version of the pipeline.

Clustering to identify OTUs is implemented with a tunable input parameter, but only one threshold may be chosen which is applied to all taxa. It is well established that different fungal taxa diversify at different rates and so the default of 97% may not be applicable to all taxa (Nilsson et al. 2008). Cases where taxon-specific thresholds are important may require running the pipeline multiple times using different clustering values.

Taxonomic classification accuracy remains limited by the breadth and curation of underlying reference databases. Unknown or misannotated taxa can lead to ambiguous or incorrect matches (for examples, see Suppl. materials 1, 2).

The pipeline’s command-line interface, while flexible, may be less accessible to users without prior bioinformatics experience. While many parameters have default values, these are provided only for convenience and not as an endorsement that any specific default value is 'correct' for any given dataset. The reason the command-line parameters exist is so that the pipeline can be fine-tuned to accommodate the large amount of variation in real-world datasets. We encourage users to think critically about what value for each parameter is best suited for each specific use-case.

The salvage workflow, though designed to prevent overinflated diversity estimates, also depends on the completeness and correctness of sample metadata for salvage grouping and may affect relative abundance, occupancy and community structure estimates if used incautiously.

Finally, functional annotation is limited to fungal taxa interpretable by FUNGuild and confidence in guild assignment can vary widely depending on taxonomic resolution and coverage.

Critical thinking, expert opinion, visual inspection of the samples, a thorough understanding of the literature and a detailed understanding of the focal ecology are all important additions to this and any pipeline.

MycorrhizaTracer is an efficient and scalable tool for processing and classifying large numbers of DNA chromatograms. It delivers an end-to-end automated workflow, taking raw chromatograms through to functionally annotated classified taxa, while reducing the need for manual inspection of trace files. This allows for consistent, unbiased processing across large datasets. All critical thresholds (quality, coverage, identity cutoffs, taxonomic filters) are exposed as command-line arguments, allowing users to tune the pipeline for different barcoding regions, target taxa and sample quality levels without modifying the code. The optional inclusion of user-supplied species lists for regional taxa ensures that classification results are ecologically plausible and aligned with expert expectations, though such lists should be used with caution as they can inadvertently exclude unexpected, introduced, poorly documented or species whose ranges are expanding. This is particularly useful for localised biodiversity assessments. The inclusion of FUNGuild (Nguyen et al. 2016) allows ecological interpretation of the taxonomic output for fungi. This helps bridge taxonomy and function, streamlining downstream ecological analyses when sequences can be assigned taxonomy at least at genus level. The pipeline generates detailed output files, including intermediate FASTA/FASTQ, raw BLAST results, clustering summaries and full FUNGuild outputs. This transparency supports reproducibility and provides users with checkpoints for quality assurance or custom downstream workflows. In benchmark tests, 2,357 samples (14,134 chromatograms) were processed through the core programme with optional clustering and salvaging in under 20 minutes on a consumer-grade laptop (MacBook Air, Apple M3, 16 GB RAM) utilising six CPU cores, demonstrating high efficiency compared to processing chromatograms by hand, even without access to a high-performance computer.