The newest Oxford Nanopore R10.4.1 full-length 16S rRNA sequencing enables the accurate resolution of species-level microbial community profiling

In the current research work, using a custom workflow, Novaseq, Pacbio, ONT R9.4.1, and ONT R10.4.1 nanopore full-length (V1-V9) 16S ribosomal RNA-sequencing were performed on the Zymo and Synthesis communities (Fig. 1A). To evaluate the taxonomic assignments in R10.4.1, LAST and Emu for long-reads and DADA2 for Novaseq classification tools were compared using the NCBI and SILVA databases. For mock samples, we focused on error profiles, the proportion of correctly classified reads, L1 distance, true positive (TP) taxa, false positive (FP) taxa, precision, recall, and F1 scores.

For the environmental samples, we focused on the performances of R10.4.1 at the species level and the comparison of ONT and PB, as well as correlations and differences from different diversity indices. We want to know what advantages R10.4.1 has over R9. Then can R10.4.1 reach a reliable taxonomic assignment at the species level compared to PB? Globally, we showed that Oxford Nanopore R10.4.1 is capable and can be used for profiling the microbiota.

Currently, 12 species have been selected to construct a synthetic agriculture community. Considering the variation in copy numbers of the 16S rRNA gene across different bacterial genomes, we employed quantitative real-time PCR (qRT-PCR) for its quantification to obtain the abundance of copy numbers. The genomic DNA of the strains was measured in triplicate by using an SYRB Green qPCR, and the efficiency was found to be 97.57%. Each gDNA was quantified using a Qubit fluorometer. The average copy number of Serratia marcescens, Enterobacter ludwigii, Achromobacter sp., Streptomyces coelicolor, Streptomyces lividans, Sinorhizobium fredii, Bradyrhizobium diazoefficiens, Mesorhizobium jarvisii, Pseudomonas putida, Achromobacter sp. H380, Salmonella enterica, and Raoultella sp. X13 are 9.07 × 10⁷, 9.74 × 10⁷, 2.77 × 10⁷, 1.41 × 10⁷, 1.48 × 10⁷, 2.14 × 10⁷, 5.45 × 10⁶, 1.30 × 10⁷, 7.29 × 10⁷, 1.63 × 10⁷, 6.79 × 10⁷, and 9.61 × 10⁷ per ng, respectively (Table S1).

The ratio of each extracted genomic DNA in a synthetic community was determined by different proportions of mass to obtain the copy number of 16S rRNA genes (Table S1). The largest DNA value in the S1 sample was 1,000 times the smallest value, and the mass value was dispersed in a gradient. The S1 sample consists of S. marcescens (39.77%), E. ludwigii (38.44%), Achromobacter sp. (9.716%), S. coelicolor (2.467%), S. lividans (1.946%), S. fredii (1.880%), B. diazoefficiens (0.2391%), M. jarvisii (0.05683%), P. putida (1.279%), Achromobacter sp. H380 (3.574%), S. enterica (0.5955%), and Raoultella sp. X13 (0.04217%), with a 943-fold difference in 16S copy number. In the S2 sample, a 20-fold difference in genomic DNA was maintained, including S. marcescens (1.564%), E. ludwigii (1.679%), Achromobacter sp. (0.4776%), S. coelicolor (0.2426%), S. lividans (0.2551%), S. fredii (7.395%), B. diazoefficiens (1.881%), M. jarvisii (4.470%), P. putida (25.16%), Achromobacter sp. H380 (0.2811%), S. enterica (23.42%), and Raoultella sp. X13 (33.17%), with a 136-fold difference in 16S copy number. Some species have a higher similarity in 16S sequences with the same genus, such as S. lividans and S. coelicolor, which share 99.93% identity, and Achromobacter sp. and Achromobacter sp. H380, which share 98.92% similarity (Fig. 1B).

The Novaseq platform generated an average of 46,932 (SD ±5,879) high-quality pair-end reads per sample, with Q20 >96% and Q30 >90% for all the samples, respectively (Table S2). For the PB platform, 32,858 (mean) ± 2,543 (SD) clean CCS tags for environmental samples and 13,575 (mean) ± 743 (SD) for mock samples were generated (Table S3).

The ONT R9.4.1 platform displayed 42,840 (mean) ± 1,066 (SD) clean reads (Q score ≥7) per sample, with an average Q score of 12.7 (Table S4). The ONT R10.4.1 platform obtained 57,252 (mean) ± 1,309 (SD) clean reads (Q score ≥10) per sample. The average Q score for raw data for the ONT R10.4.1 platform was “18.77” (raw accuracy: 98.67%), while for data without primer, it was “22.41” (accuracy: 99.42%) (Table 1).

