A method for strain-level classification of metagenomic data based on pan-genome graph

Abstract

The application discloses a method for strain-level classification of metagenomic data based on a pan-genome graph, which first captures the genomes of multiple strains of each species as a reference through the technology of the pan-genome graph, instead of using only the representative genome of each species as a reference, expands the strain diversity, and enhances the alignment accuracy of metagenomic sequencing reads and the estimation of composition at the strain level; then, by using the corresponding short read alignment tool Giraffe and the long read alignment tool GraphAligner in the alignment of NGS and TGS metagenomic data with the pan-genome graph, two types of data can be efficiently processed at the same time; finally, by further optimizing the path abundance of the species-level classification results based on the pan-genome graph, multi-species strain-level classification can be realized. The application can solve the technical problem that the strain-level classification capability is limited due to the strategy of a single reference genome in the prior art.

Description

Technical Field

[0001] This invention belongs to the field of metagenomic classification and analysis technology, and more specifically, relates to a method for classifying metagenomic data at the strain level based on pangenome maps. Background Technology

[0002] Microorganisms are widely distributed across various environments on Earth, such as oceans, soil, and the gastrointestinal tract, playing an indispensable role in maintaining ecological balance and host health. Microbial communities typically consist of a variety of microorganisms, exhibiting uneven species abundance and high diversity. Due to rapid mutation and horizontal gene transfer within microorganisms, different strains may emerge within the same species, further increasing the biodiversity and complexity of microbial communities. Previous studies have shown that different strains may exhibit different phenotypes or perform different biological functions in the environment. For example, the *Escherichia coli* strain *E. coli* O104 is pathogenic, while common *E. coli* strains are harmless. Furthermore, the distribution of strains varies across different parts of the human body. Therefore, strain-level classification of microbial communities is crucial for understanding their dynamics.

[0003] To reveal the structure and function of microbial communities, it is necessary to estimate their composition and abundance at both the species and strain levels. Metagenomic sequencing data preserves most of the genetic information of microbial communities, enabling strain-level classification; metagenomic sequencing data includes two types: second-generation (NGS) short reads and third-generation (TGS) long reads.

[0004] Classification methods at the species level and above used in NGS can be divided into three categories: the first category is marker-based methods, such as MetaPhlAn4; the second category is DNA-to-protein methods, such as Kaiju; and the third category is DNA-to-DNA methods, such as Kraken2 and CLARK. Although the aforementioned methods are based on NGS, methods such as Kraken2 and CLARK can also be used in TGS. Single-species strain-level classification methods specifically for NGS include, for example, DNA-to-DNA based StrainScan. Classification methods at the species level and above specifically for TGS can also be divided into three categories: marker-based Melon methods, DNA-to-protein based MEGAN-LR methods, and DNA-to-DNA based MetaMaps methods. MetaMaps is also used for strain-level classification in TGS. It is worth noting that Centrifuge is a method that is applicable to both NGS and TGS and can perform classification at both the species and strain levels.

[0005] However, the aforementioned existing methods generally suffer from some significant technical problems:

[0006] First, the methods mentioned above usually select a complete genome as a representative of each species. Since this single reference genome strategy cannot capture the genomic diversity of different strains within a species, it limits the ability of these methods to perform strain-level classification.

[0007] Secondly, these methods, such as MetaPhlAn4, Kaiju, Kraken2, CLARK, StrainScan, Melon, MEGAN-LR, and MetaMaps, cannot efficiently process both NGS and TGS metagenomic data simultaneously. When processing TGS data, methods developed based on NGS data often experience performance degradation due to incomplete adaptation of algorithms and parameters. Furthermore, k-mer-based methods cannot fully utilize the advantages of TGS data.

[0008] Third, these methods, such as MetaPhlAn4, Kaiju, Kraken2, CLARK, Melon, and MEGAN-LR, can only perform classification at the species level, and are difficult to achieve at the strain level or even the multi-species strain level. Furthermore, the StranScan method can only perform classification at the single-species strain level, which also has limitations. Summary of the Invention

[0009] To address the aforementioned shortcomings or improvement needs of existing technologies, this invention provides a method for classifying metagenomic data at the strain level based on a pangenome graph (PanTax). Its purpose is to solve the technical problems of limited strain-level classification capabilities due to the single reference genome strategy of existing methods; the difficulty in efficiently processing both NGS and TGS data types simultaneously due to their different characteristics; and the difficulty in performing strain-level classification, or even multi-species strain-level classification, due to insufficient development and high strain similarity.

[0010] To achieve the above objectives, according to one aspect of the present invention, a method for performing strain-level analysis of metagenomic data based on a pan-genome map is provided, comprising the following steps:

[0011] (1) Obtain multiple species and their corresponding complete genomes from public databases, and construct a reference database based on all the obtained complete genomes.

[0012] (2) For each species in the reference database obtained in step (1), construct the pangenome map corresponding to the species using all the reference genomes corresponding to the species in the reference database, and merge the pangenome maps corresponding to all species to obtain the reference pangenome map.

[0013] (3) Obtain the metagenomic data input by the user, and compare the reference pangenome map obtained in step (2) with the metagenomic data to obtain multiple alignment results for each sequencing read in the metagenomic data.

[0014] (4) Based on all alignment results of each sequencing read in the metagenomic data obtained in step (3), perform species-level classification and binning to obtain the final species composition corresponding to the metagenomic data.

[0015] (5) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), perform species-level taxonomic spectrum analysis on the species to obtain the relative abundance of the species.

[0016] (6) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), the species is subjected to path abundance optimization (PAO) processing according to the pangenome map corresponding to the species obtained in step (2) to obtain the final strain composition of the metagenomic data and the absolute abundance of each strain in the final strain composition.

[0017] (7) Based on the final species composition corresponding to the metagenomic data obtained in step (4), perform abundance adjustment and normalization calculation on the final strain composition of the metagenomic data obtained in step (6) and the absolute abundance of each strain in the final strain composition to obtain the final strain composition of the metagenomic data and the relative abundance of each strain in the final strain composition.

[0018] Preferably, step (1) includes the following sub-steps:

[0019] (1-1) Download and obtain multiple species and multiple complete genomes corresponding to each species from public databases;

[0020] (1-2) Remove plasmids from all the complete genomes obtained in step (1-1) to obtain multiple complete genomes after plasmid removal;

[0021] (1-3) For each species obtained in step (1-1), the average nucleotide similarity (ANI) between any pair of complete genomes corresponding to that species is calculated using the fastANI tool, and is used as the distance score between that pair of complete genomes corresponding to that species.

[0022] (1-4) For each species obtained in step (1-1), a heuristic clustering algorithm based on connected graphs is used to cluster the distance scores between all pairs of complete genomes corresponding to the species obtained in step (1-3) to obtain multiple complete genome sets. The complete genome set containing the most complete genomes is selected from all complete genome sets, and a predetermined number of complete genomes are selected from the selected complete genome set as multiple reference genomes corresponding to the species. All reference genomes corresponding to all species constitute a reference database.

[0023] Preferably, in steps (1-4), if the clustering result is less than a preset first clustering threshold, it indicates that the pair of complete genomes corresponding to the clustering result are strains of different species; if the clustering result is greater than or equal to the preset first clustering threshold and less than a preset second clustering threshold, it indicates that the pair of complete genomes corresponding to the clustering result are different strains of the same species; if the clustering result is greater than the preset second clustering threshold, it indicates that the pair of complete genomes corresponding to the clustering result are the same strain.

[0024] The clustering process in steps (1-4) involves obtaining clustering results that are greater than or equal to the preset first clustering threshold and less than the preset second clustering threshold as multiple complete genome sets obtained from the clustering process.

[0025] Preferably, step (2) includes the following sub-steps:

