Metagenome species level classification method based on Debrueine diagram

Abstract

The invention discloses a metagenome species level classification method based on a Debrueine graph, which takes a linear genome and a generic genome as a complementary reference framework, and adopts a two-stage process: firstly, executing approximate member query by utilizing HIBF, and screening out candidate species and strains from a database in a high-recall and low-overhead manner; and then ccDBG is constructed around the candidates, and the attribution of each read length is finely judged on the strain scale. The method has the rapid screening capability of linear reference and high-resolution capture of a generic genome on intra-population differences, so that the problems of fuzzy and confusion of species classification in a traditional process are remarkably relieved.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of macro-genome classification analysis, and more particularly relates to a macro-genome species level classification method based on DeBruijn graphs (abbreviated as DBG). BACKGROUND

[0002] In macro-genomics, species level classification is still the primary goal of most studies: on the one hand, species composition determines the community diversity pattern and ecological function; on the other hand, it is also the premise of strain analysis and functional annotation. If the species level determination deviates, the subsequent inference will lose its meaning.

[0003] The existing macro-genome species level classification methods include three kinds: one is a marker gene-based method, such as MetaPhlAn4, mOTUs3, the core principle of which is to use a set of carefully selected, highly conserved at a specific classification level and with sufficient resolving power single-copy genes as markers. It performs species identification and abundance quantification by aligning sequencing data to these marker gene databases; one is a DNA-to-protein method, such as Kaiju, DIAMOND, the core principle of which is to translate sequencing reads into possible amino acid sequences under the reading frame, and then perform fast alignment with protein reference databases; one is a DNA-to-DNA method, the core principle of which is to directly align or match sequencing reads with genomic DNA reference databases. The DNA-to-DNA method can be further divided into k-mer-based tools, such as Kraken2, Bracken, Centrifuge, Ganon2, KMCP, Sylph, and mapping-based tools, such as MetaMaps, MORA.

[0004] However, the above several existing macro-genome species level classification methods all have some defects that cannot be ignored: First, the above three methods all have the defect that the recall and precision of species classification results are difficult to balance; Second, the classification performance of the above three methods is affected by sequencing depth. With the increase of sequencing depth, the classification performance decreases; or with the decrease of sequencing depth, the classification performance decreases. SUMMARY

[0005] In view of the above defects or improvement needs of the prior art, the present application provides a macro-genome species level classification method based on DeBruijn graphs, which aims to solve the technical problems of the existing macro-genome species level classification methods that the recall and precision of species classification results are difficult to balance, and the classification performance is affected by sequencing depth.

[0006] To achieve the above objectives, according to one aspect of the present invention, a metagenomic species-level classification method based on de Bruin diagrams is provided, comprising the following steps: (1) Obtain the complete genome, and construct the HIBF index of the reference database based on the complete genome and using the hierarchical staggered Bloom filter HIBF; (2) Obtain the metagenomic sequencing data input by the user, and based on the HIBF index of the reference database obtained in step (1), obtain the strain corresponding to each read length in the metagenomic sequencing data, and perform false positive filtering on all the strains corresponding to the obtained read lengths, so as to obtain multiple initial screening species, the strains corresponding to each initial screening species, and the read lengths corresponding to the strains. (3) Use the expectation-maximization EM algorithm to obtain the relative abundance between each primary screening species and its corresponding primary screening strain obtained in step (2), and perform threshold filtering on all the obtained relative abundance to obtain multiple candidate species and multiple candidate strains. (4) Based on the multiple candidate strains obtained in step (3) and the complete genomes of multiple strains obtained in step (1), a color and compressed de Bruin diagram ccDBG is constructed using the GGCAT tool. The ccDBG diagram records multiple k-mers and the strains to which each k-mer belongs, where the value of k is 19 or higher. (5) Obtain the metagenomic sequencing data input by the user, use the GGCAT tool to split each read in the metagenomic sequencing data into k-mers to obtain the set of k-mers corresponding to the read, and obtain multiple k-mers that match the set of k-mers corresponding to the read in the ccDBG graph constructed in step (4), as well as the strains to which each k-mer belongs. All strains corresponding to the read in the ccDBG graph are the strains associated with the read. Finally, organize the above information into a read-strain mapping table. Specifically, establish a read-strain mapping table based on the number of all reads in the metagenomic sequencing data and the number of strains associated with them. Each row of the read-strain mapping table is a different read in the metagenomic sequencing data, and the columns are the strains associated with each read in the ccDBG. (6) Obtain the final metagenomic species based on the read length-strain mapping table obtained in step (5); (7) Obtain the first result obtained in step (6) Relative abundance of the final metagenomic species As a classification result, among ∈[1, the total number of all final metagenomic species obtained in step (7)].

[0007] Preferably, step (1) specifically comprises the following steps: firstly, obtaining a complete genome; secondly, obtaining all strains and species in the complete genome; thirdly, establishing a database information table according to all strains and species in the complete genome, wherein the database information table comprises five entries, i.e., the name of each strain, the classification number of each strain, the name of each species, the classification number of each species, and the path of the obtained complete genome; and finally, constructing a HIBF index of a reference database based on the established database information table and the complete genome and using a hierarchical interlaced Bloom filter (HIBF).