To evaluate the accuracy of the R10.4.1 flowcell on the amplicon after improvements, error profiling was performed on the long-read sequencing platforms (Table 2). We calculated the mismatch error rates and indel error rates of samples, meanwhile obtaining a mean error rate among all the samples.

CCS reads exhibited a mean error rate of 0.4658%, including a 0.2656% mismatch error rate, a 0.136% insertion error rate, and a 0.0636% deletion error rate, respectively. R9.4.1 reads have the highest error rate of 7.16% with 94.62% accuracy compared to raw reads. Compared to R9.4.1, the total error rate of R10.4.1 reads was only 3.16%. Moreover, R10.4.1 had a significantly lower error rate, with the mismatch error rate reduced to less than half and the insertion error rate reduced to one-third. In the case of deletion error rates, R10.4.1 performed better than R9.4.1. A deletion error rate of 0.6188% was recorded in the case of R10.4.1, whereas the deletion error rate for R9.4.1 was recorded at 2.6635%. As the synthetic community has greater sequence diversity, the error rates of ONT reads in the synthetic community are lower than those in Zymo. To demonstrate the existence of errors, we chose sequences of B. diazoefficiens from different platforms with a single copy of the 16S gene. The Integrative Genomics Viewer (IGV) results showed that R10.4.1 performed well and had lower error rates as compared to other platforms (Fig. S1).

The similarity of 16S impacts the accuracy of species identification. At least 82.18% of the reads on the PB platform had a similarity of more than 99% (sample S1, Fig. 1C). However, the similarity of reads obtained via R9.4.1 is mostly between 92% and 95%. Moreover, for R10.4.1, more than 97% of identity was accounted for by Zymo (43.89%), S1 (77.10%), and S2 (67.26%), respectively.

To quantify the results obtained from R10.4.1 in comparison to other platforms, the following matrices were used to evaluate the performance at both the genus and species levels. For the Zymo sample, PB- and ONT-based mapping methods accurately detected the eight TP species; however, Novaseq 16S was not suitable for species-level detection, which only detected three TP species (Fig. S2; Table S5). For PB data, LAST_Silva was used to obtain the least FPs with the lowest L1 distance, while Emu_NCBI possessed the second lowest FPs. For ONT data, R10.4.1 with “Emu_Silva” produced the lowest FPs. R9.4.1 was slightly better(Lowest FPs) than R10.4.1 when using the NCBI database across different methods. Using LAST, the R10.4.1 platform gets a smaller L1 distance than the R9.4.1 platform across different databases. The comparison of R10.4.1’s FP is as follows: Emu_SILVA < LAST_NCBI < Emu_NCBI < LAST_SILVA. Overall, PB outperformed ONT for the Zymo sample.

The S1 sample was an intentionally challenging community containing several microbes of the same genus but different species with similar 16S sequences (Fig. 2; Table S6). At the species level, using the LAST_NCBI method, the ONT data obtained the complete recall (12/12). However, one replicate of PB data detected 11 TPs (not detected Raoultella sp. X13, with an abundance of 0.0004). Novseq data cannot achieve species-level identification. Software, and databases contribute to the major differences in community abundance using long-reads platforms. When compared to LAST, the Emu method gets less “Other” results. Also, the number of Other found by the NCBI database is less than what the SILVA database reports. The R10.4.1 reads produced the lowest L1 distances (Emu_NCBI: 0.348), followed by R9.4.1 (Emu_NCBI: 0.364). Notably, PB produced the largest L1 distance (Emu_NCBI: 0.435). But the Emu method yields fewer FPs as well as lesser TPs (Table S6; Emu_silva for PB: lost 5 TPs; Emu_Silva for R10.4.1 and R9.4.1: lost 3 TPs; Emu_Self cannot get the full TPs). In addition, among the methods that achieve complete recall, R10.4.1 obtains the lowest L1 distance (0.525). At the species level, using the NCBI database results in a smaller L1 distance than using the Silva database. However, at the genus level, both databases exhibit similar L1 distances. When using the NCBI database, Emu achieves a smaller L1 distance at the species level, while LAST exhibits the smallest L1 distance at the genus level.