[0026] (2-1) For each species in the reference database obtained in step (1), if the species has more than one reference genome in the reference database, the PGGB tool is used to construct the pan-genome map corresponding to the species for all reference genomes; if the species has only one reference genome in the reference database, the reference genome is divided into fragments of 1024 bp in length, and each fragment is used as a node to construct the pan-genome map, thereby obtaining the pan-genome map corresponding to the species.

[0027] (2-2) Merge the pangenome maps of all species obtained in step (2-1) to obtain a reference pangenome map.

[0028] In step (3), for second-generation short-read metagenomic data, the Giraffe tool is used to compare the two; specifically, the vg tool is first used to build an index for the reference pan-genome map to accelerate the comparison process, and then the Giraffe tool is used to compare the short-read metagenomic data with the reference pan-genome map to obtain multiple comparison results for each sequencing read in the short-read metagenomic data;

[0029] In step (3), for third-generation long-read metagenomic data, the GraphAligner tool is used to compare the two to obtain multiple alignment results for each sequencing read in the long-read metagenomic data.

[0030] Step (4) includes the following sub-steps:

[0031] (4-1) For each sequencing read in the metagenomic data obtained in step (3), the optimal alignment result of the sequencing read is obtained based on all alignment results of the sequencing read.

[0032] (4-2) For each sequencing read in the metagenomic data obtained in step (3), the species corresponding to the sequencing read is obtained according to the best alignment result of the sequencing read obtained in step (4-1). The species corresponding to all sequencing reads in the metagenomic data constitute the first species composition corresponding to the metagenomic data (which has potential false positives).

[0033] (4-3) Perform false positive species filtering on the first species composition corresponding to the metagenomic data obtained in step (4-2) to obtain the second species composition as the final species composition corresponding to the metagenomic data.

[0034] Preferably, in step (4-1), if a sequencing read has multiple alignment results, only the best alignment result is retained to ensure that each sequencing read corresponds to one species. Specifically, in the alignment results of short reads, each sequencing read has only one alignment result, so this alignment result is the best alignment result. In the alignment results of long reads, the best alignment result is determined from all alignment results by the alignment quality MAPQ value, the number of matched bases, identity, and sequencing read length.

[0035] Step (4-3) specifically involves obtaining the MAPQ values ​​of multiple sequencing reads corresponding to each species in the first species composition from the best alignment results of multiple sequencing reads corresponding to that species. If at least one sequencing read has a MAPQ value that reaches the highest score of 60, and at least one-tenth of the sequencing reads have a MAPQ value greater than 2, then the species is identified as true; otherwise, the species is identified as a false positive and then filtered.

[0036] Preferably, step (5) includes the following sub-steps:

[0037] (5-1) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), the absolute abundance of the species is obtained based on the best alignment result of each sequencing read corresponding to the species obtained in step (4-1).

[0038] This step uses the following formula to calculate the absolute abundance c of the i-th species in the final species composition. i :

[0039]

[0040] Where i∈[1, the total number of species in the final species composition T], n i It is the total number of sequencing reads corresponding to the i-th species, j∈[1, n] i ], n ij It is the j-th sequencing read corresponding to the i-th species, |n ij | is the number of bases in the j-th sequencing read corresponding to the i-th species, m i It is the total number of reference genomes corresponding to the i-th species, k∈[1, m] i ], g ik It is the k-th reference genome corresponding to the i-th species, |g ik | is the number of bases in the k-th reference genome corresponding to the i-th species.

[0041] (5-2) Normalize the absolute abundance of each species in the final species composition obtained in step (5-1) to obtain the relative abundance of that species.

[0042] This step specifically uses the following formula to calculate the relative abundance a of the i-th species in the final species composition. i :

[0043]

[0044] Preferably, step (6) includes the following sub-steps:

[0045] (6-1) Set the counter i = 1;

[0046] (6-2) Determine whether i is greater than the total number of species T corresponding to the final species composition of the metagenomic data. If so, the process ends; otherwise, proceed to step (6-3).

[0047] (6-3) For the i-th species in the final species composition, calculate the observed abundance of each node in the pangenome graph corresponding to the species based on the best alignment results of each sequencing read corresponding to the species obtained in step (4-1).

[0048] (6-4) For the i-th species in the final species composition, the false positive strains in the species are filtered according to the pan-genome map of the species obtained in step (2-1) and the optimal alignment results of each sequencing read of the species obtained in step (4-1) to obtain the strain composition of the species after the first false positive filtering.

[0049] (6-5) For the i-th species in the final species composition, based on the observed abundance of each node in the pangenome graph corresponding to the species obtained in step (6-3), perform the first PAO solution on the strain composition after the first false positive filtering obtained in step (6-4) for the species, and retain the strains with absolute abundance greater than zero according to the solution results, so as to obtain the first strain composition corresponding to the species and the absolute abundance of each strain in the first strain composition.

[0050] (6-6) For the i-th species in the final species composition, the false positive strains in the species are filtered by the optimal alignment results of each sequencing read corresponding to the species obtained in step (4-1) and the first strain composition and the absolute abundance of each strain in the first strain composition obtained in step (6-5) to obtain the strain composition of the species after the second false positive filtering.

[0051] (6-7) For the i-th species in the final species composition, based on the observed abundance of each node in the pangenome graph corresponding to the species obtained in step (6-3), perform a second PAO solution on the strain composition corresponding to the species after the second false positive filtering obtained in step (6-6) (this process is exactly the same as the first PAO solution process mentioned above, and will not be repeated here). Based on the solution results, retain strains with an absolute abundance greater than zero to obtain the second strain composition and the absolute abundance of each strain in the second strain composition. These two are respectively used as the final strain composition of the species and the absolute abundance of each strain in the final strain composition.

[0052] (6-8) Set i = i + 1 and return to step (6-2).

[0053] Preferably, for each node in the pangenome map corresponding to the species, step (6-3) involves calculating the observed abundance of node v using the following formula.

[0054]

[0055] Where v∈V, V is all nodes in the pangenome graph corresponding to this species, |v| is the number of bases in node v, and R v It refers to all sequencing reads covering this node corresponding to this species, r∈[1,R]. v ], R vr It is the r-th sequencing read covering this node for this species, |R vr | is the number of bases in the r-th sequencing read covering this node for this species;

[0056] Step (6-4) specifically involves first constructing a three-node group v from every three consecutive (sequentially appearing) nodes along the genomic path of each strain in the pan-genome graph corresponding to the i-th species obtained in step (2-1). triplet ;

[0057] Then, obtain the set V of three-node groups consisting of all three-node groups covered by all sequencing reads corresponding to the species in the pan-genome map corresponding to the species obtained in step (4-1). reads ;

[0058] Subsequently, for each strain's genomic pathway in the pan-genome map corresponding to this species, all three-node groups along that strain's genomic pathway are obtained. For each three-node group corresponding to that strain's genomic pathway, if this three-node group only appears in that strain's genomic pathway and not in other strains' genomic pathways, then this three-node group is the strain-specific three-node group corresponding to that strain's genomic pathway. All strain-specific three-node groups corresponding to that strain's genomic pathway constitute the strain-specific three-node group set V. strain ;

[0059] Subsequently, for each strain's genomic path in the pan-genome map corresponding to this species, the proportion f of the strain-specific three-node set corresponding to that genomic path appearing in all sequencing reads corresponding to this species is calculated using the following formula. strain :

[0060]

[0061] Among them, f strain ∈[0,1], where 1 indicates that the set of strain-specific three-node groups corresponding to the genomic pathway of this strain is covered by all sequencing reads corresponding to this species;

[0062] Finally, if the obtained ratio f strain If the value is less than 0.3, the strain corresponding to the genomic pathway of this strain is identified as a false positive strain and filtered out; otherwise, the strain is identified as true and is included in the strain composition after the first false positive filtering.

[0063] Preferably, the PAO objective function in step (6-5) is:

[0064]