[0008] Preferably, step (2) comprises the following sub-steps: (2-1) obtaining the user-input metagenomic sequencing data and performing decomposition processing on each read in the metagenomic sequencing data using a k-mer method to obtain a k-mer set corresponding to the read; (2-2) inputting the k-mer set corresponding to each read obtained in step (2-1) into the HIBF index of the reference database constructed in step (1) to perform an approximate membership query to obtain a mapping relationship between the read and each strain in the HIBF index of the reference database and a k-mer matching result; (2-3) obtaining the length L_i of the i-th read, the length k of the k-mer set corresponding to the i-th read, and the total number of k-mers in the i-th read , and obtaining the score of the j-th strain in the HIBF index of the reference database on the i-th read ; wherein i∈[1, the total number of reads in the metagenomic sequencing data], j∈[1, the total number of strains in the HIBF index of the reference database], represents the number of k-mer matching results between the i-th read and the j-th strain; (2-4) for the i-th read, retaining all strains in the HIBF index of the reference database that have a mapping relationship with the i-th read and have a score exceeding 0.75 as high-score strains, obtaining the maximum value among the scores of all high-score strains, and deleting the high-score strains having a score less than or equal to the maximum value*0.9 to obtain co-leading strains, using all co-leading strains as primary screening strains, and using the species corresponding to all co-leading strains as primary screening species.

[0009] Preferably, step (3) comprises the following sub-steps: (3-1) using the primary screening species obtained in step (2), the corresponding primary screening strains and reads thereof, and the EM algorithm to solve and perform normalization processing to obtain the relative abundance between each primary screening species and the corresponding primary screening strains thereof; (3-2) Based on the relative abundance of each primary screening species and its corresponding primary screening strain obtained in step (3-1), all primary screening species are filtered, i.e., those with a relative abundance greater than 10 are retained. -7 The initial screening species were used as candidate species; when the total number of candidate species exceeded 1000, the three strains with the highest relative abundance of all candidate species were retained as candidate strains; when the number of candidate species was less than or equal to 1000, the strains with a relative abundance greater than 10 were retained. -7 All the initial screening strains were used as candidate strains.

[0010] Preferably, step (6) includes the following sub-steps: (6-1) Obtain the maximum value of each row from the read length-strain mapping table obtained in step (5), and set the other values ​​in the row to 0, keeping only the maximum value. If there are multiple maximum values ​​in a row, keep all the maximum values, thus obtaining the read length-strain mapping table after the first round of processing. (6-2) Divide all values ​​in the read-strain mapping table obtained in step (6-1) after the first round of processing by... To obtain the updated read-strain mapping table; where This indicates the length of each read, expressed in bases (bp). (6-3) Take the maximum value of each row in the updated read-strain mapping table obtained in step (6-2) as the score corresponding to that row, and calculate the 30th percentile of the scores corresponding to all rows. Then the 30th percentile Cut off to interval Within this range, discard entries in the updated read-strain mapping table that score less than or equal to the 30th percentile. All rows are processed to obtain the read-strain mapping table after the second round of processing; (6-4) Delete all rows in the read length-strain mapping table obtained in the second round of processing obtained in step (6-3) that do not meet the first preset condition, wherein the first preset condition is that all strains in the row correspond to the same species in the database information table obtained in step (1), so as to obtain the read length-strain mapping table after the third round of processing. (6-5) Determine whether the metagenomic sequencing data obtained in step (5) is paired-end sequencing data or single-paired-end sequencing data. If it is paired-end sequencing data, first delete all rows in the read length-strain mapping table after the third round of processing obtained in step (6-4) that do not meet the second preset condition. The second preset condition is that all strains in the row correspond to the same species in the database information table obtained in step (1). Then delete the columns and rows with all zeros in the read length-strain mapping table after the third round of processing to obtain the read length-strain mapping table after the fourth round of processing. Then proceed to step (6-6). If it is single-paired-end sequencing data, delete the columns and rows with all zeros in the read length-strain mapping table after the third round of processing to obtain the read length-strain mapping table after the fourth round of processing. Then proceed to step (6-6). (6-6) Set all the remaining strains in the read length-strain mapping table after the fourth round of processing as candidate species corresponding to all species in the database information table obtained in step (1), and obtain the number of supported read lengths for each candidate species; (6-7) Based on the number of supporting reads for each candidate species obtained in step (6-6), the read-strain mapping table obtained in the fourth round of processing in step (6-5) is further filtered to obtain the metagenomic species classification results of the metagenomic sequencing data.

[0011] Preferably, the number of supported reads for a candidate species is equal to the number of rows with values ​​greater than zero in the corresponding column of the read-strain mapping table after the fourth round of processing, corresponding to the strain of the candidate species in the database information table obtained in step (1).

[0012] Preferably, step (6-7) specifically involves: First, based on the read length-strain mapping table obtained after the fourth round of processing in step (6-6), the background rate for data-driven processing is obtained. : in, This represents the number of rows in the read-strain mapping table after the fourth round of processing, i.e., the number of reads remaining after the filtering processes from steps (6-1) to (6-6). The total number of reads in the metagenomic sequencing data input by the user. The number of candidate species obtained in step (6-6); constant 10 6 This is used for numerical scaling on the million-scale, which facilitates calculation and threshold setting; Then, based on the background rate The signal-to-noise ratio for each candidate species obtained in step (6-6) ,in This indicates the number of supporting reads for the candidate species; Finally, all signal-to-noise ratios are preserved. Candidate species, as the final metagenomic species.