Among the expected species in sample S2, all of the long-read platforms obtained complete recall at both the genus and species levels by using LAST_NCBI. Using the above methods, PB performs better (L1 distance: 0.198 at the genus level and 0.211 at the species level) than ONT R10.4.1 (0.552 at the genus level and 0.624 at the species level). The L1 distance obtained from the NCBI database is smaller than that from the Silva database. When using the Emu method, none of the long-read platforms can achieve complete recall. In addition, PB gets far fewer FPs than ONT data (Table S7). The L1 distances of the three long-read methods analyzed using LAST_NCBI are as follows: PB < R10.4.1 < R9.4.1 (Fig. S3; Table S7). At both the species and genus levels, using the NCBI database results in a smaller L1 distance than using the Silva database. Similar to sample S1, the Emu method lost TPs more than LAST.

Correctly assigned rates from the NCBI database were higher than those from the SILVA database in the synthetic communities, but in Zymo, an opposite trend was observed (Fig. 3A and B). While using the LAST_NCBI method, the correctly assigned rates increase with increased accuracy. Interestingly, while using Emu_NCBI at species level, R9.4.1 performs better than R10.4.1 (Fig. 3B).

Using the NCBI and SILVA databases, we aligned the generated sequencing data using the Emu and LAST aligners (Fig. 3C and D). Emu successfully assigns all generated data to species and genus levels and determines the target species based on the best hit according to the algorithm. However, while using LAST, the sequences with identity less than 90% utilized the lowest common ancestor (LCA) method to determine species, and the assigned ratio of the NCBI is greater than that of the SILVA database, which grows as the accuracy improves. Particularly for synthetic communities, the assigned rate of PB at the species level was greatly decreased. These results indicate that the ability to classify sequences to the genus or species level depends on sample type, method, and database.

Mock communities can be useful for evaluating the performance of different sequencing platforms, but they lack richness. Therefore, we analyzed the microbial compositions produced by different platforms using real samples (water and soil). All samples were rarified based on the minimal number of reads (26,000 reads). As with mock, Emu successfully assigned all reads to species and genus levels. Although the rate of assigned reads by Emu was greater than that of LAST, the number of annotated species and genus was fewer. Furthermore, the number of species and genus identified in ONT R10.4.1 was higher than those identified via PB (Table 3). Next, we used the LAST aligner for further comparisons. For YW samples, 699 genera and 1,052 species were identified by PB, while 946 genera and 1,562 species were identified by R10.4.1, and 973 species and 816 genera were identified by R9.4.1. The results of read classifications were in the following order: PB > R10.4.1 > R9.4.1 for species and R10.4.1 > PB > R9.4.1 for the genus.

For NW samples, PB identified 281 genera and 405 species, while for ONT, R10.4.1 identified 566 genera and 836 species, and R9.4.1 identified the least number of species and genus. The number of species identified by R10.4.1 was much higher than that of R9.4.1. Regardless of species or genus, the rate of PB read classifications was greater than that of R10.4.1. For HS samples, R10.4.1 reads identified more species and genus than PB, although the ratio of read classifications was lower than PB. In short, the rate of R10.4.1 read classifications and the number of species and genus identified were better than those identified by R9.4.1 (Fig. S4). However, R10.4.1 obtained more species than PB.

For overall species composition results, we first compared the correlations between LAST and Emu aligners. The correlations between the two aligners were poor (Table S8). We analyzed the average abundance of the top 15 species across different platforms, aggregated them, and then calculated the correlations.

For YW samples, the average abundance of the top 15 species shared by PB and R10.4.1 was higher than that of PB and R9.4.1 at both species and genus levels (Fig. S5). Comparing the number of top 15 common species in PB and R10.4.1, the LAST results were more accurate than Emu, and the correlation between LAST and Emu was not significant. Next, we focused on the results of LAST aligners. When the LAST aligner was used, the R10.4.1 and PB had a Pearson correlation coefficient (PCC) of 0.51 (P = 0.013) for species level and 0.83 (P = 3.1e⁻⁶) for genus level. Moreover, for NW samples, the average abundance of the common top 15 species shared by PB and R10.4.1 was no less than that of PB and R9.4.1 at both species and genus levels (Fig. S6). PB and R10.4.1 have a PCC of 0.796 (P = 9.2e⁻⁶) at species levels and 0.811 (P = 1.4e⁻⁵) at genus levels, respectively. Furthermore, for HS samples, there were 10 common species and 9 genus detected by R10.4.1 and PB, but only 6 common species and 8 genus were detected by R9.4.1 and PB. PCC at the species level was 0.99 (P = 1.3e⁻¹⁹) and 0.985 (P = 6.9e⁻¹⁷) at the genus level (Fig. 4). These results suggest that the major species of R10.4.1 and PB have a similar distribution, especially in HS samples.