[0065] Where V represents all nodes in the pangenome map corresponding to this species, and P represents the genomic pathways of all strains that constructed the pangenome map corresponding to this species. It is the observed abundance of node v∈V in the pangenome graph, a pIt is the absolute abundance of the strain corresponding to the strain genomic pathway p∈P, which is the decision variable that PAO needs to solve.

[0066] Step (6-6) specifically involves the following steps: First, for each strain in the first strain composition corresponding to the species obtained in step (6-5), the strain genome path corresponding to the strain in the pangenome map corresponding to the species is obtained, and used as the strain genome path corresponding to the strain in the first strain composition.

[0067] Then, for each strain's genomic path in the first strain composition corresponding to the species, the average coverage depth of each strain-specific three-node group in the strain-specific three-node group set corresponding to the genomic path of the strain is calculated based on the optimal alignment results of each sequencing read corresponding to the species obtained in step (4-1). Specifically, this process uses the following formula to calculate the strain-specific three-node group v. triplet Average coverage depth c(v) triplet ):

[0068]

[0069] Where s∈S, S is all sequencing reads covering the strain-specific three-node group corresponding to this species, |s is the number of bases in the sequencing read s covering the strain-specific three-node group corresponding to this species, and |v triplet | is the strain's specific three-node group v triplet The number of bases.

[0070] Subsequently, for each strain's genomic pathway in the first strain composition, the estimated absolute abundance of the strain's genomic pathway is calculated based on the average coverage depth of each strain-specific three-node group in the strain-specific three-node group set corresponding to that genomic pathway. Specifically, this process uses the following formula to calculate the estimated absolute abundance a of the strain's genomic pathway p. p,triplet :

[0071]

[0072] Among them, |V strain The strain-specific three-node set V corresponds to the genomic pathway p of this strain. strain The total number of strain-specific three-node groups.

[0073] Subsequently, for each strain in the first strain composition, the absolute abundance a of that strain is calculated using the following formula. p The estimated absolute abundance of strains corresponding to the strain's genomic pathway (a) p,triplet Differences between them:

[0074]

[0075] Finally, if the difference d is obtained strain If the value is greater than 0.45, the strain corresponding to the genomic pathway of this strain is identified as a false positive strain and filtered out; otherwise, the strain is identified as true and is included in the strain composition after the second false positive filtering.

[0076] Preferably, step (7) includes the following sub-steps:

[0077] (7-1) Based on the final species composition corresponding to the metagenomic data obtained in step (4) and the absolute abundance of each species in the final species composition obtained in step (5-1), the absolute abundance of each strain in the final strain composition of the metagenomic data obtained in step (6) is adjusted to obtain the adjusted absolute abundance of the strain; in this step, for each species in the final species composition corresponding to the metagenomic data obtained in step (4), if the strain with the largest absolute abundance in the species exceeds the absolute abundance corresponding to the species, the absolute abundance of all strains in the species is scaled to make the sum of the absolute abundance of all strains in the species equal to the absolute abundance of the species.

[0078] (7-2) For the final strain composition of the metagenomic data, the adjusted absolute abundance of all strains in the final strain composition obtained in step (7-1) is normalized to obtain the relative abundance of each strain in the final strain composition.

[0079] In summary, compared with the prior art, the above-described technical solutions conceived by this invention can achieve the following beneficial effects:

[0080] (1) Because the present invention employs step (2), it captures the genomes of multiple strains of each species as references through the pan-genome map technology, rather than using only the representative genome of each species as a reference; the present invention enhances the alignment accuracy of metagenomic sequencing reads and the estimation of composition at the strain level by expanding strain diversity, thus solving the technical problem that the existing methods have limited strain-level classification capabilities due to the strategy of using a single reference genome.

[0081] (2) Because the present invention employs step (3), it uses the corresponding short read alignment tool Giraffe and long read alignment tool GraphAligner respectively in the alignment of NGS and TGS metagenomic data with pan-genome maps. Therefore, it can solve the technical problem that existing methods such as MetaPhlAn4, Kaiju, Kraken2, CLARK, StrainScan, Melon, MEGAN-LR and MetaMaps cannot efficiently process both NGS and TGS metagenomic data at the same time due to the different characteristics of NGS and TGS data.

[0082] (3) Since the present invention adopts steps (4) to (7), it can obtain the classification results at the species level by comparing metagenomic data with the reference pangenome, and further optimize the path abundance based on the pangenome map to obtain the classification results at the strain level for all species. Therefore, it can solve the technical problem that existing methods such as MetaPhlAn4, Kaiju, Kraken2, CLARK, Melon, and MEGAN-LR can only classify at the species level, but it is difficult to achieve strain level classification or even multi-species strain level classification. Attached Figure Description

[0083] Figure 1 This is a schematic diagram of the method for performing strain-level analysis of metagenomic data based on pangenome maps according to the present invention;

[0084] Figure 2 These are the species-level analysis results from the zymo NGS dataset;

[0085] Figure 3 These are the results of strain-level analysis of the zymo NGS dataset;

[0086] Figure 4 This is a flowchart of the method for performing strain-level analysis of metagenomic data based on pangenome maps according to the present invention. Detailed Implementation

[0087] 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.

[0088] The basic idea of ​​this invention is to address the limitations of existing methods in data processing and strain-level classification from three aspects. First, it captures the genomes of multiple strains of each species as references using pan-genome mapping technology, rather than just using a representative genome of each species. This expansion of strain diversity enhances the alignment accuracy of metagenomic sequencing reads and the estimation of strain-level composition. Second, by using the corresponding short-read alignment tool Giraffe and long-read alignment tool GraphAligner in the alignment of NGS and TGS metagenomic data with the pan-genome map, respectively, both types of data can be processed efficiently simultaneously. Finally, by further optimizing the path abundance of the species-level classification results based on the pan-genome map, multi-species strain-level classification can be achieved.

[0089] like Figure 1 and Figure 4 As shown, this invention provides a method for performing strain-level analysis of metagenomic data based on pangenome maps, comprising the following steps:

[0090] (1) Obtain multiple species and their corresponding complete genomes from public databases, and construct a reference database based on all the obtained complete genomes.

[0091] This step is a data preprocessing step and is optional for the user. Users can also specify the genome to be used to build the reference database.

[0092] This step includes the following sub-steps:

[0093] (1-1) Download and obtain multiple species and multiple complete genomes corresponding to each species from public databases;

[0094] (1-2) Remove plasmids from all the complete genomes obtained in step (1-1) to obtain multiple complete genomes after plasmid removal;

[0095] (1-3) For each species obtained in step (1-1), the average nucleotide similarity (ANI) between any pair of complete genomes corresponding to that species is calculated using the fastANI tool, and is used as the distance score between that pair of complete genomes corresponding to that species.

[0096] (1-4) For each species obtained in step (1-1), a heuristic clustering algorithm based on connected graphs is used to cluster the distance scores between all pairs of complete genomes corresponding to the species obtained in step (1-3) to obtain multiple complete genome sets. The complete genome set containing the most complete genomes is selected from all complete genome sets, and a predetermined number of complete genomes (the predetermined number is less than or equal to 10, preferably 10; if the number of complete genomes in the complete genome set is less than the predetermined number, the actual number is taken) are selected from the selected complete genome set as multiple reference genomes corresponding to the species. All reference genomes corresponding to all species constitute a reference database.

[0097] Specifically, in this step, if the clustering result is less than a preset first clustering threshold (which is 95%), it indicates that the pair of complete genomes corresponding to this clustering result are strains from different species; if the clustering result is greater than or equal to the preset first clustering threshold and less than a preset second clustering threshold (which is 99.9%), it indicates that the pair of complete genomes corresponding to this clustering result are different strains from the same species; if the clustering result is greater than the preset second clustering threshold, it indicates that the pair of complete genomes corresponding to this clustering result are from the same strain (the difference in genome structure is extremely small and can be almost ignored). This step's clustering process obtains clustering results greater than or equal to the preset first clustering threshold and less than the preset second clustering threshold as multiple sets of complete genomes obtained through clustering.