[0013] Preferably, step (7) is to obtain the first step using the following formula. Relative abundance of the final metagenomic species : in Indicates the first The number of supporting reads for a final metagenomic species. Indicates the first The representative genome length of a final metagenomic species, This represents all the final metagenomic species obtained in step (6).

[0014] Preferably, the first Representative genome length of a final metagenomic species It is obtained through the following process: First, the total number of all bases in each strain in the database information table obtained in step (1) is summed, and the summation result is taken as the genome length of that strain; then, the number of bases in the database information table is obtained. The final metagenomic species is represented in the database information table along with all corresponding strains; finally, the arithmetic mean of the genome lengths of all strains is taken as the result. The representative genome length of a final metagenomic species.

[0015] In summary, compared with the prior art, the above-described technical solutions conceived by this invention can achieve the following beneficial effects: (1) This invention employs steps (2-6), using the linear genome and pangenome map as complementary reference frames, and adopts a two-stage process: First, the hierarchical staggered Bloom filter (HIBF) is used to perform approximate membership query to screen out candidate species from the database with high recall and low overhead; then, ccDBG is constructed around these candidates to determine the attribution of each read at the strain scale; this design combines the rapid screening capability of linear reference with the high-resolution capture of intrapopulation differences by pangenome, thereby significantly solving the classification ambiguity and confusion problems in traditional processes, and obtaining more accurate metagenomic classification results than existing classification methods while ensuring high recall. Therefore, it can solve the technical problem of difficulty in balancing recall and accuracy in existing methods. (2) Due to the adoption of step (6), this invention proposes a dynamic threshold calculation formula suitable for filtering false positive species in metagenomics. Only when the number of reads supporting a species exceeds the dynamic threshold will it be retained as a result. The calculation process considers the sample community complexity, filtering intensity, and data volume, which can achieve high accuracy and robustness in species classification. It can obtain accurate classification and relative abundance regardless of the change in sample sequencing depth, thus solving the technical problem that classification performance is affected by sequencing depth. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the metagenomic species-level classification method based on the de Bruin diagram of this invention; Figure 2 This presents the species-level analysis results from NGS simulated datasets at different sequencing depths. The x-axis represents the gradient change from ultra-low to high coverage, with UC / LC / MC / HC labels corresponding to ultra-low, low, medium, and high coverage nodes on the axis, respectively. Figure 2 (a) shows the changes in recall, precision, F1 score, and Bray-Curtis distance of sim100 with coverage gradient. Figure 2 (b) shows the changes in recall, precision, F1 score, and Bray-Curtis distance of sim1000 with coverage gradient; Figure 3 These are the species-level analysis results from NGS simulated datasets under different community complexities. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0018] The basic idea of ​​this invention is to provide a metagenomic species-level classification method based on de Bruin diagrams. This method uses linear genomes and pan-genomes as complementary reference frames and employs a two-stage process: First, HIBF is used to perform approximate membership queries to screen candidate species and strains from the database with high recall and low overhead. Then, ccDBG is constructed around these candidates to finely determine the attribution of each read at the strain level. This design combines the rapid screening capability of linear references with the high-resolution capture of intrapopulation differences by pan-genomes, thereby significantly alleviating the problems of species classification ambiguity and confusion in traditional processes.

[0019] like Figure 1As shown, this invention provides a metagenomic species-level classification method based on de Bruin diagrams, comprising the following steps: (1) Obtain the complete genome, and construct the HIBF index of the reference database based on the complete genome and using the hierarchical interleaved Bloom filter (HIBF); Specifically, this step involves: first, obtaining the complete genome; then, obtaining all strains and species within the complete genome; subsequently, establishing a database information table based on all strains and species within the complete genome, including five entries: the name of each strain, the taxonomic number of each strain, the name of each species, the taxonomic number of each species, and the path where the obtained complete genome is located; finally, constructing a HIBF index for the reference database based on the established database information table and the complete genome, using a hierarchical interleaved Bloom filter (HIBF). The advantage of step (1) is that HIBF effectively addresses highly imbalanced reference set structures by splitting excessively large datasets and merging excessively small datasets across multiple levels, while simultaneously improving cache locality through an interleaved storage strategy. This invention employs a window minimum hash scheme when building the index and reuses this scheme during queries, significantly reducing the index size while maintaining matching semantics.

[0020] (2) Obtain the metagenomic sequencing data input by the user, and based on the HIBF index of the reference database obtained in step (1), obtain the strain corresponding to each read length (in English, read, representing sequencing sequence) in the metagenomic sequencing data, and perform false positive filtering on all strains corresponding to the obtained read lengths, so as to obtain multiple initial screening species, the strain corresponding to each initial screening species, and the read length corresponding to the strain.