Boxplots of alpha diversity across three platforms are shown in Fig. 5. R10.4.1 can obtain a higher PB for Chao1 and feature. In addition, PB and R10.4.1 lie closer on the Shannon and Simpson indices, and their significance is lower than PB and R9.4.1 (Fig. 5A through C). Especially, PB and R10.4.1 do not show statistical significance (Shannon: P = 0.76; Simpon: P = 0.076) for YW samples, but there were significant differences among PB and R10.4.1.

The results of rarefaction curves showed that the Emu aligner was easier to saturate because of the identification of fewer species (Fig. 5D). By using LAST aligners, R10.4.1 is often harder to saturate. For NW samples, species saturation can be reached at 26,000 reads, but for HS and YW samples, more reads (~40,000) are required.

Finally, we compared different platforms based on the Bray–Curtis matrix. Current results suggest that species-level community composition was significantly diverse across different platforms (HS: F = 13.32, P = 0.001; NW: F = 9.80, P = 0.001; YW: F = 4.93, P = 0.007) and aligners (HS: F = 70.35, P = 0.001; NW: F = 31.16, P = 0.001; YW: F = 43.72, P = 0.001) based on a permutational multivariate analysis of variance (PERMANOVA) and visualized by Principal Coordinates Analysis (PCoA) for the unfiltered data (Fig. 6A through C). The aligners separated along the x-axis (HS: 72.97%; NW: 55.13%; YW: 66.65%), and the platforms separated along the y-axis (HS: 16.53%; NW: 24.7%; YW: 15.44%). The effect of aligners on species composition is remarkable.

In the current research work, we have for the first time compared the ONT R10.4.1 amplicon with the PB amplicon in 16S rRNA. First, we evaluate the performance of R10.4.1 on mock samples. In the synthetic community, 77.10% of the R10.4.1 reads were more than 97% identical to the reference (Fig. 1C). ONT R10.4.1 achieved a good recall at species levels (Tables S5 to S7). In addition, for environmental samples, R10.4.1 can obtain more species than PB (Table 3); meanwhile, PB and R10.4.1 are closer on the Shannon and Simpson indices (Fig. 5). Among the major species, R10.4.1 and PB platforms reveal a similar microbiome. In the current research, we tried to show that both classifiers and databases affect species composition (Fig. 6). Finally, we demonstrated that the ONT R10.4.1 flowcell for full-length 16S enabled species-level taxonomic identification for environmental samples.