[0098] (2) For each species in the reference database obtained in step (1), construct the pangenome map corresponding to the species using all the reference genomes corresponding to the species in the reference database, and merge the pangenome maps corresponding to all species to obtain the reference pangenome map.

[0099] This step includes the following sub-steps:

[0100] (2-1) For each species in the reference database obtained in step (1), if the species has more than one reference genome in the reference database, the PGGB (PanGenome Graph Builder) tool is used to construct the pan-genome graph corresponding to the species for all reference genomes. If the species has only one reference genome in the reference database, the reference genome is divided into fragments of 1024 bp in length, and each fragment is used as a node to construct the pan-genome graph, thereby obtaining the pan-genome graph corresponding to the species.

[0101] In this step, the pangenome graph is a variation graph, denoted as G = (V, E, P), where V are nodes, E are edges, and P represents the genomic paths of all strains formed by all paths in the pangenome graph G. G is a bidirectional sequence graph, where each node v ∈ V represents a base sequence seq(v) of indefinite length; each edge e ij =(v i ,v j )∈E represents the node v i and v j The sequences seq(v) represent respectively. i ) and seq(v j The edge appears as a continuous subsequence in a certain strain genome path that constitutes the graph; each edge is contained in the strain genome path P. Let... Represents node v i If the complementary sequence is e, then e ij reverse edge Defined as In a mutation graph, nodes and edges can be traversed in either forward or reverse direction. Each strain's genome, as part of the definition of the mutation graph G, is represented as a specific strain genome path p = (v1,...,v...). l )∈P, corresponding to a strain. If a strain genome corresponds to only one path in the pangenome map, then that path is the strain genome path p; if a strain genome corresponds to multiple paths in the pangenome map, then all paths of the strain genome are taken as the strain genome path p.

[0102] The advantage of this step (2-1) is that the pangenome mapping technology can capture the genomes of multiple strains of each species as references, rather than just using the representative genome of each species as a reference.

[0103] (2-2) Merge the pangenome maps corresponding to all species obtained in step (2-1) (the node order needs to be re-labeled during the merging process, and no new nodes and edges are introduced) to obtain a reference pangenome map.

[0104] The advantage of this step (2-2) is that by comparing metagenomic data with the pangenome map once, the alignment results of sequencing reads of all species in the metagenomic data can be obtained.

[0105] (3) Obtain the metagenomic data input by the user, and compare the reference pangenome map obtained in step (2) with the metagenomic data to obtain multiple alignment results for each sequencing read in the metagenomic data.

[0106] In this step, for second-generation short-read metagenomic data, the Giraffe tool is used for alignment. Specifically, the vg tool is first used to build an index for the reference pan-genome map to accelerate the alignment process. Then, the Giraffe tool is used to align the short-read metagenomic data with the reference pan-genome map to obtain multiple alignment results for each sequencing read in the short-read metagenomic data.

[0107] In this step, for third-generation long-read metagenomic data, the GraphAligner tool is used to align the two. This tool does not require indexing and directly aligns the long-read metagenomic data with the reference pangenome map to obtain multiple alignment results for each sequencing read in the long-read metagenomic data.

[0108] The advantage of this step (3) is that the corresponding short-read alignment tool Giraffe and long-read alignment tool GraphAligner are used for second-generation short-read and third-generation long-read metagenomic data respectively, so that the present invention can process two types of metagenomic data simultaneously and efficiently.

[0109] (4) Perform species-level taxonomic binning based on all alignment results of each sequencing read in the metagenomic data obtained in step (3) to obtain the final species composition corresponding to the metagenomic data.

[0110] This step includes the following sub-steps:

[0111] (4-1) For each sequencing read in the metagenomic data obtained in step (3), the optimal alignment result of the sequencing read is obtained based on all alignment results of the sequencing read.

[0112] In this step, if a sequencing read has multiple alignment results, only the best alignment result is retained, thus ensuring that each sequencing read corresponds to one species.

[0113] Specifically, in short read alignments, each sequencing read has only one alignment result, so this alignment result is the best alignment result; while in long read alignments, the best alignment result is determined from all alignment results by using the Mapping Quality (MAPQ) value, the number of matching bases, identity, and sequencing read length.

[0114] (4-2) For each sequencing read in the metagenomic data obtained in step (3), the species corresponding to the sequencing read is obtained according to the best alignment result of the sequencing read obtained in step (4-1). The species corresponding to all sequencing reads in the metagenomic data constitute the first species composition corresponding to the metagenomic data (which has potential false positives).

[0115] (4-3) Perform false positive species filtering on the first species composition corresponding to the metagenomic data obtained in step (4-2) to obtain the second species composition as the final species composition corresponding to the metagenomic data.

[0116] Specifically, for each species in the first species composition, the MAPQ values ​​of multiple sequencing reads corresponding to that species are obtained from the best alignment results of multiple sequencing reads corresponding to that species. Among all the sequencing reads corresponding to that species, if at least one sequencing read has a MAPQ value that reaches the highest score of 60 (this indicates that at least one sequencing read corresponding to that species is statistically significant to match that species), and at least one-tenth of the sequencing reads have a MAPQ value greater than 2, then the species is identified as true; otherwise, the species is identified as a false positive and then filtered out.

[0117] (5) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), perform species-level taxonomic profiling to obtain the relative abundance of the species.

[0118] This step includes the following sub-steps:

[0119] (5-1) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), the absolute abundance of the species is obtained based on the best alignment result of each sequencing read corresponding to the species obtained in step (4-1).

[0120] This step uses the following formula to calculate the absolute abundance c of the i-th species in the final species composition. i :

[0121]

[0122] Where i∈[1, the total number of species in the final species composition T], n i It is the total number of sequencing reads corresponding to the i-th species, j∈[1, n] i ], n ij It is the j-th sequencing read corresponding to the i-th species, |n ij | is the number of bases in the j-th sequencing read corresponding to the i-th species, m iIt is the total number of reference genomes corresponding to the i-th species, k∈[1, m] i ], g ik It is the k-th reference genome corresponding to the i-th species, |g ik | is the number of bases in the k-th reference genome corresponding to the i-th species.

[0123] (5-2) Normalize the absolute abundance of each species in the final species composition obtained in step (5-1) to obtain the relative abundance of that species.

[0124] This step specifically uses the following formula to calculate the relative abundance a of the i-th species in the final species composition. i :

[0125]

[0126] (6) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), the species is subjected to path abundance optimization (PAO) processing according to the pangenome map corresponding to the species obtained in step (2) to obtain the final strain composition of the metagenomic data and the absolute abundance of each strain in the final strain composition.

[0127] This step includes the following sub-steps:

[0128] (6-1) Set the counter i = 1 (where i is the species pointer in the final species composition corresponding to the metagenomic data, and i = 1 refers to the first species in the final species composition corresponding to the metagenomic data).

[0129] (6-2) Determine whether i is greater than the total number of species T corresponding to the final species composition of the metagenomic data. If so, the process ends; otherwise, proceed to step (6-3).

[0130] (6-3) For the i-th species in the final species composition, calculate the observed abundance of each node in the pangenome graph corresponding to the species based on the best alignment results of each sequencing read corresponding to the species obtained in step (4-1).

[0131] For each node in the pangenome map corresponding to this species, this step uses the following formula to calculate the observed abundance of node v.

[0132]

[0133] Where v∈V, V is all nodes in the pangenome graph corresponding to this species, |v| is the number of bases in node v, and R vIt refers to all sequencing reads covering this node corresponding to this species, r∈[1,R]. v ], R vr It is the r-th sequencing read covering this node for this species, |R vr | represents the number of bases in the r-th sequencing read covering this node for this species.