[0021] This step includes the following sub-steps: (2-1) Obtain the metagenomic sequencing data input by the user, and use the k-mer method to decompose each read in the metagenomic sequencing data to obtain the k-mer set corresponding to the read; (2-2) Input the set of k-mers corresponding to each read length obtained in step (2-1) into the HIBF index of the reference database constructed in step (1) to perform an approximate member query to obtain the mapping relationship and k-mer matching results between the read length and each strain in the HIBF index of the reference database. (2-3) Obtain the length L_i of the i-th read length obtained in step (2-1), the length k of the k-mer set corresponding to the i-th read length, and obtain the total number of k-mers in the i-th read length. Obtain the score of the j-th strain in the HIBF index of the reference database on the i-th read. Where i∈[1, total number of reads in metagenomic sequencing data], j∈[1, total number of strains in the HIBF index of the reference database], This represents the number of k-mer matches between the i-th read and the j-th strain; (2-4) For the i-th read, all strains in the HIBF index of the reference database that have a mapping relationship with it and have a score greater than 0.75 are retained as high-scoring strains. The maximum score among all high-scoring strains is obtained, and high-scoring strains with scores less than or equal to the maximum score * 0.9 are deleted. The result is used as the co-leading strain. All co-leading strains are used as the initial screening strains, and the species corresponding to all co-leading strains are used as the initial screening species. (3) Use the Expectation-Maximization (EM) algorithm to obtain the relative abundance between each primary screening species and its corresponding primary screening strain obtained in step (2), and perform threshold filtering on all the obtained relative abundance to obtain multiple candidate species and multiple candidate strains. This step includes the following sub-steps: (3-1) Based on the primary screening species, their corresponding primary screening strains and their read lengths obtained in step (2), and using the EM algorithm to solve and normalize them, the relative abundance between each primary screening species and its corresponding primary screening strains is obtained. (3-2) Based on the relative abundance of each primary screening species and its corresponding primary screening strain obtained in step (3-1), all primary screening species are filtered, i.e., those with a relative abundance greater than 10 are retained. -7 The initial screening species were used as candidate species; when the total number of candidate species exceeded 1000, the three strains with the highest relative abundance of all candidate species were retained as candidate strains; when the number of candidate species was less than or equal to 1000, the strains with a relative abundance greater than 10 were retained. -7 All initially screened strains were used as candidate strains; The advantage of step (3-2) is that it can obtain a small-scale and reliable strain set, thereby effectively reducing the computational cost of ccDBG construction in the subsequent step (4) and effectively reducing the background noise in the subsequent step (6).

[0022] (4) Based on the multiple candidate strains obtained in step (3) and the complete genomes of multiple strains obtained in step (1), use the GGCAT tool (specifically the building module) to construct colored and compacted de Bruijn graphs (ccDBG for short). The ccDBG graph records multiple k-mers and the strains to which each k-mer belongs, where the value of k is 19 or higher, preferably 31. (5) Obtain the metagenomic sequencing data input by the user, and use the GGCAT tool (specifically the query module) to split each read in the metagenomic sequencing data into k-mers to obtain the set of k-mers corresponding to the read. In the ccDBG graph constructed in step (4), obtain multiple k-mers that match the set of k-mers corresponding to the read, as well as the strains to which each k-mer belongs. All strains corresponding to the read in the ccDBG graph are the strains associated with the read. Finally, organize the above information into a read-strain mapping table. Specifically, establish a read-strain mapping table based on the number of all reads in the metagenomic sequencing data and the number of strains associated with them. Each row of the read-strain mapping table represents a different read in the metagenomic sequencing data, and the columns represent the strains associated with each read in the ccDBG. For example, suppose a read of length ACTGCTGCTG, k=3, then the set of k-mers corresponding to this read is {ACT, CTG, TGC, GCT, CTG, TGC, GCT, CTG}. In the ccDBG graph, all k-mers matching this set are ACT, CTG, and TGC. Specifically, the strains corresponding to k-mer ACT in the ccDBG graph are A and B, and the strains corresponding to k-mer CTG in the ccDBG graph are B and C. If TGC belongs to strain C in the ccDBG graph, then the strains associated with this read length are A, B, and C. In the established read length-strain mapping table, the value of read1 in the "strain A" column is 1 because there is only one matching k-mer {ACT}; the value of read1 in the "strain B" column is 2 because there are two matching k-mers {ACT,CTG}; the value of read1 in the "strain C" column is 2 because there are two matching k-mers {CTG,TGC}.

[0023] In general, in the read-strain mapping table established in this step, for each row (i.e. each read), the number of strains associated with that read will be recorded in its corresponding columns.

[0024] (6) Obtain the final metagenomic species based on the read length-strain mapping table obtained in step (5).

[0025] This step includes the following sub-steps: (6-1) Obtain the maximum value of each row from the read length-strain mapping table obtained in step (5), and set the other values ​​in the row to 0, keeping only the maximum value. If there are multiple maximum values ​​in a row, keep all the maximum values, thus obtaining the read length-strain mapping table after the first round of processing. (6-2) Divide all values ​​in the read-strain mapping table obtained in step (6-1) after the first round of processing by... To obtain the updated read-strain mapping table; where This indicates the length of each read (in bases per bp).

[0026] (6-3) Take the maximum value of each row in the updated read-strain mapping table obtained in step (6-2) as the score corresponding to that row, and calculate the 30th percentile of the scores corresponding to all rows. Then the 30th percentile Cut off to interval Within this range, discard entries in the updated read-strain mapping table that score less than or equal to the 30th percentile. All rows are processed to obtain the read-strain mapping table after the second round of processing; The purpose of calculating the 30th percentile of the scores for all rows in step (6-3) is that each row represents a read. Some low-quality reads, while associated with certain strains, have a small number of matching strains, meaning the read has low similarity to the strain. Filtering each row using the 30th percentile can eliminate results from low-quality reads. Furthermore, this 30th percentile... Cut off to interval This is to avoid filtering too strongly or too weakly.