The “mock sample” includes two types of communities: “commercial community” and “laboratory synthetic community.” The ZymoBIOMICSTM Microbial Community DNA Standard (https://www.zymoresearch.com/zymobiomics-community-standard, D6305) was employed as the reference sample. The mock community contains the genomic DNA of eight bacterial species and two fungi. The synthetic community contains 12 species commonly isolated in the agricultural environment, obtained from the State Key Laboratory of Agricultural Microbiology, Huazhong Agricultural University. Streptomyces coelicolor A (3)2 and S. lividans ZX1 were cultured at 28℃ in Murashig and Skoog medium. Pseudomonas putida TK2440, Achromobacter sp. H380, Serratia marcescens, Raoultella sp. X13, and Enterobacter ludwigii were cultured in Luria-Bertani (LB) liquid medium at 28℃. Salmonella enterica (ATCC14028) was cultured overnight at 37°C in LB broth. Mesorhizobium huakuii 7653R and Sinorhizobium fredii HH103 were cultured for 3 days at 28°C in a Trypticase-Yeast extract medium. Bradyrhizobium diazoefficiens USDA110 was cultured in Yeast Mannitol Agar medium.

The concentration of purified genomic DNA was initially determined using Qubit 3.0 (Thermo Scientific) according to the manufacturer’s instructions. The 341 F/805 R primers were used to amplify the 16S rRNA gene for qRT-PCR (61). Reactions were run in volumes of 20 µL, containing 10 µL of MonAmp Fast SYBR Green qPCR Mix [MonAmp SYBRR Green qPCR Mix (Low ROX), MQ10201S], 8.2 µL of nuclease-free water, 0.4 µL of forward primers (10 µM), 0.4 µL of reverse primers (10 µM), and 1 µL of DNA template. PCR amplification was carried out with the following protocol: 30 s at 95°C and 40 cycles of 15 s at 95°C, 15 s at 56°C, and 45 s at 72°C. The Bio-Rad CFX Manager software was used to perform data analysis.

The synthetic communities were mixed according to the Max/Min DNA content ratio (1:20 and 1:1,000) to construct two mock communities (named S1 and S2). The theoretical bacterial compositions are provided in Table S1.

Water samples were collected from Lake Nanhu (114.37°E, 30.48°N, Wuhan, China) “NW” and the bank of the Xinwuli Yangtze River Green Shipping Comprehensive Service area “YW” (114.25°E, 30.51°N, Wuhan, China). Three replicates are spaced at intervals of 5 m. Soil samples “HS” were collected in three replicates from the foot of Shizi Mountain (114.37°E, 30.48°N, Huazhong Agricultural University, Wuhan, China). Samples were collected in July 2022. The soils were stored in the freezer (−20°C). Water samples were filtered using sterile 0.22-μm MCE membranes and a Pall vacuum manifold (Wuhan Century Trusty Technology, HX-PV91).

DNA was extracted from 0.3 g samples using the Qiagen DNeasy Powersoil Pro Kit No. 47014 (Qiagen). Nanodrop ONE (Thermo Scientific) was used for DNA quality control, and Qubit 4.0 (Thermo Scientific) was used for quantification with the Invitrogen dsDNA HS Assay Kit. Water samples “100 mL” were mixed thoroughly and filtered using sterile 0.22-μm MCE filter membranes (Shanghai Xingya Purification Materials Factory) and a Pall vacuum manifold (Wuhan Century Trusty Technology, HX-PV91). The 16S sequences of B. diazoefficiens were verified by sanger sequencing.

For the Novaseq platform, the 16S universal primer 341 F/805 R (7) targets the V3-V4 region of the 16S gene. The PCR program was performed under the following conditions: denaturation for 3 min at 95°C; 25 cycles of denaturation at 95°C for 30 s, at 55°C for 30 s, and at 72°C for 15 s; and final extension at 72°C for 5 min. The purified DNA samples were quantified using the Qubit 4.0 fluorometer (Invitrogen, Thermo Fisher Scientific, Oregon, USA). Then, the sequencing libraries were built using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England Biolabs, MA, USA), following the manufacturer’s instructions. Finally, the library was sequenced on an Illumina NovaSeq 6000 platform at Wuhan Benagen Technology Company Limited (Wuhan, China). The 27 F/1492 R(62)targets the full-length 16S rRNA gene and was only used on the ONT and PB platforms. For the PB platform, the KOD One PCR Master Mix (TOYOBO Life Science) was used to perform 25 cycles of PCR amplification, with initial denaturation at 95°C for 2 min, followed by 25 cycles of denaturation at 98°C for 10 s, annealing at 55°C for 30 s, extension at 72°C for 1 min 30 s, and a final step at 72°C for 2 min. The PCR amplicons were purified with Agencourt AMPure XP Beads (Beckman Coulter, Indianapolis, IN, USA) and quantified using the Qubit dsDNA HS Assay Kit and Qubit 4.0 Fluorometer (Invitrogen, Thermo Fisher Scientific, Oregon, USA). SMRTbell libraries were prepared from the amplified DNA using the SMRTbell Express Template Prep Kit 2.0 according to the manufacturer’s instructions (Pacific Biosciences). Purified SMRTbell libraries from the pooled and barcoded samples were sequenced on a single PacBio Sequel II 8 M cell using the Sequel II Sequencing Kit 2.0. For the ONT platform, the first PCR protocol was 30 s at 94°C, 12 cycles of 30 s at 94°C, 20 s at 60°C, 65°C for 2 min, and 10 min at 65°C. The second PCR barcoded primer in 10 cycles. The barcoded amplicons were purified using the AMPure XP beads (Beckman Coulter, Brea, CA, USA) as per Nanopore’s instructions. Samples were then quantified using the Qubit 4.0 fluorometer (Life Technologies, Carlsbad, USA) and pooled in equimolar concentrations. A pooled DNA mixture was constructed using 1 µg of purified PCR product from each sample for library construction. The library was built using the SQK-LSK110 (for R9.4.1 flow cell) and SQK-LSK114 (for R10.4.1 flow cell) ligation kits (Oxford Nanopore Technologies, Oxford, UK) following the manufacturer’s instructions. The purified library was loaded onto R9.4.1 and R10.4.1 flowcells and then sequenced using a PromethION sequencer (Oxford Nanopore Technologies, Oxford, UK) with 48-h runs at Wuhan Benagen Technology Company Limited (Wuhan, China).