[0134] (6-4) For the i-th species in the final species composition, the false positive strains in the species are filtered according to the pan-genome map of the species obtained in step (2-1) and the optimal alignment results of each sequencing read of the species obtained in step (4-1) to obtain the strain composition of the species after the first false positive filtering.

[0135] Specifically, this step involves first constructing a triplet node from every three consecutive (sequentially appearing) nodes along the genomic pathway of each strain in the pan-genome graph corresponding to the i-th species obtained in step (2-1). triplet ;

[0136] Then, obtain the set V of three-node groups consisting of all three-node groups covered by all sequencing reads corresponding to the species in the pan-genome map corresponding to the species obtained in step (4-1). reads ;

[0137] Subsequently, for each strain's genomic pathway in the pan-genome map corresponding to this species, all three-node groups along that strain's genomic pathway are obtained. For each three-node group corresponding to that strain's genomic pathway, if this three-node group only appears in that strain's genomic pathway and not in other strains' genomic pathways, then this three-node group is the strain-specific three-node group corresponding to that strain's genomic pathway. All strain-specific three-node groups corresponding to that strain's genomic pathway constitute the strain-specific three-node group set V. strain ;

[0138] Subsequently, for each strain's genomic path in the pan-genome map corresponding to this species, the proportion f of the strain-specific three-node set corresponding to that genomic path appearing in all sequencing reads corresponding to this species is calculated using the following formula. strain :

[0139]

[0140] Among them, f strain ∈[0,1], where 1 indicates that the set of strain-specific three-node groups corresponding to the genomic pathway of this strain is covered by all sequencing reads corresponding to this species;

[0141] Finally, if the obtained ratio f strainIf the value is less than 0.3, the strain corresponding to the genomic pathway of this strain is identified as a false positive strain and filtered out; otherwise, the strain is identified as true and is included in the strain composition after the first false positive filtering.

[0142] (6-5) For the i-th species in the final species composition, based on the observed abundance of each node in the pangenome graph corresponding to the species obtained in step (6-3), perform the first PAO solution on the strain composition after the first false positive filtering obtained in step (6-4) for the species, and retain the strains with absolute abundance greater than zero according to the solution results, so as to obtain the first strain composition corresponding to the species and the absolute abundance of each strain in the first strain composition.

[0143] In this step, PAO is a linear programming problem. It predicts the absolute abundance of each strain in a species by minimizing the difference between the observed abundance of each node in the pangenome graph corresponding to that node and the absolute abundance of the strain corresponding to that node. The PAO objective function is defined as follows:

[0144]

[0145] Where V represents all nodes in the pangenome map corresponding to this species, and P represents the genomic pathways of all strains that constructed the pangenome map corresponding to this species. It is the observed abundance of node v∈V in the pangenome graph, a p It is the absolute abundance of the strain corresponding to the strain genomic pathway p∈P, which is the decision variable that PAO needs to solve.

[0146] (6-6) For the i-th species in the final species composition, the false positive strains in the species are filtered by the optimal alignment results of each sequencing read corresponding to the species obtained in step (4-1) and the first strain composition and the absolute abundance of each strain in the first strain composition obtained in step (6-5) to obtain the strain composition of the species after the second false positive filtering.

[0147] Specifically, this step involves firstly obtaining the strain genome path corresponding to the strain in the first strain composition of the species obtained in step (6-5), and using it as the strain genome path corresponding to the strain in the pangenome map of the species.

[0148] Then, for each strain in the first strain composition corresponding to the species, the average coverage depth of each strain-specific three-node group in the strain-specific three-node group set corresponding to the species is calculated based on the best alignment result of each sequencing read corresponding to the species obtained in step (4-1).

[0149] Specifically, this process uses the following formula to calculate the strain-specific three-node group v. triplet Average coverage depth c(v) triplet ):

[0150]

[0151] Where s∈S, S is all sequencing reads covering the strain-specific three-node group corresponding to this species, |s is the number of bases in the sequencing read s covering the strain-specific three-node group corresponding to this species, and |v triplet | is the strain's specific three-node group v triplet The number of bases.

[0152] Subsequently, for each strain in the first strain composition, the estimated absolute abundance of strains in the strain genome path is calculated based on the average coverage depth of each strain-specific three-node group in the strain-specific three-node group set corresponding to the strain genome path.

[0153] Specifically, this process uses the following formula to calculate the estimated absolute abundance a of the strain's genomic pathway p. p,triplet :

[0154]

[0155] Among them, |V strain The strain-specific three-node set V corresponds to the genomic pathway p of this strain. strain The total number of strain-specific three-node groups.

[0156] Subsequently, for each strain in the first strain composition, the absolute abundance a of that strain is calculated using the following formula. p The estimated absolute abundance of strains corresponding to the strain's genomic pathway (a) p,triplet Differences between them:

[0157]

[0158] Finally, if the difference d is obtained strain If the value is greater than 0.45, the strain corresponding to the genomic pathway of this strain is identified as a false positive strain and filtered out; otherwise, the strain is identified as true and is included in the strain composition after the second false positive filtering.

[0159] (6-7) For the i-th species in the final species composition, based on the observed abundance of each node in the pangenome graph corresponding to the species obtained in step (6-3), perform a second PAO solution on the strain composition corresponding to the species after the second false positive filtering obtained in step (6-6) (this process is exactly the same as the first PAO solution process mentioned above, and will not be repeated here). Based on the solution results, retain strains with an absolute abundance greater than zero to obtain the second strain composition and the absolute abundance of each strain in the second strain composition. These two are respectively used as the final strain composition of the species and the absolute abundance of each strain in the final strain composition.

[0160] (6-8) Set i = i + 1 and return to step (6-2);

[0161] The advantage of this step (6) is that, for each species, PAO treatment in the pan-genome can accurately resolve the strain composition and relative abundance of strains.

[0162] (7) Based on the final species composition corresponding to the metagenomic data obtained in step (4), perform abundance adjustment and normalization calculation on the final strain composition of the metagenomic data obtained in step (6) and the absolute abundance of each strain in the final strain composition to obtain the final strain composition of the metagenomic data and the relative abundance of each strain in the final strain composition.

[0163] This step includes the following sub-steps:

[0164] (7-1) Based on the final species composition corresponding to the metagenomic data obtained in step (4) and the absolute abundance of each species in the final species composition obtained in step (5-1), the absolute abundance of each strain in the final strain composition of the metagenomic data obtained in step (6) is adjusted to obtain the adjusted absolute abundance of the strain.

[0165] In this step, for each species in the final species composition corresponding to the metagenomic data obtained in step (4), if the strain with the highest absolute abundance in the species exceeds the absolute abundance of the species, the absolute abundance of all strains in the species is scaled so that the sum of the absolute abundance of all strains in the species is equal to the absolute abundance of the species.

[0166] (7-2) For the final strain composition of the metagenomic data, the adjusted absolute abundance of all strains in the final strain composition obtained in step (7-1) is normalized to obtain the relative abundance of each strain in the final strain composition.

[0167] This invention is a method for classifying metagenomic data at the strain level based on pangenome maps. Its basic task is to perform taxonomic classification at different taxonomic levels, including species and strain levels. The input is metagenomic sequencing reads (including short and long reads), and the output is the composition and relative abundance of the metagenomic genome at different levels.

[0168] After performing all the above steps, the final species-level classification result of the present invention is: the final species composition in step (4) and the relative abundance of each species in the final species composition in step (5);

[0169] The final strain-level classification result of this invention is: the final strain composition of step (7) and the relative abundance of each strain in the final strain composition.

[0170] It is worth noting that this invention focuses on solving classification problems at the metagenomic species and strain levels, as these problems remain quite challenging. Theoretically, this invention is also applicable to classification at levels above the species level. Furthermore, strain-level classification (i.e., steps (6) and (7)) is an optional analysis, implemented based on species-level classification.