[0027] (6-4) Delete all rows in the read length-strain mapping table obtained in the second round of processing obtained in step (6-3) that do not meet the first preset condition, wherein the first preset condition is that all strains in the row correspond to the same species in the database information table obtained in step (1), so as to obtain the read length-strain mapping table after the third round of processing. (6-5) Determine whether the metagenomic sequencing data obtained in step (5) is paired-end sequencing data or single-end sequencing data. If it is paired-end sequencing data, first delete all rows in the read length-strain mapping table after the third round of processing obtained in step (6-4) that do not meet the second preset condition. The second preset condition is that all strains in the row correspond to the same species in the database information table obtained in step (1). Then delete the columns and rows with all zeros in the read length-strain mapping table after the third round of processing to obtain the read length-strain mapping table after the fourth round of processing. Then proceed to step (6-6). If it is single-end sequencing data, delete the columns and rows with all zeros in the read length-strain mapping table after the third round of processing to obtain the read length-strain mapping table after the fourth round of processing. Then proceed to step (6-6). In sub-step (6-4) above, only the read-strain mapping table after the second round of processing is retained. The purpose of all values ​​being distributed in rows of the same species is that each row represents a read. If a read is associated with strains from multiple species, this situation may arise from poor read quality, poor genome quality used to construct ccDBG, or very similar genomes between strains of two species. In any case, it means that it is difficult to determine which species the read belongs to, so the read is considered interfering and is therefore deleted. Similarly, in sub-step (6-5) above, paired-end sequencing data shows reads in pairs, which appear in two rows in the table. A pair of reads must come from the same species. However, if a pair of reads is associated with strains of different species in the read-strain mapping table, it means that one of the reads in the pair is problematic. However, it is impossible to determine which one is problematic in the table, so the pair of reads is considered interfering and is therefore deleted as a whole.

[0028] (6-6) Set all the remaining strains in the read length-strain mapping table after the fourth round of processing as candidate species corresponding to all species in the database information table obtained in step (1), and obtain the number of supported read lengths for each candidate species; Specifically, the number of supported reads for a candidate species is equal to the number of rows with values ​​greater than zero in the corresponding column of the read-strain mapping table after the fourth round of processing, corresponding to the strain of the candidate species in the database information table obtained in step (1).

[0029] (6-7) Based on the number of supporting reads for each candidate species obtained in step (6-6), the read-strain mapping table obtained in the fourth round of processing in step (6-5) is further filtered to obtain the metagenomic species classification results of the metagenomic sequencing data. Specifically, this step involves first obtaining the background rate for data-driven processing based on the read length-strain mapping table obtained in step (6-6) after the fourth round of processing. : in, This represents the number of rows in the read-strain mapping table after the fourth round of processing, i.e., the number of reads remaining after the filtering processes from steps (6-1) to (6-6). The total number of reads in the metagenomic sequencing data input by the user. The number of candidate species obtained in step (6-6); constant 10 6 This is used for numerical scaling on the million-scale, which facilitates calculation and threshold setting; Then, based on the background rate The signal-to-noise ratio (SNR) for each candidate species obtained in step (6-6) ,in This indicates the number of supporting reads for the candidate species; Finally, all signal-to-noise ratios are preserved. Candidate species, as the final metagenomic species.

[0030] The advantage of step (6-6) is that it innovatively and reasonably combines the amount of metagenomic sample data, sample community complexity, and the strength of the initial filtering, dynamically calculates the filtering threshold for each sample, and avoids excessive penalties on rare real species while effectively eliminating noise-driven false positive species, thus ensuring the recall rate. This allows for accurate final metagenomic species to be obtained at different sequencing depths and community complexities.

[0031] (7) Obtain the first result obtained in step (6) Relative abundance of the final metagenomic species As a classification result, among ∈[1, the total number of all final metagenomic species obtained in step (7)].

[0032] Specifically, this step uses the following formula to obtain the first... Relative abundance of the final metagenomic species : in Indicates the first The number of supporting reads for a final metagenomic species. Indicates the first The representative genome length of a final metagenomic species, This represents all the final metagenomic species obtained in step (6).

[0033] No. Representative genome length of a final metagenomic species It is obtained through the following process: First, the total number of all bases in each strain in the database information table obtained in step (1) is summed, and the summation result is taken as the genome length of that strain; then, the number of bases in the database information table is obtained. The final metagenomic species is represented in the database information table along with all corresponding strains; finally, the arithmetic mean of the genome lengths of all strains is taken as the result. The representative genome length of a final metagenomic species.

[0034] The advantage of this step (7) is that it can reduce the bias caused by differences in genome size, make the sum of the composition of the obtained sample equal to 1 across species, and allow direct comparison between different datasets and sequencing platforms.

[0035] Experimental results The following will describe the experiments performed according to the specific embodiments of the present invention, and a comparative analysis will be conducted with existing methods, mainly reflected in the following four indicators: (1) Precision and Recall Precision and recall are defined as follows: at a certain taxonomic level (such as species level or strain level), precision assesses the proportion of taxa correctly predicted by the method out of all taxa predicted by the method, while recall assesses the proportion of taxa correctly predicted by the method out of the actual taxa present in the metagenomic sequencing data.