For Novaseq data, an average of 50,000–60,000 reads were generated for each of the samples in the present study. Cutadapt (v.3.5) software was used to identify and remove primer sequences and perform length filtering to obtain clean sequences.

For ONT data, the raw fast5 files were base-called using Guppy version 6.2.7 to generate fastq files. Raw reads were demultiplexed and adapters trimmed using Porechop (version 0.2.4, https://github.com/rrwick/porechop). Next, all reads shorter than 1.2 kbp and longer than 1.8 kbp were removed with Nanofilt (version 2.5.0, https://github.com/wdecoster/nanofilt).

For PB data, the raw reads generated from sequencing were filtered and demultiplexed using the SMRT Link software (version 8.0) with minPasses ≥5 and minPredictedAccuracy ≥0.9, to obtain the CCS reads. Subsequently, Lima (version 1.7.0) was employed to assign the CCS sequences to the corresponding samples based on their barcodes. CCS reads containing no primers and those reads beyond the length range (1,200–1,800 bp) were discarded through the recognition of forward and reverse primers and quality filtering using the Cutadapt quality control process (version 2.7).

The 16S rRNA gene of eight bacterial species was constructed in a database for “Zymo-specific.” Thus, 16S from the 12 species in our synthetic community created a database for “Synthesis-specific.” These two databases are named “Self-database.” Then, the NCBI 16S database and the SILVA 138 database were used. Then, the percentage identity matrix of 16S sequences was calculated using Clustal2.1 (63).

In this study, three methods were used for classification. Firstly, for the identification of feature sequences (ASVs), DADA2 (36) was used to filter raw data, splice paired-end reads, and remove chimeric sequences with QIIME2 software. Subsequently, ASVs were assigned to taxonomic groups for species annotation by using an SK-learn algorithm. This method was used for Novaseq and PB data. Secondly, the LAST program (parameters: -s 2 T 0 -a 1 -b 1 -q 1 -e 45 -m 500) was used to map clean reads to the database. For sequences whose identity exceeds 90%, the taxonomy was established based on the best matches; otherwise, the LCA of the best matching taxonomy estimation was used. Thirdly, Emu was used for taxonomy assignment with minimap2 by default. The latter two methods were used for PB data and ONT data, respectively.

Reads were used in minimap2 against the reference sample. The error profile was determined by count-errors.py (https://github.com/arq5x/nanopore-scripts/blob/master/count-errors.py), calculating mismatches, insertions, and deletions by MN and CIGAR tags (64). IGV was used for visualization (65). The identity distribution of reads was described based on the results of LAST against the Self-reference sequence. Reads were mapped to the Self database for accurate classification of each read by individual methods. Then, reads in either the SILVA or the NCBI results with Self-reference were considered to be correctly classified reads. The proportion of correctly classified reads from different classifiers and different databases was counted, respectively.

L1 distance is essentially the linear error and is calculated using the equation ${\displaystyle s\in \mathrm{S}\left|\mathrm{E}\mathrm{s}-\mathrm{I}\mathrm{s}\right|}$. where set S is the union of all species in the database and the threshold values, and Es and Is are the estimated and observed relative abundances for species “s.” “TP” represents true positives, “FP” represents false positives, and “FN” represents false negatives. Precision is defined as the percentage of TP that exists in the sample: ${\displaystyle (TP)/(TP+FP)}$. Recall shows how many true positives were picked up: ${\displaystyle (TP)/(TP+FN)}$. The F-score is calculated by taking the harmonic mean of the two values ${\displaystyle 2\cdot (\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\times \mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l})/(\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}+\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l})}$. All mock results are trimmed to include only taxa with an abundance of ≥0.01%, before the calculation of performance metrics.

PCC (r) and P-value of significance between microbial compositions obtained by different methods were calculated using the Python package “scipy.stats.pearsonr.” Heat maps were generated using the R package “pheatmap.” An α-diversity index was calculated with a vegan package (66) based on rarefied abundance counts. A t-test was used for calculating whether these samples are significantly different with a P-value of <0.05 (^*P < 0.05; ^**P < 0.01; ^***P < 0.001; ^****P < 0.0001). The Bray-Curtis dissimilarity matrix was also calculated. PCoA and PERMANOVA were used to analyze statistical differences between platforms and aligners based on the Bray-Curtis matrix. The statistical analysis and graphical visualizations were performed in R version 4.0.3.