[0171] 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:

[0172] (1) Precision and Recall

[0173] 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 data.

[0174] The calculation formula is as follows:

[0175]

[0176] 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 data but is predicted by the method, and FN represents a taxonomic unit that exists in the metagenomic data but is not predicted.

[0177] (2)AUPR

[0178] Plot the curves of precision and recall as a function of species abundance thresholds and calculate the area under the curve (AUPR).

[0179] For methods that cannot recall all actually present taxa in the metagenomic data, the portion from their highest recall to 100% recall will be assigned an AUPR score of 0. Conversely, methods that achieve 100% recall will not be further penalized in their AUPR score for additional false positive taxa.

[0180] (3) L2 distance

[0181] The calculation formula is as follows:

[0182]

[0183] Here, U represents the intersection of the taxonomic units predicted by the method and the actual taxonomic units existing in the metagenomic data. For a certain taxonomic unit u in U, it is represented by... and α i Let i represent the predicted relative abundance and the actual relative abundance of classification unit i, respectively.

[0184] Example 1

[0185] The experimental data were prepared as follows: sim-low NGS dataset. Sixty complete genomes from the publicly available NCBI RefSeq database were selected, representing 30 species with two strains per species. The CAMISIM tool was used to simulate 2G of second-generation short-read NGS data.

[0186] 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 and relative abundance of each species in the metagenomic data of the present invention, and the strain composition and relative abundance of each strain in the metagenomic data.

[0187] For each of the other methods, the experimental steps of this method are as follows:

[0188] (1) Obtain the reference database constructed in this example according to the present invention, and perform data preprocessing to obtain the reference database for species-level experiments and strain-level experiments of the method.

[0189] In this step, the data preprocessing is as follows: For species-level experiments, for each species, the longest complete genome of that species is selected from the reference database constructed in this example as the reference genome, and then the reference database for species-level experiments is constructed. Note that MetaPhlAn4 uses the reference database constructed by this method as the reference database for species-level experiments. For strain-level experiments, the reference database constructed in this example is directly used as the reference database for strain-level experiments.

[0190] (2) Obtain the metagenomic data input by the user, use the reference database for species-level experiments, and perform species-level experiments according to the specific implementation method of this method. If the method outputs sequencing read alignment results, perform species-level taxonomic analysis based on the sequencing read alignment results to obtain the species composition corresponding to the metagenomic data and the abundance of each species in the species composition. Otherwise, use the species composition corresponding to the metagenomic data and the relative abundance of each species in the species composition directly output by this method.

[0191] In this step, the species-level taxonomic analysis is specifically performed as follows: First, for each sequencing read alignment result of this method, the species corresponding to that sequencing read is obtained. The species corresponding to all sequencing reads in the metagenomic data constitute the species composition of the metagenomic data. Then, for each species in the species composition, all sequencing read alignment results corresponding to that species are obtained, and the relative abundance of that species is calculated using step (5) of this invention.

[0192] (3) If the method can perform experiments at the strain level, obtain the metagenomic data input by the user, use the reference database for experiments at the strain level, and perform experiments at the strain level according to the specific implementation of the method, if the method outputs sequencing alignment results, then perform strain-level taxonomic analysis based on the sequencing alignment results to obtain the strain composition corresponding to the metagenomic data and the relative abundance of each strain in the strain composition.

[0193] In this step, the process for strain-level taxonomic analysis is the same as that for species-level taxonomic analysis, except that species are replaced with strains, which will not be described again here.

[0194] Table 1-1: Species-level experimental results of the sim-low NGS dataset

[0195]

[0196] Table 1-2: Experimental results at the strain level in the sim-low NGS dataset

[0197]

[0198]

[0199] In this experiment, the best results for each metric are indicated in bold. At the species level, PanTax of this invention performed best in recall, AUPR, and L2 distance. In precision, it was only slightly lower than MetaPhlAn4, but still significantly better than other methods. At the strain level, PanTax significantly outperformed other methods in precision and L2 distance, while performing comparably in recall and AUPR.

[0200] Example 2

[0201] The experimental data were prepared as follows: sim-low ONT R10.4 dataset. Sixty complete genomes from the publicly available NCBIRefSeq database were selected, representing 30 species with two strains per species. The CAMISIM tool was used to simulate 2G of three generations of long-read ONT R10.4 data.

[0202] 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 and relative abundance of each species in the metagenomic data of the present invention, and the strain composition and relative abundance of each strain in the metagenomic data.

[0203] 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 data of other methods, as well as the strain composition and relative abundance of each strain in the metagenomic data.

[0204] Table 2-1: Species-level experimental results on the sim-low ONT R10.4 dataset

[0205]

[0206] Table 2-2: Experimental results at the strain level in the sim-low ONT R10.4 dataset

[0207]

[0208] In this experiment, the best results for each metric are indicated in bold. At the species level, PanTax demonstrated the best performance in precision, recall, AUPR, and L2 distance. At the strain level, PanTax significantly outperformed other methods in precision, recall, and L2 distance, while its AUPR performance was only slightly weaker than Centrifuge.

[0209] Example 3

[0210] The experimental data was prepared as follows: The sim-high HiFi dataset consisted of 1000 complete genomes selected from the publicly available NCBI RefSeq database. This included 373 species, with 1–5 strains per species. The CAMISIM tool was used to simulate 20GB of three generations of long-read HiFi data.

[0211] 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 and relative abundance of each species in the metagenomic data of the present invention, and the strain composition and relative abundance of each strain in the metagenomic data.

[0212] 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 data of other methods, as well as the strain composition and relative abundance of each strain in the metagenomic data.

[0213] Table 3-1: Species-level experimental results of the sim-high HiFi dataset

[0214]

[0215] Table 3-2: Experimental results at the strain level in the sim-high HiFi dataset

[0216]

[0217] In this experiment, the best results for each metric are indicated in bold. This experiment demonstrates the performance of each method on complex datasets. At the species level, PanTax of this invention shows the best performance in precision, recall, AUPR, and L2 distance. At the strain level, PanTax shows the best performance in precision; its L2 distance is comparable to Centrifuge and significantly better than MetaMaps; and its recall and AUPR performance are only slightly weaker than Centrifuge but better than MetaMaps.

[0218] Example 4

[0219] The experimental data was prepared as follows: the Zymo NGS dataset. This dataset is Illumina data from the ZymoBIOMICS microbial community standard dataset (catalog number: D6330), and the data type is second-generation short-read data. This simulated community contains 10 microorganisms: 8 bacteria and 2 yeasts. Since this experiment mainly focuses on the results at the bacterial strain level, the two yeast species were ignored in the implementation.

[0220] The specific implementation steps of this embodiment are performed according to the specific implementation method of the present invention, with only step (1) being adjusted. The adjusted step (1) is as follows:

[0221] (1) Obtain multiple species and their corresponding complete genomes from public databases, and construct a reference database based on all the obtained complete genomes.

[0222] This step includes the following sub-steps:

[0223] (1-1) Given the known species present in the dataset, for each species present in the dataset, download and obtain multiple complete genomes corresponding to that species from a public database;

[0224] (1-2) Remove plasmids from all the complete genomes obtained in step (1-1) to obtain multiple complete genomes after plasmid removal;

[0225] (1-3) For each species obtained in step (1-1), the average nucleotide similarity (ANI) between any pair of complete genomes corresponding to that species is calculated using the fastANI tool, and is used as the distance score between that pair of complete genomes corresponding to that species.

[0226] (1-4) For each species obtained in step (1-1), the distance scores between all pairs of complete genomes corresponding to the species obtained in step (1-3) are clustered using a single-linkage hierarchical clustering algorithm to obtain a complete genome set. A predetermined number of complete genomes (preferred number is less than or equal to 80, preferably 80; if the number of complete genomes in the complete genome set is less than the predetermined number, the actual number is used) are selected from the selected complete genome set as multiple reference genomes corresponding to the species. All reference genomes corresponding to all species constitute a reference database.