[0036] The calculation formula is as follows: In this context, TP represents a taxonomic unit correctly predicted by the method, FP represents a taxonomic unit that is not present in the metagenomic sequencing data but is predicted by the method, and FN represents a taxonomic unit that exists in the metagenomic sequencing data but is not predicted.

[0037] (2) F1 score The F1-score is the harmonic mean of precision and recall, and it is a core metric for comprehensively evaluating the accuracy and completeness of prediction results. When a model exhibits an imbalance between precision and recall (e.g., high precision but low recall, or vice versa), using the F1-score is more valuable than using either metric alone. Its value ranges from [0,1], with values ​​closer to 1 indicating better overall model performance.

[0038] (3) Bray–Curtis distance (BC_distance) The calculation formula is as follows: in, The representation is the union of the taxonomic units predicted by the method and the actual taxonomic units present in the metagenomic sequencing data. For A certain classification unit in, using and Representing classification units The actual and predicted relative abundance.

[0039] Example 1 The experimental data were prepared as follows: First, sequencing depth was divided into four levels: 0.05× average coverage was defined as ultra-low coverage (UC), 0.5× as low coverage (LC), 2× as medium coverage (MC), and 5× as high coverage (HC). Second, two datasets were generated according to classification complexity: sim100 (45 species / 100 strains) and sim1000 (373 / 1000) NGS datasets. To further characterize the scalability of the method, continuous depth gradient experiments were conducted: the data size of sim100 ranged from 0.01Gb to 1Gb (corresponding to average coverage of 0.05×–5×), and the data size of sim1000 ranged from 0.1Gb to 10Gb (0.05×–5×). To reflect the difference between real samples and the reference database, a 0.5% genomic mutation rate was introduced into the above datasets. The CAMISIM tool was used to simulate second-generation short-read NGS data.

[0040] The specific implementation steps of this embodiment are performed in accordance with the specific implementation method of the present invention, so as to obtain the species composition corresponding to the metagenomic sequencing data of the present invention and the relative abundance of each species in the species composition.

[0041] For each of the other methods, the experimental steps of this method are as follows: (1) Obtain the reference dataset constructed in this example, and construct the reference dataset index to obtain the reference database of this method.

[0042] Note that MetaPhlAn4 uses a reference database that it constructs itself as the reference database for species-level experiments using this method.

[0043] (2) Obtain the metagenomic sequencing data input by the user, use a reference database, and perform species-level classification according to the specific implementation method of the method to obtain the species composition corresponding to the metagenomic sequencing data and the relative abundance of each species in the species composition.

[0044] Table 1-1: Species-level experimental results of the sim100-UCNGS dataset Table 1-2: Species-level experimental results of the sim100-LCNGS dataset Table 1-3: Species-level experimental results of the sim100-MCNGS dataset Table 1-4: Species-level experimental results of the sim100-HCNGS dataset Table 1-5: Species-level experimental results of the sim1000-UCNGS dataset Table 1-6: Species-level experimental results of the sim1000-LCNGS dataset Table 1-7: Species-level experimental results of the sim1000-MCNGS dataset Table 1-8: Species-level experimental results of the sim1000-MCNGS dataset In this experiment, the best results for each metric are indicated in bold. Themis, the invention described above, performed best in all F1-scores. Regarding recall, while not the best, it was sometimes slightly worse than four tools that prioritized high recall: Ganon, Bracken, Kraken2, and Centrifuger. These four tools achieved slightly higher recall, but at the cost of significant precision. In terms of precision, it consistently performed best on low-complexity datasets like sim100. On high-complexity datasets like sim1000, while not always the best, it was sometimes slightly worse than two tools that prioritized high precision: MetaPhlAn4 and KMCP. These two tools achieved slightly higher precision, but at the cost of significant recall. BC distance performs best in ultra-low and low coverage datasets, and slightly worse than KMCP in medium and high coverage datasets. KMCP is more suitable for datasets with high coverage; to ensure accuracy, it often removes rare species with relatively low abundance. This results in many rare species not being detected, but high-abundance species being accurately identified, leading to relatively accurate relative abundance predictions. Therefore, BC distance may be inflated. Furthermore, it can be seen that only Themis maintains consistently high performance across ultra-low and high coverage datasets, while other tools show significant variations with changes in sample data coverage. In summary, Themis is the only tool that can consistently achieve high performance in species-level classification across datasets of varying complexity, regardless of coverage.

[0045] To better illustrate the results of this embodiment, a simple R script was used to visualize the results and compare them with the distribution of real species and strains in the dataset, such as... Figure 2 As shown, where Figure 2 'a' represents the result of sim100. Figure 2 b represents the result of sim1000.

[0046] Example 2 The experimental data were prepared as follows: five datasets were generated according to classification complexity: sim100 (45 species / 100 strains), sim500 (195 / 500), sim1000 (373 / 1000), sim2000 (805 / 2000), and sim3000 (1243 / 3000). To reflect the differences between real samples and the reference database, a 0.5% genomic mutation rate was introduced into the above datasets. The CAMISIM tool was used to simulate second-generation short-read NGS data with an ultra-low coverage (UC) of 0.05× for datasets of different complexity.

[0047] The specific implementation steps of this embodiment are performed in accordance with the specific implementation method of the present invention, so as to obtain the species composition corresponding to the metagenomic sequencing data of the present invention and the relative abundance of each species in the species composition.