[0227] Specifically, this step uses clustering to obtain clustering results with less than 99.9% of the data, in order to obtain multiple complete genome sets.

[0228] To better illustrate the results of this embodiment, a simple Python script was used to visualize the results obtained by this invention and compare them with the distribution of real species and strains in the dataset, such as... Figure 2 and Figure 3 As shown.

[0229] 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 data of other methods, as well as the strain composition and relative abundance of each strain in the metagenomic data.

[0230] Table 4: Experimental results at the strain level in the zymo NGS dataset

[0231] method accuracy Recall rate AUPR L2 distance PanTax 0.875 0.875 0.733 0.202 Centrifuge 0.014 0.875 0.858 0.494

[0232] This embodiment correctly predicted strains from 7 out of 8 species (one strain from each species). The strain from the remaining species was not in the reference database and therefore could not be found; however, a very similar strain was found, with an ANI of 99.9284%, indicating that the two strains were essentially identical. At the strain level, PanTax significantly outperformed Centrifuge in precision and L2 distance, while showing similar recall. However, it underperformed Centrifuge in AUPR because, as mentioned earlier, the abundance of strains not in the reference database was ranked higher, which did not actually affect the results.

[0233] These benchmark results demonstrate that PanTax significantly outperforms existing state-of-the-art methods, primarily due to its significantly higher precision and recall (at both the species and strain levels), while maintaining comparable or better performance in other areas. Furthermore, PanTax can efficiently process both second- and third-generation metagenomic data simultaneously.

[0234] 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 method for classifying metagenomic data at the strain level based on pangenome maps, characterized in that, Includes the following steps: (1) Obtain multiple species and their corresponding complete genomes from public databases, and construct a reference database based on all the obtained complete genomes; step (1) includes the following sub-steps: (1-1) Download and obtain multiple species and multiple complete genomes corresponding to each species from public databases; (1-2) Remove plasmids from all the complete genomes obtained in step (1-1) to obtain multiple complete genomes after plasmid removal; (1-3) For each species obtained in step (1-1), the average nucleotide similarity (ANI) between any pair of complete genomes corresponding to that species is calculated using the fastANI tool, and is used as the distance score between that pair of complete genomes corresponding to that species. (1-4) For each species obtained in step (1-1), a heuristic clustering algorithm based on connected graphs is used to cluster the distance scores between all pairs of complete genomes corresponding to the species obtained in step (1-3) to obtain multiple complete genome sets. The complete genome set containing the most complete genomes is selected from all complete genome sets, and a predetermined number of complete genomes are selected from the selected complete genome set as multiple reference genomes corresponding to the species. All reference genomes corresponding to all species constitute a reference database. In steps (1-4), if the clustering result is less than the preset first clustering threshold, it means that the pair of complete genomes corresponding to the clustering result are strains of different species. If the clustering result is greater than or equal to the preset first clustering threshold and less than the preset second clustering threshold, it means that the pair of complete genomes corresponding to the clustering result are different strains of the same species. If the clustering result is greater than the preset second clustering threshold, it means that the pair of complete genomes corresponding to the clustering result are from the same strain; The clustering process in steps (1-4) involves obtaining clustering results that are greater than or equal to the preset first clustering threshold and less than the preset second clustering threshold as multiple complete genome sets obtained from the clustering process. (2) For each species in the reference database obtained in step (1), construct the pangenome map corresponding to the species using all the reference genomes corresponding to the species in the reference database, and merge the pangenome maps corresponding to all species to obtain the reference pangenome map; (3) Obtain the metagenomic data input by the user, and compare the reference pangenome map obtained in step (2) with the metagenomic data to obtain multiple alignment results for each sequencing read in the metagenomic data; (4) Based on all alignment results of each sequencing read in the metagenomic data obtained in step (3), perform species-level classification and binning to obtain the final species composition corresponding to the metagenomic data; (5) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), perform species-level taxonomic spectrum analysis on the species to obtain the relative abundance of the species. (6) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), the species is subjected to path abundance optimization (PAO) processing according to the pangenome map corresponding to the species obtained in step (2) to obtain the final strain composition of the metagenomic data and the absolute abundance of each strain in the final strain composition. (7) Based on the final species composition corresponding to the metagenomic data obtained in step (4), perform abundance adjustment and normalization calculation on the final strain composition of the metagenomic data obtained in step (6) and the absolute abundance of each strain in the final strain composition to obtain the final strain composition of the metagenomic data and the relative abundance of each strain in the final strain composition.

2. The method for classifying metagenomic data at the strain level based on pangenome maps according to claim 1, characterized in that, Step (2) includes the following sub-steps: (2-1) For each species in the reference database obtained in step (1), if the species has more than one reference genome in the reference database, the PGGB tool is used to construct the pan-genome map corresponding to the species for all reference genomes; if the species has only one reference genome in the reference database, the reference genome is divided into fragments of 1024 bp in length, and each fragment is used as a node to construct the pan-genome map, thereby obtaining the pan-genome map corresponding to the species. (2-2) Merge the pangenome maps of all species obtained in step (2-1) to obtain a reference pangenome map; In step (3), for second-generation short-read metagenomic data, the Giraffe tool is used to compare the two; specifically, the vg tool is first used to build an index for the reference pan-genome map to accelerate the comparison process, and then the Giraffe tool is used to compare the short-read metagenomic data with the reference pan-genome map to obtain multiple comparison results for each sequencing read in the short-read metagenomic data; In step (3), for third-generation long-read metagenomic data, the GraphAligner tool is used to compare the two to obtain multiple alignment results for each sequencing read in the long-read metagenomic data; Step (4) includes the following sub-steps: (4-1) For each sequencing read in the metagenomic data obtained in step (3), the optimal alignment result of the sequencing read is obtained based on all alignment results of the sequencing read; (4-2) For each sequencing read in the metagenomic data obtained in step (3), the species corresponding to the sequencing read is obtained according to the best alignment result of the sequencing read obtained in step (4-1). The species corresponding to all sequencing reads in the metagenomic data constitute the first species composition corresponding to the metagenomic data. The first species composition has potential false positives. (4-3) Perform false positive species filtering on the first species composition corresponding to the metagenomic data obtained in step (4-2) to obtain the second species composition as the final species composition corresponding to the metagenomic data.

3. The method for classifying metagenomic data at the strain level based on pangenome maps according to claim 2, characterized in that, In step (4-1), if a sequencing read has multiple alignment results, only the best alignment result is retained to ensure that each sequencing read corresponds to one species. Specifically, in the alignment results of short reads, each sequencing read has only one alignment result, so this alignment result is the best alignment result. In the alignment results of long reads, the best alignment result is determined from all alignment results by the alignment quality MAPQ value, the number of matching bases, identity, and sequencing read length. Step (4-3) specifically involves obtaining the MAPQ values ​​of multiple sequencing reads corresponding to each species in the first species composition from the best alignment results of multiple sequencing reads corresponding to that species. If at least one sequencing read has a MAPQ value that reaches the highest score of 60, and at least one-tenth of the sequencing reads have a MAPQ value greater than 2, then the species is identified as true; otherwise, the species is identified as a false positive and then filtered.

4. The method for classifying metagenomic data at the strain level based on pangenome maps according to claim 3, characterized in that, Step (5) includes the following sub-steps: (5-1) For each species in the final species composition corresponding to the metagenomic data obtained in step (4), the absolute abundance of the species is obtained based on the best alignment result of each sequencing read corresponding to the species obtained in step (4-1). This step uses the following formula to calculate the absolute abundance c of the i-th species in the final species composition. i : Where i∈[1, the total number of species in the final species composition T], n i It is the total number of sequencing reads corresponding to the i-th species, j∈[1, n] i ], n ij It is the j-th sequencing read corresponding to the i-th species, |n ij | is the number of bases in the j-th sequencing read corresponding to the i-th species, m i It is the total number of reference genomes corresponding to the i-th species, k∈[1, m] i ], g ik It is the k-th reference genome corresponding to the i-th species, |g ik | is the number of bases in the k-th reference genome corresponding to the i-th species; (5-2) Normalize the absolute abundance of each species in the final species composition obtained in step (5-1) to obtain the relative abundance of that species. This step specifically uses the following formula to calculate the relative abundance a of the i-th species in the final species composition. i :

5. The method for classifying metagenomic data at the strain level based on a pangenome map according to claim 4, characterized in that, Step (6) includes the following sub-steps: (6-1) Set the counter i = 1; (6-2) Determine whether i is greater than the total number of species T corresponding to the final species composition of the metagenomic data. If so, the process ends; otherwise, proceed to step (6-3). (6-3) For the i-th species in the final species composition, calculate the observed abundance of each node in the pangenome graph corresponding to the species based on the best alignment results of each sequencing read corresponding to the species obtained in step (4-1); (6-4) For the i-th species in the final species composition, the false positive strains in the species are filtered according to the pan-genome map of the species obtained in step (2-1) and the best alignment results of each sequencing read of the species obtained in step (4-1) to obtain the strain composition of the species after the first false positive filtering. (6-5) For the i-th species in the final species composition, based on the observed abundance of each node in the pangenome graph corresponding to the species obtained in step (6-3), perform the first PAO solution on the strain composition corresponding to the species after the first false positive filtering obtained in step (6-4), and retain the strains with absolute abundance greater than zero according to the solution results, so as to obtain the first strain composition corresponding to the species and the absolute abundance of each strain in the first strain composition. (6-6) For the i-th species in the final species composition, the false positive strains in the species are filtered by the optimal alignment results of each sequencing read corresponding to the species obtained in step (4-1) and the first strain composition and the absolute abundance of each strain in the first strain composition obtained in step (6-5) to obtain the strain composition of the species after the second false positive filtering. (6-7) For the i-th species in the final species composition, based on the observed abundance of each node in the pangenome graph corresponding to the species obtained in step (6-3), perform a second PAO solution on the strain composition corresponding to the species after the second false positive filtering obtained in step (6-6), and retain strains with absolute abundance greater than zero according to the solution results to obtain the second strain composition and the absolute abundance of each strain in the second strain composition. The two are respectively used as the final strain composition of the species and the absolute abundance of each strain in the final strain composition. (6-8) Set i = i + 1 and return to step (6-2).

6. The method for classifying metagenomic data at the strain level based on pangenome maps according to claim 5, characterized in that, For each node in the pangenome map corresponding to this species, step (6-3) involves calculating the observed abundance of node v using the following formula. Where v∈V, V is all nodes in the pangenome graph corresponding to this species, |v| is the number of bases in node v, and R v It refers to all sequencing reads covering this node corresponding to this species, r∈[1,R]. v ], R vr It is the r-th sequencing read covering this node for this species, |R vr | is the number of bases in the r-th sequencing read covering this node for this species; Step (6-4) specifically involves first constructing a three-node group v from every three consecutive nodes along the genomic pathway of each strain in the pan-genome graph corresponding to the i-th species obtained in step (2-1). triplet ; Then, obtain the set V of three-node groups consisting of all three-node groups covered by all sequencing reads corresponding to the species in the pan-genome map corresponding to the species obtained in step (4-1). reads ; Subsequently, for each strain's genomic pathway in the pan-genome map corresponding to this species, all three-node groups along that strain's genomic pathway are obtained. For each three-node group corresponding to that strain's genomic pathway, if this three-node group only appears in that strain's genomic pathway and not in other strains' genomic pathways, then this three-node group is the strain-specific three-node group corresponding to that strain's genomic pathway. All strain-specific three-node groups corresponding to that strain's genomic pathway constitute the strain-specific three-node group set V. strain ; Subsequently, for each strain's genomic path in the pan-genome map corresponding to this species, the proportion f of the strain-specific three-node set corresponding to that genomic path appearing in all sequencing reads corresponding to this species is calculated using the following formula. strain : Among them, f strain ∈[0,1],1] indicates that the set of strain-specific three-node groups corresponding to the genomic pathway of this strain is covered by all sequencing reads corresponding to this species; Finally, if the obtained ratio f strain If the value is less than 0.3, the strain corresponding to the genomic pathway of this strain is identified as a false positive strain and filtered out; otherwise, the strain is identified as true and is included in the strain composition after the first false positive filtering.

7. The method for classifying metagenomic data at the strain level based on a pangenome map according to claim 6, characterized in that, The objective function for PAO in step (6-5) is: Where V represents all nodes in the pangenome map corresponding to this species, and P represents the genomic pathways of all strains that constructed the pangenome map corresponding to this species. It is the observed abundance of node v∈V in the pangenome graph, a p It is the absolute abundance of the strain corresponding to the strain genomic pathway p∈P, which is the decision variable that PAO needs to solve; Step (6-6) specifically involves the following steps: First, for each strain in the first strain composition corresponding to the species obtained in step (6-5), the strain genome path corresponding to the strain in the pangenome map corresponding to the species is obtained, and used as the strain genome path corresponding to the strain in the first strain composition. Then, for each strain's genomic path in the first strain composition corresponding to the species, the average coverage depth of each strain-specific three-node group in the strain-specific three-node group set corresponding to the genomic path of the strain is calculated based on the optimal alignment results of each sequencing read corresponding to the species obtained in step (4-1). Specifically, this process uses the following formula to calculate the strain-specific three-node group v. triplet Average coverage depth c(v) triplet ): Where s∈S, S is all sequencing reads covering the strain-specific three-node group corresponding to this species, |s| is the number of bases in the sequencing read s covering the strain-specific three-node group corresponding to this species, and |v triplet | is the strain's specific three-node group v triplet The number of bases; Subsequently, for each strain's genomic pathway in the first strain composition, the estimated absolute abundance of the strain's genomic pathway is calculated based on the average coverage depth of each strain-specific three-node group in the strain-specific three-node group set corresponding to that genomic pathway. Specifically, this process uses the following formula to calculate the estimated absolute abundance a of the strain's genomic pathway p. p,triplet : Among them, |V strain The strain-specific three-node set V corresponds to the genomic pathway p of this strain. strain The total number of strain-specific three-node groups; Subsequently, for each strain in the first strain composition, the absolute abundance a of that strain is calculated using the following formula. p The estimated absolute abundance of strains corresponding to the strain's genomic pathway (a) p,triplet Differences between them: Finally, if the difference d is obtained strain If the value is greater than 0.45, the strain corresponding to the genomic pathway of this strain is identified as a false positive strain and filtered out; otherwise, the strain is identified as true and is included in the strain composition after the second false positive filtering.

8. The method for classifying metagenomic data at the strain level based on pangenome maps according to claim 7, characterized in that, Step (7) includes the following sub-steps: (7-1) Based on the final species composition corresponding to the metagenomic data obtained in step (4) and the absolute abundance of each species in the final species composition obtained in step (5-1), adjust the absolute abundance of each strain in the final strain composition of the metagenomic data obtained in step (6) to obtain the adjusted absolute abundance of the strain. In this step, for each species in the final species composition corresponding to the metagenomic data obtained in step (4), if the strain with the highest absolute abundance in the species exceeds the absolute abundance of the species, the absolute abundance of all strains in the species is scaled so that the sum of the absolute abundance of all strains in the species is equal to the absolute abundance of the species. (7-2) For the final strain composition of the metagenomic data, the adjusted absolute abundance of all strains in the final strain composition obtained in step (7-1) is normalized to obtain the relative abundance of each strain in the final strain composition.

Images