[0048] The specific implementation steps of other methods are performed according to Example 1 (not repeated here) to obtain the species composition and relative abundance of each species in the metagenomic sequencing data of other methods.

[0049] Table 2-1: Species-level experimental results of the sim100-UCNGS dataset Table 2-2: Species-level experimental results of the sim500-UCNGS dataset Table 2-3: Species-level experimental results of the sim1000-UCNGS dataset Table 2-4: Species-level experimental results of the sim2000-UCNGS dataset Table 2-5: Species-level experimental results of the sim3000-UCNGS dataset In this experiment, the optimal results for each metric are indicated in bold. At extremely low coverage, the Themis of this invention performs best in both F1-score and BC distance. Furthermore, its performance remains robust regardless of changes in cluster complexity. To better illustrate the results of this embodiment, a simple R script is used to visualize the experimental results and compare them with the actual species and strain distributions in the dataset, such as... Figure 3 As shown.

[0050] These benchmark results show that Themis outperforms existing state-of-the-art methods by a wide margin, primarily by achieving higher F1-scores and BC distances at both ultra-low and high sequencing depths, while maintaining comparable or better performance in other aspects.

[0051] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A metagenomic species-level classification method based on de Bruin diagrams, characterized in that, Includes the following steps: (1) Obtain the complete genome, and construct the HIBF index of the reference database based on the complete genome and using the hierarchical staggered Bloom filter HIBF; (2) Obtain the metagenomic sequencing data input by the user, and based on the HIBF index of the reference database obtained in step (1), obtain the strain corresponding to each read length in the metagenomic sequencing data, and perform false positive filtering on all the strains corresponding to the obtained read lengths, so as to obtain multiple initial screening species, the strains corresponding to each initial screening species, and the read lengths corresponding to the strains. (3) Use the expectation-maximization EM algorithm to obtain the relative abundance between each primary screening species and its corresponding primary screening strain obtained in step (2), and perform threshold filtering on all the obtained relative abundance to obtain multiple candidate species and multiple candidate strains. (4) Based on the multiple candidate strains obtained in step (3) and the complete genomes of multiple strains obtained in step (1), a color and compressed de Bruin diagram ccDBG is constructed using the GGCAT tool. The ccDBG diagram records multiple k-mers and the strains to which each k-mer belongs, where the value of k is 19 or higher. (5) Obtain the metagenomic sequencing data input by the user, use the GGCAT tool to split each read in the metagenomic sequencing data into k-mers to obtain the set of k-mers corresponding to the read, and obtain multiple k-mers that match the set of k-mers corresponding to the read in the ccDBG graph constructed in step (4), as well as the strains to which each k-mer belongs. All strains corresponding to the read in the ccDBG graph are the strains associated with the read. Finally, organize the above information into a read-strain mapping table. Specifically, establish a read-strain mapping table based on the number of all reads in the metagenomic sequencing data and the number of strains associated with them. Each row of the read-strain mapping table is a different read in the metagenomic sequencing data, and the columns are the strains associated with each read in the ccDBG. (6) Obtain the final metagenomic species based on the read length-strain mapping table obtained in step (5); (7) Obtain the first result obtained in step (6) Relative abundance of the final metagenomic species As a classification result, among ∈[1, the total number of all final metagenomic species obtained in step (7)].

2. The metagenomic species-level classification method based on de Bruin diagrams according to claim 1, characterized in that, Step (1) specifically involves: first, obtaining the complete genome; then, obtaining all strains and species in the complete genome; subsequently, establishing a database information table based on all strains and species in the complete genome, including five entries: the name of each strain, the classification number of each strain, the name of each species, the classification number of each species, and the path where the obtained complete genome is located; finally, constructing a HIBF index for the reference database based on the established database information table and the complete genome, and using a hierarchical staggered Bloom filter (HIBF).

3. The metagenomic species-level classification method based on de Bruin diagrams according to claim 1 or 2, characterized in that, Step (2) includes the following sub-steps: (2-1) Obtain the metagenomic sequencing data input by the user, and use the k-mer method to decompose each read in the metagenomic sequencing data to obtain the k-mer set corresponding to the read; (2-2) Input the set of k-mers corresponding to each read length obtained in step (2-1) into the HIBF index of the reference database constructed in step (1) to perform an approximate member query to obtain the mapping relationship and k-mer matching results between the read length and each strain in the HIBF index of the reference database. (2-3) Obtain the length L_i of the i-th read length obtained in step (2-1), the length k of the k-mer set corresponding to the i-th read length, and obtain the total number of k-mers in the i-th read length. Obtain the score of the j-th strain in the HIBF index of the reference database on the i-th read. Where i∈[1, total number of reads in metagenomic sequencing data], j∈[1, total number of strains in the HIBF index of the reference database], This represents the number of k-mer matches between the i-th read and the j-th strain; (2-4) For the i-th read, retain all strains in the HIBF index of the reference database that have a mapping relationship with it and whose scores are greater than 0.75 as high-scoring strains. Obtain the maximum score among all high-scoring strains, and multiply the scores of all high-scoring strains that are less than or equal to the maximum score by 1. High-scoring strains (0.9) are deleted, and the resulting strains are used as co-leading strains. All co-leading strains are used as primary screening strains, and the species corresponding to all co-leading strains are used as primary screening species.

4. The metagenomic species-level classification method based on de Bruin diagrams according to any one of claims 1 to 3, characterized in that, Step (3) includes the following sub-steps: (3-1) Based on the primary screening species, their corresponding primary screening strains and their read lengths obtained in step (2), and using the EM algorithm to solve and normalize them, the relative abundance between each primary screening species and its corresponding primary screening strains is obtained. (3-2) Based on the relative abundance of each primary screening species and its corresponding primary screening strain obtained in step (3-1), all primary screening species are filtered, i.e., those with a relative abundance greater than 10 are retained. -7 The initial screening species were used as candidate species; when the total number of candidate species exceeded 1000, the three strains with the highest relative abundance of all candidate species were retained as candidate strains; when the number of candidate species was less than or equal to 1000, the strains with a relative abundance greater than 10 were retained. -7 All the initial screening strains were used as candidate strains.

5. The metagenomic species-level classification method based on de Bruin diagrams according to claim 4, characterized in that, Step (6) includes the following sub-steps: (6-1) Obtain the maximum value of each row from the read length-strain mapping table obtained in step (5), and set the other values ​​in the row to 0, keeping only the maximum value. If there are multiple maximum values ​​in a row, keep all the maximum values, thus obtaining the read length-strain mapping table after the first round of processing. (6-2) Divide all values ​​in the read-strain mapping table obtained in step (6-1) after the first round of processing by... To obtain the updated read-strain mapping table; where This indicates the length of each read, expressed in bases (bp). (6-3) Take the maximum value of each row in the updated read-strain mapping table obtained in step (6-2) as the score corresponding to that row, and calculate the 30th percentile of the scores corresponding to all rows. Then the 30th percentile Cut off to interval Within this range, discard entries in the updated read-strain mapping table that score less than or equal to the 30th percentile. All rows are processed to obtain the read-strain mapping table after the second round of processing; (6-4) Delete all rows in the read length-strain mapping table obtained in the second round of processing obtained in step (6-3) that do not meet the first preset condition, wherein the first preset condition is that all strains in the row correspond to the same species in the database information table obtained in step (1), so as to obtain the read length-strain mapping table after the third round of processing. (6-5) Determine whether the metagenomic sequencing data obtained in step (5) is paired-end sequencing data or single-paired-end sequencing data. If it is paired-end sequencing data, first delete all rows in the read length-strain mapping table after the third round of processing obtained in step (6-4) that do not meet the second preset condition. The second preset condition is that all strains in the row correspond to the same species in the database information table obtained in step (1). Then delete the columns and rows with all zeros in the read length-strain mapping table after the third round of processing to obtain the read length-strain mapping table after the fourth round of processing. Then proceed to step (6-6). If it is single-paired-end sequencing data, delete the columns and rows with all zeros in the read length-strain mapping table after the third round of processing to obtain the read length-strain mapping table after the fourth round of processing. Then proceed to step (6-6). (6-6) Set all the remaining strains in the read length-strain mapping table after the fourth round of processing as candidate species corresponding to all species in the database information table obtained in step (1), and obtain the number of supported read lengths for each candidate species; (6-7) Based on the number of supporting reads for each candidate species obtained in step (6-6), the read-strain mapping table obtained in the fourth round of processing in step (6-5) is further filtered to obtain the metagenomic species classification results of the metagenomic sequencing data.

6. The metagenomic species-level classification method based on de Bruin diagrams according to claim 5, characterized in that, The number of supported reads for a candidate species is equal to the number of rows with values ​​greater than zero in the corresponding column of the read-strain mapping table after the fourth round of processing, corresponding to the strain of the candidate species in the database information table obtained in step (1).

7. The metagenomic species-level classification method based on de Bruin diagrams according to claim 6, characterized in that, Steps (6-7) are as follows: First, based on the read length-strain mapping table obtained after the fourth round of processing in step (6-6), the background rate for data-driven processing is obtained. : ； in, This represents the number of rows in the read-strain mapping table after the fourth round of processing, i.e., the number of reads remaining after the filtering processes from steps (6-1) to (6-6). The total number of reads in the metagenomic sequencing data input by the user. The number of candidate species obtained in step (6-6); constant 10 6 This is used for numerical scaling on the million-scale, which facilitates calculation and threshold setting; Then, based on the background rate The signal-to-noise ratio for each candidate species obtained in step (6-6) ,in This indicates the number of supporting reads for the candidate species; Finally, all signal-to-noise ratios are preserved. Candidate species, as the final metagenomic species.

8. The metagenomic species-level classification method based on de Bruin diagrams according to claim 7, characterized in that, Step (7) is to obtain the first step using the following formula. Relative abundance of the final metagenomic species : ； in Indicates the first The number of supporting reads for a final metagenomic species. Indicates the first The representative genome length of a final metagenomic species, This represents all the final metagenomic species obtained in step (6).

9. The metagenomic species-level classification method based on de Bruin diagrams according to claim 8, characterized in that, No. Representative genome length of a final metagenomic species It is obtained through the following process: First, the total number of all bases in each strain in the database information table obtained in step (1) is summed, and the summation result is taken as the genome length of that strain; then, the number of bases in the database information table is obtained. The final metagenomic species is represented in the database information table along with all corresponding strains; finally, the arithmetic mean of the genome lengths of all strains is taken as the result. The representative genome length of a final metagenomic species.

Images