Species identification model, method of establishing same and use thereof

Abstract

This invention discloses a species identification model, its establishment method, and its application, belonging to the field of bioinformatics analysis technology. Based on the evolutionary relationships of species in conserved regions, a strategy for inferring the existence of species through species comparison spectra is employed. Using a known sequence distribution matrix, linear non-negative least squares is used to process the comparison results, estimating the species composition and abundance of the samples. Combining sequence detection data with a machine learning-based species composition inference model, this species identification model effectively reduces false positives in bioinformatics caused by technical noise, better reflects the actual species abundance, and thus improves the accuracy and usability of bioinformatics workflow results.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics analysis technology, and in particular to a species identification model, its establishment method, and its application. Background Technology

[0002] Metagenomic high-throughput sequencing (mNGS) is a detection technology that can directly obtain complete nucleic acid sequence information from clinical samples without relying on microbial culture. Compared with traditional culture methods, mNGS can overcome the limitations of pathogen isolation, the inability to culture some microorganisms, or excessively long culture cycles, providing a new technical approach for rapid and broad-spectrum detection of infectious diseases.

[0003] However, due to the complexity of sequencing and data analysis processes, the reliability of mNGS results is still affected by a variety of factors, making it difficult to accurately attribute sequences from different sources in the sample to the real microbial community, thus making it difficult to correctly infer its species composition.

[0004] The causes of this bias can be mainly divided into two categories: environmental noise and technical noise. Environmental noise mainly comes from environmental microorganisms, reagent contamination, or background bacteria introduced during sample preparation. Although the nucleic acids of these microorganisms are real, they are not related to actual infection. Technical noise, on the other hand, refers to the erroneous detection of species that do not physically possess nucleic acids after the analysis process.

[0005] False positives caused by technical noise are usually due to the following reasons: 1) Mutations caused by PCR amplification or sequencing errors: High-throughput sequencing may produce PCR amplification errors or sequencing errors, leading to base mutations, which may result in non-specific alignments and thus incorrect detection of species that do not exist. 2) Inaccurate or incomplete databases: If the database contains incorrect annotations or does not fully cover the genomes of certain species, it may lead to misclassification or false positives. 3) Species similarity: Species with closely related evolutionary relationships (such as microorganisms of the same family or genus) often share a large number of genomic sequences, causing specific fragments to be aligned to multiple species simultaneously.

[0006] Due to the combined influence of one or more of the above factors, a large number of false positive species are often detected after bioinformatics analysis, and even specific alignment sequences of false positive species may appear.

[0007] To improve the reliability of mNGS results, various bioinformatics false positive filtering methods have been developed, including the following categories: 1) Specific sequence filtering: This method identifies species-specific fragments (k-mers or reads) to determine the actual existence of a species, avoiding false positives caused by non-specific alignments. 2) Read count threshold filtering: This method sets a minimum read count or relative abundance threshold; species below this threshold are considered noise and excluded. 3) Intragenus relationship filtering: For species within the same genus, only the top-ranked species in terms of read count are retained, while the rest are considered background noise and removed. Alternatively, the read counts of species within the same genus are compared; if a species' read count is significantly lower than that of the top species in the same genus, it is considered a false positive.

[0008] However, some closely related species have fewer distinguishing specific sequences, which may be incorrectly filtered out; some species' specific sequences are caused by technical noise and may be incorrectly retained. Furthermore, some pathogens have low abundance in sequencing data even in the infectious state; simply filtering by abundance or read count may lead to the incorrect filtering of low-abundance pathogens. On the other hand, false positives caused by multiple alignments of high-abundance species sometimes also have high read counts, which cannot be correctly filtered out using this method. Considering the significant differences among species in actual infections, not all genera have similar first n sequences to indicate their actual presence; adopting a uniform logic will lead to the incorrect filtering or retention of some species, and the above methods cannot account for the mutual influence between species within genera. Estimating species abundance based on specific sequences leads to non-specific sequence allocation, which severely affects its accuracy.

[0009] The aforementioned issues mean that current bioinformatics analysis processes cannot accurately and reliably resolve false positive problems. Summary of the Invention

[0010] To address the issue that the aforementioned bioinformatics analysis process cannot accurately and reliably resolve false positives, this invention provides a method for establishing a species identification model. The species identification model established using this method can effectively reduce bioinformatics false positives caused by the aforementioned technical noise, better reflect the actual species abundance, and thus improve the accuracy and usability of bioinformatics process results.

[0011] On the one hand, the present invention provides a method for establishing a species identification model, comprising the following steps:

[0012] Establish a standard database: Collect genome sequences of various species to form a standard database;

[0013] Establish a reference matrix: Based on the genome sequences of the reference alignment species, simulate and generate detection reads, align the simulated reads to the standard database, and obtain the sequence number distribution of each species in the list of species to be identified, which is the alignment vector of each species. Combine these to form a reference matrix.

[0014] Sample simulation: Based on the genome sequence of each species in the list of species to be identified, test reads are generated to form a reads sampling library for that species; some species in the list of species to be identified are randomly selected, and reads are extracted from the corresponding reads sampling library to form a simulated sample. Bioinformatics analysis is performed on the simulated sample. For each detected species, if the species is the input species of the simulated sample, it is classified as positive; otherwise, it is classified as negative.

[0015] Establish a weighted sequence matrix: Obtain all reads aligned to each positive species. Based on the multiple alignment results in the standard database, calculate the weighted sequence vector for each species according to its contribution to the species. Arrange the final weighted sequence counts of all aligned species in species order to obtain the weighted sequence vector of the detected species. Then, merge the weighted sequence vectors of all species by column to obtain a weighted sequence matrix. Each column of the weighted sequence matrix is ​​the weighted sequence vector of each species, and the rows are all aligned species.

[0016] Feature acquisition: Based on the weighted sequence vectors of each species detected in the simulated samples, observation vectors are obtained. Using the observation vectors and the corresponding species matrices in the reference matrix as input, the absolute abundance vector of each species is calculated using the non-negative least squares method. And vectorize the absolute abundance of each species Normalization is performed to obtain the relative abundance vector of the samples. And obtain the relative abundance of each species s. , as feature X;

[0017] Model establishment: Using the characteristics X of each species as input and the result Y of whether the species is included in the simulated sample as output, the XGBoost model is used for training to obtain the species identification model.

[0018] Most current species abundance estimation methods are essentially based on the number of specific alignment reads or k-mer specific matches for each species to infer its relative abundance. However, some methods further optimize the accuracy by using statistical modeling or threshold filtering.

[0019] The inventors propose a strategy that infers the existence of species based on species evolutionary relationships in species-conserved regions and through species comparison spectra. Based on a known species reference sequence distribution matrix, a linear model is constructed between the sample detection sequence distribution vector and the reference matrix. This model is then processed using linear non-negative least squares to estimate the species composition and abundance of the samples. By combining sequence detection data with a species composition inference model derived from machine learning, this approach can effectively reduce false positives in bioinformatics caused by the aforementioned technical noise, better reflect the actual species abundance, and thus improve the accuracy and usability of bioinformatics workflow results.

[0020] In some implementations, the reference comparison species are selected from at least two or more species in the standard database. It is understood that the reference comparison species are selected from species that are clinically significant or commonly detected, i.e., important species requiring differentiation and identification.

[0021] In some implementations, the species in the list of species to be identified are selected from at least two or more of the reference comparison species. It is understood that the species to be identified are selected based on the species required for identification.

[0022] In some implementations, the method for simulating the generation of detection reads in the step of establishing the reference matrix is ​​as follows: The theoretical detection region of the genome is obtained, and with a mutation rate of 0.002 ± 0.002 and an insertion / deletion mutation ratio of 0.1 ± 0.1, reads positively correlated with the length of the theoretical detection region are generated using a sliding window of 40-250 bp and a step size of 1. It is understood that the sliding window design can be consistent with the expected read length.

[0023] In some implementations, the method for establishing the reference matrix in the step of establishing the reference matrix is ​​as follows: the simulated reads are aligned to the standard database to obtain the sequence number of each reference alignment species aligned to all species in the standard database, and the sequence number is normalized to obtain the reference matrix.

[0024] In some implementations, the sample simulation step involves obtaining the reads sampling library by: obtaining the genome sequence of each species in the list of species to be identified; generating reads by simulating a mutation rate of 0.002±0.002 and an insertion / deletion ratio of 0.1±0.1, with a sliding window of 40-250 bp and a step size of 1.

[0025] In some implementations, during the sample simulation step, 6-8 species are randomly selected from the list of species to be identified. For each species, a natural number is randomly drawn from 10 to 100,000, and reads of the corresponding number are drawn from the corresponding species' reads sampling library to form simulated samples with different read gradients. Preferably, the natural number is selected from 10, 100, 1000, 10000, and 100000.

[0026] In some implementations, the weighted sequence number is obtained in the step of establishing the weighted sequence matrix by the following method: counting the number of species m that a certain read can be compared with, and including the read with a weight of 1 / m in each compared species, and summing the weights of all reads that can be compared with the species that are positive, which is the weighted sequence number.

[0027] In some implementations, the absolute abundance vector is used in the feature acquisition step. Calculated using the following method:

[0028] Constructing a linear relationship:

[0029] Where b is the observation vector; S is the reference matrix; and p is the abundance vector to be estimated. Indicates unexplained residuals;

[0030] The species absolute abundance vector is solved using the nonnegative least squares method. .

[0031] In some embodiments, during the feature acquisition step, the matrix corresponding to the species in the reference matrix is ​​a submatrix containing only the corresponding species. It is understood that either the full matrix or a submatrix can be used to achieve the purpose of this invention, but using a submatrix can speed up the computation.

[0032] In some implementations, the model training parameters are set as follows during the model building step:

[0033]

[0034] In some implementations, during the model building step, the training data is divided into a training set and a validation set, with the ratio of the training set to the validation set being 7-9:1-3. The training set is used to train the model, and the validation set is used to optimize and / or evaluate the model's performance.

[0035] In some implementations, in the steps of establishing the reference matrix and sample simulation, the simulation generates detection reads that simulate detection reads using mNGS sequencing or capture metagenomics detection methods;

[0036] The feature X also includes the cosθ, cluster_size and SAR_rate parameters for each species;

[0037] The cosθ is obtained by the following method: using the reference matrix S for each species Perform fitting to obtain the fitted vector. The fitted vector is obtained according to the following formula. The cosine of the angle between the observed vector and the species, cosθ: cosθ = (b·b_fit) / (||b|| ||b_fit||), where ;

[0038] The cluster_size is the number of species contained in the cluster.

[0039] The SAR_rate is calculated using the following formula:

[0040] SAR_rate = number of species-specific sequences aligned / total number of sequences aligned to the species.

[0041] In some implementations, between the step of establishing the weighted sequence matrix and the step of obtaining the feature, there is also an inter-species clustering step, which is: calculating the unifrac distance between each pair of columns in the weighted sequence matrix, and clustering them into multiple clusters under each preset unifrac distance according to several preset distances;

[0042] In the feature acquisition step, the observation vector is obtained by extracting the columns corresponding to all species under each cluster in the weighted sequence matrix, and summing them column by column; the absolute abundance vector... The data was fitted using clusters obtained from clustering; the relative abundance of each species s was... The relative abundance proportion of this species in the resulting clusters after clustering. .

[0043] In some implementations, the preset unifrac distances for inter-species clustering are 0.05, 0.1, and 0.2. For each species *s* included in the reference matrix in a simulated sample, the preset unifrac distances of 0.05, 0.1, and 0.2 are taken respectively. cosθ, cluster_size, and SAR_rate are used as features X.

[0044] In the early stages of model building, the inventors attempted to select various different feature combinations as inputs, including common ones such as p_adjusted. However, experiments showed that using the above features as model training inputs yielded excellent results.

[0045] On the other hand, the present invention provides a method for establishing a species identification model, comprising the following steps:

[0046] Establish a standard database: Collect genome sequences of various species, which are sequences containing the 16S V1-V9 variable regions, and compile them into a standard database in taxon format;

[0047] Establish a reference matrix: Sort all sequences in the standard database and compare them from beginning to end. If the consistency of the total length of two sequences exceeds a preset range, the shorter sequence is excluded as redundancy, and the longest sequence is retained as the sequence to be included in the reference matrix. Based on the sequences included in the reference matrix, simulated reads are generated. The simulated reads are compared with usearch according to a threshold to obtain the alignment sequence vectors of each taxon in the standard database. The sequence vector of each taxon is used as a column of the matrix to form the reference matrix ref_matrix.

[0048] Establish a weighted reference matrix: For each taxon in the reference matrix, all simulated reads are weighted according to their contribution to the taxon based on their multiple alignment results in the reference database, resulting in a weighted sequence vector for each taxon. The weighted sequence vector of each taxon is used as a column of the matrix to form the weighted reference matrix; each column of the weighted sequence matrix is ​​the weighted sequence vector of each taxon, and the rows are all taxons that can be aligned.

[0049] Sample simulation: For each taxon corresponding to all sequences included in the reference matrix, simulate and generate detection reads to form a reads sampling library for that taxon. Randomly select some taxons and extract reads from the corresponding taxon's reads sampling library to form a simulated sample. Perform bioinformatics analysis on the simulated sample. For each species detected, if the species is the species corresponding to the input taxon of the simulated sample, classify it as positive; otherwise, classify it as negative.

[0050] Feature acquisition: Obtain the alignment sequence vector and weighted alignment sequence vector of the simulated samples, and fit them with the reference matrix and the weighted reference matrix respectively using the non-negative least squares method to obtain the feature of each detected taxon i. and After normalization, we obtain and Then, based on the correspondence between taxon and species, the value of each detected species s is obtained. , , and The and These are all the taxons belonging to this species. and The sum of, the and These are all the taxons belonging to this species. and The sum of the values; and calculate Wrate, Urate, and Wrate according to the following method. and :

[0051] The Wrate is the weighted sequence number of the species / the number of sequences aligned to the species;

[0052] The Urate is the number of species-specific sequences divided by the number of sequences aligned to that species.

[0053] The for ×1420 / number of aligned sequences;

[0054] The for ×1420 / number of weighted alignment sequences;

[0055] The , Wrate, Urate and Compositional feature X;

[0056] Model establishment: Using the characteristics X of each species as input and the result Y of whether the species is included in the simulated sample as output, the XGBoost model is used for training to obtain the species identification model.

[0057] In some specific embodiments of the above scheme, in the step of establishing the reference matrix, the preset range of the full-length consistency of the sequence is ≥99%.

[0058] In some specific embodiments of the above scheme, in the step of establishing the reference matrix, the method of simulating the generation of detection reads is as follows: obtain the theoretical detection region of the genome, with a mutation rate of 0.002±0.002 and an insertion and / or deletion ratio of 0.1±0.1, and with a sliding window of 40-250bp and a step size of 1, simulate the generation of reads with a number of reads positively correlated with the length of the theoretical detection region.

[0059] In some specific implementations of the above scheme, the threshold for the usearch alignment in the step of establishing the reference matrix is ​​0.97.

[0060] In some specific embodiments of the above scheme, in the sample simulation step, the reads sampling library is obtained by the following method: obtaining the genomic sequence of each taxon corresponding to all sequences included in the reference matrix, generating reads by simulating a mutation rate of 0.002±0.002, an insertion / deletion ratio of 0.1±0.1, and a sliding window of 40-250bp with a step size of 1.

[0061] In some specific embodiments of the above scheme, in the step of establishing a weighted sequence matrix, the weighted sequence number is obtained by the following method: count the number m of taxons that a certain read can be compared with, and include the read with a weight of 1 / m in each taxon that is compared with, and sum the weights of all reads in the positive taxon with a weight of 1 / m, which is the weighted sequence number.

[0062] In some specific embodiments of the above scheme, in the sample simulation step, 3-8 of the taxon corresponding to the sequence included in the reference matrix are randomly selected. For each taxon, a natural number is randomly drawn from 10-100000. Reads of the corresponding number of reads are drawn from the reads sampling library of the corresponding taxon to form simulated samples with different read gradients. Preferably, the natural number is selected from: 10, 100, 1000, 10000 and 100000.

[0063] In some specific embodiments of the above scheme, in the feature acquisition step, the absolute abundance vector Calculated using the following method:

[0064] Constructing a linear relationship:

[0065] Where b is the observation vector; S is the reference matrix; and p is the abundance vector to be estimated. Represents unexplained residuals; the observation vector includes an alignment sequence vector and a weighted alignment sequence vector, and the reference matrix includes a reference matrix and a weighted reference matrix;

[0066] The absolute abundance vector of Taxon is solved using the nonnegative least squares method. and For each taxon i, its absolute abundance is... and .

[0067] In some specific embodiments of the above scheme, the model training parameters are set as follows in the model building step:

[0068]

[0069] In some specific embodiments of the above scheme, in the model building step, the training data is divided into a training set and a validation set, with the ratio of the training set to the validation set being 7-9:1-3. The training set is used to train the model, and the validation set is used to optimize and / or evaluate the model performance.

[0070] On the other hand, the present invention also provides a species identification model established by the above-described method.

[0071] On the other hand, the present invention also provides a species identification method, comprising the following steps: taking the sequencing data of the sample to be identified for bioinformatics analysis, obtaining the feature X' of the sample to be identified according to the above-mentioned method of establishing a weighted sequence matrix and obtaining features, taking the feature X' as input, substituting it into the above-mentioned species identification model, and outputting the result Y' of the sample to be identified, thus obtaining the identification result.

[0072] On the other hand, the present invention also provides an electronic device including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the species identification method as described above.

[0073] On the other hand, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the species identification method as described above.

[0074] On the other hand, the present invention also provides a computer program product comprising a computer program, characterized in that the computer program, when executed by a processor, implements the species identification method as described above.

[0075] Based on common knowledge in the field, the above-mentioned preferred conditions can be combined arbitrarily to obtain various preferred embodiments of the present invention.

[0076] The reagents and raw materials used in this invention are all commercially available.

[0077] The positive and progressive effects of this invention are as follows:

[0078] This invention provides a method for establishing a species identification model. The species identification model established using this method can effectively reduce the false positives in bioinformatics caused by the aforementioned technical noise, better reflect the actual species abundance, and thus improve the accuracy and usability of bioinformatics process results.

[0079] Furthermore, the species identification model of this invention for sample species identification does not rely on the detection of specific alignment sequences, nor on artificially defined indicators such as sequence number threshold, sequence number multiple threshold, and species ranking. It has the advantages of small fluctuations and high stability, and can also more accurately estimate the abundance of each species in the sequencing data. Attached Figure Description

[0080] Figure 1 This is a schematic diagram of the performance of the model in Example 1.

[0081] Figure 2 This is a schematic diagram of the performance of the model in Example 2.

[0082] Figure 3 This is a schematic diagram of the performance of the model in Example 3.

[0083] Figure 4 This is a schematic diagram of the performance of the model in Example 5.

[0084] Figure 5 This is a schematic diagram of the performance of the model in Example 6.

[0085] Figure 6 This is a schematic diagram illustrating the abundance estimation performance of the model in Example 6.

[0086] Figure 7 This is a schematic diagram of the abundance estimation performance of Bracken in Example 6. Detailed Implementation

[0087] The present invention is further illustrated below by way of embodiments, but the invention is not limited to the scope of the embodiments described herein. Experimental methods in the following embodiments that do not specify specific conditions were performed according to conventional methods and conditions, or as selected according to the product instructions.

[0088] Example 1

[0089] A species identification model for identifying species in metagenomic data captured by infectious pathogen detection was established using the following method:

[0090] I. Standard Database

[0091] It includes microbial species detected in human-related samples, covering common clinical pathogens, colonizing bacteria, opportunistic pathogens, and common environmental microorganisms, totaling 25,863 species. High-quality genome sequences of the aforementioned species are collected to construct a standard database.

[0092] II. Establishing a Reference Matrix

[0093] 1563 pathogenic microorganisms that are of significant importance or commonly detected in clinical testing were selected as species to be included in the reference matrix, which are the reference comparison species. In this embodiment, species existing in CRP and species with high historical detection rates were selected, and their genome sequences were obtained.

[0094] Based on the probe set used in the metagenomic capture detection process, the theoretical capture region of the genome of the aforementioned species is obtained. Using a mutation rate of 0.002, an insertion / deletion mutation ratio of 0.1, a sliding window of 40 bp, and a step size of 1, simulated reads are generated that are positively correlated with the length of the theoretical capture region of the probe set used for detection. For a single species, n reads are simulated based on its genome, where n is a value positively correlated with the length of the theoretical capture region of that species. These are then normalized to a matrix equivalent to 10,000 sequences. For example, if the theoretical capture region length of species A is 'a' and the theoretical capture region length of species B is 'b', then the simulated read count for species A is n. A The number of simulated reads n for species B B The ratio is approximately a / b.

[0095] After standard metagenomic bioinformatics workflow analysis (including quality control, host removal, plasmid removal, and alignment with reference genome databases), the number of aligned sequences for each species to the aforementioned standard databases is generated. For example, species A simulates to obtain n. A Each read is compared to a standard database to obtain the sequence number distribution of each species in the standard database, which is the alignment vector of species A. The alignment vectors of each species are then combined to form a reference matrix ref_matrix.

[0096] Understandably, theoretically all species in the standard database could be included in the species for which the reference matrix is ​​built, but this would increase the time for generating the reference matrix and subsequent running time.

[0097] An example of how to build a reference matrix (ref_matrix) is as follows: Assume the species included in the reference matrix are Escherichia_coli, Klebsiella_pneumoniae, and Mycobacterium_tuberculosis, and the species in the standard database are Escherichia_coli, Shigella*, Klebsiella_variicola, Klebsiella_pneumoniae, and Mycobacterium_tuberculosis. For the species included in the reference matrix, Escherichia_coli (theoretical capture area approximately 0.6M), Klebsiella_pneumoniae (theoretical capture area approximately 0.5M), and Mycobacterium_tuberculosis (theoretical capture area approximately 1.7M), 60,000, 50,000, and 170,000 reads are extracted respectively, and aligned to the standard database to obtain the alignment matrix shown in the table below.

[0098] Table 1. Comparison Matrix

[0099]

[0100] Normalize each species in the table above to a matrix equivalent to 10,000 input sequences to obtain the normalized reference matrix in the table below.

[0101] Table 2. Schematic diagram of the normalized reference matrix ref_matrix

[0102]

[0103] The above process is repeated 3 times, and the final reference matrix is ​​the average of the results.

[0104] III. Clinical Sample Simulation

[0105] Forty common pathogenic microorganisms detected in clinical testing were selected (i.e., species to be identified). Reads were generated by simulating a sequencing error rate of 0.002, an insertion / deletion mutation ratio of 0.1, a sliding window of 40 bp, and a step size of 1. This constitutes the reads sampling library for the species.

[0106] The simulation method for the sample is as follows: 6-8 species are randomly selected from the above 40 species. For each species, a value is randomly selected from (10, 100, 1000, 10000, 100000) as the corresponding number of reads extracted from the above reads sampling library for that species. The reads of 6-8 species with different reads gradients form a simulated sample.

[0107] Perform metagenomic bioinformatics analysis on the simulated samples as in step one. For each species detected, if the species is the input species designed for the sample, classify it as positive; otherwise, classify it as negative.

[0108] This embodiment simulated a total of 1,000 simulated samples and 373,565 species detection records, of which 7,497 species were detected as positive and 366,068 species were detected as negative.

[0109] IV. Constructing the weighted sequence matrix of the samples

[0110] 1. Obtain the weighted sequence vector

[0111] For each species detected in the simulated sample, all reads aligned to that species in the simulated sample (including specific alignment sequences and non-specific multiple alignment sequences) are counted. For these reads, a weighted sequence vector for each species is obtained by weighting them based on their multiple alignment results in the standard database.

[0112] Specifically, the method for calculating the weighted sequence vector of a specific detected species is as follows:

[0113] 1) Collect reads: Obtain all n reads in the sample that can be matched to the detected species.

[0114] 2) Calculate the weighted contribution of a single read: For each read, count the number of species m that it can be aligned to, and record the contribution of the read to each aligned species as 1 / m. For example, when a read can be aligned to 3 species in the standard database, its weighted contribution to each species is 1 / 3.

[0115] 3) Accumulated weighted contribution: For all n reads, the weighted contribution of each aligned species is accumulated (i.e. summed) to obtain the final weighted sequence number for each aligned species.

[0116] 4) Forming a weighted sequence vector: Arrange the final weighted sequence numbers of all compared species in the specified species order to obtain the weighted sequence vector of the detected species.

[0117] 2. Obtain the weighted sequence matrix.

[0118] The weighted sequence vectors of all species obtained in the above steps are merged column-wise to obtain the weighted sequence matrix ref_matrix_w.

[0119] For example, suppose a total of 15 reads were detected in the current sample, and the alignment results are as follows (the numbers in parentheses represent the weighted contribution of each read to each species):

[0120] Read 1: E (E=1)

[0121] Read 2: E, S (E=1 / 2, S=1 / 2)

[0122] Read 3: E, S (E=1 / 2, S=1 / 2)

[0123] Read 4: E, S (E=1 / 2, S=1 / 2)

[0124] Read5: E, S, Kv, Kp (E=1 / 4, S=1 / 4, Kv=1 / 4, Kp=1 / 4)

[0125] Read 6: Kp (Kp=1)

[0126] Read 7: Kp (Kp=1)

[0127] Read8: Kp, S, E (Kp=1 / 3, S=1 / 3, E=1 / 3)

[0128] Read9: Kp, Kv (Kp=1 / 2, Kv=1 / 2)

[0129] Read10: Kp, Kv (Kp=1 / 2, Kv=1 / 2)

[0130] Read 11: M (M=1)

[0131] Read 12: M (M=1)

[0132] Read 13: M (M=1)

[0133] Read 14: M (M=1)

[0134] Read 15: M (M=1)

[0135] The weighted sequence alignment vector of Escherichia coli (E) is calculated as follows: In this example, all sequences aligned to E are Read1, Read2, Read3, Read4, Read5, and Read8. These sequences can be aligned to species E, S, Kv, and Kp. The weighted contribution of each sequence of these species is summed to obtain the weighted sequence vector of E (arranged in the order of E, S, Kv, Kp): E vector=[3.08,2.08,0.25,0.58].

[0136] The specific calculation process is as follows:

[0137] E: 1 + 1 / 2 + 1 / 2 + 1 / 2 + 1 / 4 + 1 / 3 = 3.08

[0138] S: 1 / 2 + 1 / 2 + 1 / 2 + 1 / 4 + 1 / 3 = 2.08

[0139] Kv: 1 / 4 = 0.25

[0140] Kp: 1 / 4 + 1 / 3 = 0.58

[0141] Similarly, we can obtain the weighted sequence vectors of S, Kv, Kp, and M respectively, and combine them to obtain the weighted sequence matrix, as shown below.

[0142] Table 3. Weighted sequence matrix

[0143]

[0144] In the weighted sequence matrix above, each column represents the sequence alignment for each species, showing all species whose sequences can be aligned to that species.

[0145] V. Inter-species clustering of the samples.

[0146] Using the column vectors of the matrix as feature representations of the detected species, the unifrac distance between each pair of detected species is calculated. Specifically, the weighted_normalized_fp64 function of Python's unifrac library is used to calculate the unifrac distance between each pair, and DBSCAN is used to cluster the species at distances of 0.05, 0.1, and 0.2, forming multiple clusters at each distance.

[0147] VI. Feature Acquisition

[0148] 1. Training features are obtained through non-negative least squares method.

[0149] For each cluster, extract the columns of all species under that cluster in the matrix, and sum them column by column to obtain the observation vector of that cluster.

[0150] For example, assuming that at a distance of 0.2, the above samples are clustered into cluster1: E, S, Kv, Kp; and cluster2: M. The weighted sequences of each cluster are summed to obtain the observation vector for that cluster, as shown in the table below.

[0151] Table 4. Observation vectors for each cluster (b)

[0152]

[0153] Using the observed vector (b) and the submatrix (S) of the species reference matrix that lists the compared species (i.e., the species included in the reference matrix) as input, a linear relationship is constructed:

[0154]

[0155] Where b is the observation vector; S is the reference matrix; and p is the abundance vector to be estimated. This indicates an unexplained residual.

[0156] The species absolute abundance vector is solved using the nonnegative least squares method.

[0157]

[0158] The relative abundance vector is obtained after normalization. .

[0159]

[0160] Taking the observation vector of cluster1 as an example: using the observation vector ([5.99, 4.99, 3, 5.99, 0]) and the species reference matrix as input, the abundance vector arranged in the column order of the reference matrix is ​​obtained by solving the non-negative least squares method. (In the cluster1 example, among E, S, Kv, Kp, the species included in the reference matrix are [Escherichia_coli, Klebsiella_pneumoniae]), the estimated... The values ​​are [0.0005, 0.0005]. After normalization within cluster 1, the relative proportions of both species in cluster 1 are both 0.5 (0.0005 / (0.0005+0.0005)), i.e. = 0.5, = 0.5, [0.5, 0.5].

[0161] At the same time for cluster 1 By performing a fitting operation, the fitted vector of the cluster is obtained. Let be denoted as b_fit, where S is the reference matrix ref_matrix (as shown in Table 2). For example, the above cluster1... That is, [0.0005, 0.0005], through The sequence distribution of [Escherichia_coli, Shigella*, Klebsiella_variicola, Klebsiella_pneumoniae, Mycobacterium_tuberculosis] can be obtained as a fitting result, and this vector is the fitting vector b_fit.

[0162] For each species *s*, clustering thresholds of 0.05, 0.1, and 0.2 were selected for clustering, and the relative abundance proportion of the species in the cluster was obtained when clustering was performed at different thresholds. Obtain the fitted vector for this cluster. The cosine of the angle between the cluster and the observation vector, cosθ, and the number of species in the cluster, cluster_size.

[0163] Specifically, cosθ = (b·b_fit) / (||b|| ||b_fit||) (where...) )

[0164] 2. Calculate the SAR_rate of the species.

[0165] Calculate using the following formula:

[0166] SAR_rate = Number of species-specific sequences aligned / Total number of sequences aligned to this species

[0167] The aforementioned specific sequence refers to the sequence that is uniquely matched to this species.

[0168] VII. Model Establishment

[0169] For each species *s* in the reference matrix for each simulated sample, clustering thresholds of 0.05, 0.1, and 0.2 are respectively taken. The 10 statistical measures, cosθ, cluster_size, and SAR_rate, are used as the feature matrix X, and whether the species is involved is used as the result Y. The XGBoost model of the xgboost package in Python is trained and validated according to the parameter settings in the table below.

[0170] Table 5. XGBoost Model Parameters

[0171]

[0172] The training set comprised 70% of the samples, while the validation set comprised 30%. PRAUC (Area Under the Precision-Recall Curve) and F1 score were calculated. The trained model achieved a PRAUC of 0.99 on the validation set. Figure 1 ), F1 is 0.98.

[0173] Example 2

[0174] Referring to the method in Example 1, another 30 species, different from those included in the analysis in Example 1, were selected, and 3,000 samples were simulated using the same method, resulting in a total of 472,059 species detection records.

[0175] Substituting the simulated sample data above into the model established in Example 1, and analyzing the detection records of these samples using the model trained in Example 1, the PR AUC value was 0.99 ( Figure 2 ), F1 is 0.96.

[0176] Example 3

[0177] A total of 6,535 real clinical samples were selected, containing 121,045 detection records. Among them, the reported species are positive species, and the other unreported species are excluded as potentially existing species such as environmental background bacteria, contaminating bacteria, and colonizing bacteria. Other detection records are defined as negative species.

[0178] The above data was used to build a model according to the modeling method in Example 1, with 70% as the training set and 30% as the test set. The XGBoost model from the xgboost package in Python was used for training and validation. The PR AUC value was 0.99 ( Figure 3 F1 is 0.99.

[0179] Example 4

[0180] Samples from 40 species selected in Example 1, which were clinically validated in our unit and included or excluded, were selected. Sequencing data were obtained using the same detection method as in Example 1. These data were then substituted into the model established in Example 1. The model trained in Example 1 was used to analyze these samples and output positive or negative results for the identification of specific species. All results were verified to be consistent with clinical validation.

[0181] Example 5

[0182] Following the method described in Example 1, a reference matrix was formed from the genomes of the species in the database according to their complete genome simulation sequences. Simultaneously, the simulation samples were simulated using normal mNGS sequencing data. Apart from this, all other steps were consistent with the method described in Example 1.

[0183] A total of 1294 samples were simulated, with 70% used as the training set and 30% as the test set. The PR AUC value was 0.99. Figure 4 ), F1 is 0.98.

[0184] Example 6

[0185] A species identification model for 16s species identification is established using the following method:

[0186] I. Standard Database

[0187] The greengene2 database (https: / / ftp.microbio.me / greengenes_release / 2024.09 / ) was used as the standard database.

[0188] II. Establishing a Reference Matrix

[0189] Sequences containing the 16S V1-V9 variable region from the greengene2 database were screened, all sequences were sorted, and compared from beginning to end. If the full length similarity of two sequences was ≥99%, the shorter sequence was excluded as redundant, and only the longest sequence was kept as the sequence to be included in the reference matrix.

[0190] For each sequence in the processed database, simulate 60bp reads with a step size of 1, i.e. simulate the generation of detection reads. Align these reads with a threshold of 0.97 using usearch to obtain the distribution of the number of aligned sequences to each sequence in the standard database. Combine the results of all sequences to form a reference matrix ref_matrix.

[0191] For example, assuming there are 3 taxons included in the reference matrix, and the standard database also contains 3 taxons, all simulated reads (approximately 1420 reads) of taxon1, taxon2, and taxon3 are compared to the standard database using usearch, resulting in the reference matrix shown in the table below.

[0192] Table 6. Reference Matrix (ref_matrix)

[0193]

[0194] III. Establishing a weighted reference matrix

[0195] The simulated sequence and its alignment results from step two are also used to generate the weighted reference matrix, as follows:

[0196] 1. Obtain the weighted sequence vector

[0197] For each taxon in the reference matrix, a weighted sequence vector for each taxon is obtained by weighting all its simulated sequences based on its multi-pair results in the standard database.

[0198] Specifically, the method for calculating the weighted sequence vector of a given taxon is as follows:

[0199] 1) Collect reads: Get the n reads of this taxon simulation.

[0200] 2) Calculate the weighted contribution of a single read: For each read, count the number of taxons m that it can be compared with, and record the contribution of the read to each taxon as 1 / m.

[0201] 3) Accumulate weighted contribution: For all n reads, accumulate the weighted contribution of each alignment taxon to obtain the final weighted sequence number of each alignment taxon.

[0202] 4) Form a weighted sequence vector: Arrange the final weighted sequence numbers of all aligned taxons in the specified taxon order to obtain the weighted sequence vector of the specified taxon.

[0203] 2. Obtain the weighted reference matrix ref_matrix_w.

[0204] The weighted sequence vector of each taxon obtained in the above steps is used as a column of a matrix to form a weighted reference matrix ref_matrix_w.

[0205] For example, assuming there are 3 taxons in the database, this embodiment includes all taxons in the standard database in the reference matrix. All simulated reads (approximately 1420 reads) of taxon1, taxon2, and taxon3 are aligned to the standard database using usearch. The specific alignment details of the 1420 reads of taxon1 are as follows (the value in parentheses is the weighted contribution of one read to the taxon):

[0206] Read1 - Read1000: taxon1, taxon2, taxon3 (1 / 3)

[0207] Read1001 - Read1400: taxon1, taxon2 (1 / 2)

[0208] Read1401 - Read1420: taxon1 (1)

[0209] Sum the weighted contributions of each taxon that can be matched by these reads to obtain the weighted sequence vector of taxon1 (arranged in taxon1, taxon2, taxon3): taxon1 vector=[553.33, 533.33,333.33].

[0210] The specific calculation process is as follows:

[0211] taxon1: 1 / 3 * 1000 + 1 / 2 * 400 + 1 * 20 = 553.33

[0212] taxon2: 1 / 3 * 1000 + 1 / 2 * 400 = 533.33

[0213] taxon3: 1 / 3 * 1000 = 333.33

[0214] Similarly, suppose the specific alignment of the 1420 reads of taxon2 is as follows: 1000 reads are aligned to taxon1, taxon2, and taxon3; 100 reads are aligned to taxon2 and taxon3; 300 reads are aligned to taxon1 and taxon2; and 20 reads are aligned to taxon2. Similarly, suppose the specific alignment of the 1420 reads of taxon3 is as follows: 1100 reads are aligned to taxon1, taxon2, and taxon3; 100 reads are aligned to taxon2 and taxon3; and 220 reads are aligned to taxon3. Calculate their alignment weighted sequence vectors using the same weighted calculation method as for taxon1.

[0215] The final output is a matrix, namely the weighted reference matrix ref_matrix_w, where each column is the sequence alignment for each taxon, and the rows are all the taxons that can be aligned.

[0216] Table 7. Weighted reference matrix ref_matrix_w

[0217]

[0218] IV. Sample Simulation

[0219] For each taxon corresponding to all sequences included in the reference matrix, simulate 60bp reads with a sequencing error rate of 0.001 and a step size of 1 to form a reads sampling library for that taxon.

[0220] The simulation method for the sample is as follows: 5 taxons are randomly selected, and each taxon randomly selects a value from (10, 100, 1000, 10000, 100000) as the number of corresponding reads drawn from the above reads sampling library for that taxon. The reads of 5 taxons with different read gradients form a simulated sample (5 taxons may belong to the same species, that is, the number of species in the sample is less than or equal to 5).

[0221] These samples were compared using a search with a threshold of 0.97. For each species detected, if the species had a corresponding input taxon, it was classified as positive; otherwise, it was classified as negative.

[0222] That is, the judgment is made based on the correspondence between taxon and species. For example, taxon1 and taxon2 both belong to the species Escherichia coli, and taxon3 belongs to the species Klebsiella pneumoniae. If Escherichia coli is detected, the classification is positive if taxon1 or taxon2 is introduced, and negative if only taxon3 is introduced.

[0223] The absolute abundance of a species is the sum of the number of reads of all taxons of that species in the sample. This example contains a total of 1535 simulated samples and 556016 species detection records, of which 7623 are positive and 548393 are negative.

[0224] V. Feature Acquisition

[0225] 1. Obtaining the alignment sequence vector and weighted alignment sequence vector of the samples.

[0226] The alignment sequence vector is the number of alignment sequences that the sample aligns to all taxons, where the order of taxons is the same as the column names of the reference matrix ref_matrix.

[0227] The weighted alignment sequence vector represents the number of weighted alignment sequences of this sample to all taxons, where the taxons are arranged in the same order as the column names of the reference matrix ref_matrix. Specifically, its calculation method is as follows:

[0228] 1) Collect reads: Get all n reads for this sample.

[0229] 2) Calculate the weighted contribution of a single read: For each read, count the number of taxons m that it can be compared with, and record the contribution of the read to each taxon as 1 / m.

[0230] 3) Accumulate weighted contribution: For all n reads, accumulate the weighted contribution of each alignment taxon to obtain the final weighted sequence number of each alignment taxon.

[0231] 4) Forming a weighted alignment sequence vector: Arrange the final weighted sequence numbers of all alignment taxons in the specified taxon order to obtain the weighted alignment sequence vector of the sample.

[0232] For example, assuming 4260 reads are detected in this sample, the specific comparison details are as follows (the value in parentheses represents the weighted contribution of a single read to taxon):

[0233] Read1 - Read3100 (3100 items): taxon1, taxon2, taxon3 (1 / 3)

[0234] Read 3101 - Read 3800 (700 entries): taxon1, taxon2 (1 / 2)

[0235] Read3801 - Read4000 (200 entries): taxon2, taxon3 (1 / 2)

[0236] Read4001 - Read4020 (20 items): taxon1 (1)

[0237] Read4021 - Read4040 (20 items): taxon2 (1)

[0238] Read4041 - Read4260 (220 entries): taxon3 (1)

[0239] The alignment sequence vector (according to taxon1, taxon2, taxon3) is [3820, 4020, 3520]. The specific calculation process is as follows:

[0240] taxon1: 3100 + 700 + 20 = 3820

[0241] taxon2: 3100 + 700 + 200 + 20 = 4020

[0242] taxon3: 3100 + 200 + 220 = 3520

[0243] Its weighted alignment sequence vector (according to taxon1, taxon2, taxon3) is [1403, 1503, 1353]. The specific calculation process is as follows:

[0244] taxon1: 3100 * 1 / 3 + 700 * 1 / 2 + 20 * 1 = 1403

[0245] taxon2: 3100 * 1 / 3 + 700 * 1 / 2 + 200 * 1 / 2 + 20 * 1 = 1503

[0246] taxon3: 3100 * 1 / 3 + 200 * 1 / 2 + 220 * 1 = 1353

[0247] 2. and Obtain.

[0248] 1) Obtain the taxon level and

[0249] Using the observation vector (b, the alignment sequence vector or weighted alignment sequence vector mentioned above) and the reference matrix (S, the reference matrix or weighted reference matrix mentioned above) as input, a linear relationship is constructed:

[0250]

[0251] Where b is the observation vector; S is the reference matrix; and p is the abundance vector to be estimated. Indicates unexplained residuals

[0252] The absolute abundance vector of the Taxon level is solved using the nonnegative least squares method.

[0253]

[0254] The relative abundance vector is obtained after normalization. .

[0255]

[0256] Using taxon as the unit, the alignment vector (b) of each taxon aligned to the alignment sequence distribution in the standard database and the reference matrix ref_matrix (S) are taken as input. Through the above steps, the taxon hierarchy is obtained, arranged in the column name order of the reference matrix (in the example above, [taxon1, taxon2, taxon3]). (recorded as) )and (recorded as) For taxon i, its absolute abundance and relative abundance are respectively and .

[0257] Using taxon as the unit, the weighted sequence number distribution vector (b) of each taxon and the weighted reference matrix ref_matrix_w (S) are taken as input, and fitted using the same non-negative least squares method to obtain the taxon hierarchy arranged in the column name order of the reference matrix (in the example above, [taxon1, taxon2, taxon3]). (recorded as) )and (recorded as) For taxon i, its absolute abundance and relative abundance are respectively and .

[0258] 2) Obtain species-level information and

[0259] For each detected species *s* (derived from the corresponding taxon), calculate its species-level characteristics. and These are the corresponding taxons belonging to that species. and The sum of .

[0260] For example, suppose that after calculation using the above steps, a certain sample yields... The vector is [1,1,1] (in the order [taxon1,taxon2,taxon3]). Since taxon1 and taxon2 belong to the species *Escherichia_coli*, and taxon3 belongs to the species *Klebsiella_pneumoniae*, then the species *Escherichia_coli*... For 2, Klebsiella_pneumoniae The value is 1. The calculation is similar.

[0261] 3) Obtain species-level information and

[0262] For each detected species *s* (derived from the corresponding taxon), calculate its species-level characteristics. and These are the corresponding taxons belonging to that species. and The sum of .

[0263] For example, suppose that after calculation using the above steps, a certain sample yields... The vector is [1 / 3, 1 / 3, 1 / 3] (in the order [taxon1, taxon2, taxon3]). Since taxon1 and taxon2 belong to the species *Escherichia_coli*, and taxon3 belongs to the species *Klebsiella_pneumoniae*, then the species *Escherichia_coli*... The value is 0.67 (2 / 3) for Klebsiella pneumoniae. It is 0.33 (1 / 3). The calculation is similar.

[0264] 3. Other characteristics.

[0265] wrate: weighted sequence number of the species / number of sequences aligned to the species;

[0266] urate: the number of species-specific sequences / the number of sequences aligned to this species;

[0267] :for × 1420 / number of aligned sequences;

[0268] :for × 1420 / Number of weighted alignment sequences.

[0269] Note: Since the sequence length of 16S rRNA is relatively stable, the average number of reads in the reference matrix for each species is 1420, and this value is selected as the basis for calculation.

[0270] VI. Model Establishment

[0271] For each species s in the reference matrix for each sample, take its... , , wrate,urate, , These six statistics are used as the feature matrix X, and whether the species has an input sequence is used as the result Y. Referring to the modeling method in Example 1, the XGBoost model of the xgboost package in Python is used for training and validation.

[0272] Of these, 70% of the samples were used in the training phase and 30% in the validation phase. The F1 score was calculated, and the trained model achieved an F1 score of 0.94 and a PR AUC of 0.97 on the validation set. Figure 5 ).

[0273] VII. Species Composition and Abundance Estimation

[0274] Because the sequence length of 16S rRNA is relatively stable, the average number of reads in the reference matrix for each species is around 1420. × 1420 is the estimated number of reads for this species.

[0275] The results are as follows Figure 6 As shown, Figure 6 The abundance estimation performance of the model is shown in the figure. From left to right and top to bottom, the values ​​are L1 (Absolute Error Sum), Rmse (Root Mean Squared Error), Spearman_corr (Spearman Rank Correlation), and Pearson_corr (Pearson Correlation Coefficient). The figure shows that the model has low error, high ranking consistency, and high linear consistency.

[0276] Simultaneously, these simulated samples were analyzed using Bracken's best practice workflow (running parameter: -t 0), and the comparison database was the same as in this embodiment. The results are as follows: Figure 7 , Figure 7 The abundance estimation performance of Bracken is shown in the figure. From left to right and top to bottom, the values ​​are L1, Rmse, Spearman_corr, and Pearson_corr. It can be seen from the figure that the error is larger and the ordering consistency and linear consistency are poor.

[0277] The above results indicate that the estimation of species composition abundance in this embodiment is superior to Bracken in all indicators.

Claims

1. A method for establishing a species identification model, characterized by, Includes the following steps: Establish a standard database: Collect genome sequences of various species to form a standard database; Establish a reference matrix: Based on the genome sequences of the reference alignment species, simulate and generate detection reads, align the simulated reads to the standard database, and obtain the sequence number distribution of each species in the list of species to be identified, which is the alignment vector of each species. Combine these to form a reference matrix. Sample simulation: Based on the genome sequence of each species in the list of species to be identified, test reads are generated to form a reads sampling library for that species; some species in the list of species to be identified are randomly selected, and reads are extracted from the corresponding reads sampling library to form a simulated sample. Bioinformatics analysis is performed on the simulated sample. For each detected species, if the species is the input species of the simulated sample, it is classified as positive; otherwise, it is classified as negative. Establish a weighted sequence matrix: Obtain all reads aligned to each positive species. Based on the multiple alignment results in the standard database, calculate the weighted sequence vector for each species according to its contribution to the species. Arrange the final weighted sequence counts of all aligned species in species order to obtain the weighted sequence vector of the detected species. Then, merge the weighted sequence vectors of all species by column to obtain a weighted sequence matrix. Each column of the weighted sequence matrix is ​​the weighted sequence vector of each species, and the rows are all aligned species. Feature acquisition: based on the weighted sequence vector of each species detected in the simulation sample, an observation vector is obtained, and the observation vector and the matrix of the corresponding species in the reference matrix are taken as inputs to calculate the absolute abundance vector of each species by non-negative least squares method , and the absolute abundance vector of each species is normalized to obtain the sample relative abundance vector , and the relative abundance of each species s is obtained as the feature X; Model establishment: Using the characteristics X of each species as input and the result Y of whether the species is included in the simulated sample as output, the XGBoost model is used for training to obtain the species identification model.

2. The method for establishing a species identification model according to claim 1, wherein Meets one or more of the following conditions: (1) The reference comparison species are selected from at least two or more species in the standard database; (2) The species in the list of species to be identified are selected from at least two or more of the reference comparison species; (3) In the step of establishing the reference matrix, the method of simulating the generation of detection reads is as follows: obtain the theoretical detection region of the genome, according to the mutation rate of 0.002±0.002, the proportion of insertion and / or deletion in mutation is 0.1±0.1, and according to the sliding window of 40-250bp and the step size of 1, simulate the generation of reads that are positively correlated with the length of the theoretical detection region; (4) In the step of establishing the reference matrix, the method of establishing the reference matrix is ​​as follows: the simulated reads are aligned to the standard database to obtain the sequence number of each reference alignment species aligned to all species in the standard database, and the sequence number is normalized to obtain the reference matrix; (5) In the sample simulation step, the reads sampling library is obtained by the following method: obtain the genome sequence of each species in the list of species to be identified, generate reads by simulating the mutation rate of 0.002±0.002 and the mutation ratio of insertion / deletion of 0.1±0.1, and with a sliding window of 40-250bp and a step size of 1. (6) In the sample simulation step, 6-8 species are randomly selected from the list of species to be identified. For each species, a natural number is randomly drawn from 10 to 100,000. Reads of the corresponding number of reads are drawn from the reads sampling library of the corresponding species to form simulated samples with different read gradients. Preferably, the natural number is selected from 10, 100, 1000, 10000 and 100000. (7) In the step of establishing the weighted sequence matrix, the weighted sequence number is obtained by the following method: count the number of species m that a certain read can be compared with, and include the read with a weight of 1 / m in each species that is compared with, and sum the weights of all reads that can be compared with the species that are detected as positive, which is the weighted sequence number; (8) In the feature acquisition step, the absolute abundance vector By the following method: Constructing linear relationships: , where b is the observation vector; S is the reference matrix; p is the abundance vector to be estimated; represents the unexplained residual error; The species absolute abundance vector is solved using the nonnegative least squares method. ; (9) In the feature acquisition step, the matrix corresponding to the species in the reference matrix is ​​a submatrix containing only the corresponding species; (10) In the model building step, the model training settings are as follows: ； (11) In the model building step, the training data is divided into a training set and a validation set, and the ratio of the training set to the validation set is 7-9:1-3. The training set is used to train the model, and the validation set is used to optimize and / or evaluate the model performance.

3. The method for establishing a species identification model as described in claim 1 or 2, characterized in that, In the steps of establishing the reference matrix and sample simulation, the simulated reads generated are reads that are simulated using mNGS sequencing or capture metagenomics detection methods. The feature X also includes the cosθ, cluster_size and SAR_rate parameters for each species; The cosθ is obtained by the following method: using the reference matrix S for each species Perform fitting to obtain the fitted vector. The fitted vector is obtained according to the following formula. The cosine of the angle between the observed vector and the species' vector, cosθ: cosθ = (b·b_fit) / (||b|| ||b_fit||), where ; The cluster_size is the number of species contained in the cluster. The SAR_rate is calculated using the following formula: SAR_rate = number of species-specific sequences aligned / total number of sequences aligned to the species.

4. The method for establishing a species identification model as described in claim 3, characterized in that, Between the step of establishing the weighted sequence matrix and the step of obtaining the feature, there is also an inter-species clustering step, which is: calculating the unifrac distance between each pair of columns in the weighted sequence matrix, and clustering them into multiple clusters under each preset unifrac distance according to several preset distances; In the feature acquisition step, the observation vector is obtained by extracting the columns corresponding to all species under each cluster item in the weighted sequence matrix, and summing them column by column; the absolute abundance vector... The data was fitted using clusters obtained from clustering; the relative abundance of each species s was... The relative abundance proportion of this species in the resulting clusters after clustering. ; Preferably, the preset unifrac distances for inter-species clustering are 0.05, 0.1, and 0.

2. For each species *s* included in the reference matrix in a simulated sample, the preset unifrac distances of 0.05, 0.1, and 0.2 are respectively taken. cosθ, cluster_size, and SAR_rate are used as features X.

5. A method for establishing a species identification model, characterized in that, Includes the following steps: Establish a standard database: Collect genome sequences of various species, which are sequences containing the 16S V1-V9 variable regions, and compile them into a standard database in taxon format; Establish a reference matrix: Sort all sequences in the standard database and compare them from beginning to end. If the consistency of the total length of two sequences exceeds a preset range, the shorter sequence is excluded as redundancy, and the longest sequence is retained as the sequence to be included in the reference matrix. Based on the sequences included in the reference matrix, simulated reads are generated. The simulated reads are compared with usearch according to a threshold to obtain the alignment sequence vectors of each taxon in the standard database. The sequence vector of each taxon is used as a column of the matrix to form the reference matrix ref_matrix. Establish a weighted reference matrix: For each taxon in the reference matrix, all simulated reads are weighted according to their contribution to the taxon based on their multiple alignment results in the reference database, resulting in a weighted sequence vector for each taxon. The weighted sequence vector of each taxon is used as a column of the matrix to form the weighted reference matrix; each column of the weighted sequence matrix is ​​the weighted sequence vector of each taxon, and the rows are all taxons that can be aligned. Sample simulation: For each taxon corresponding to all sequences included in the reference matrix, simulate and generate detection reads to form a reads sampling library for that taxon. Randomly select some taxons and extract reads from the corresponding taxon's reads sampling library to form a simulated sample. Perform bioinformatics analysis on the simulated sample. For each species detected, if the species is the species corresponding to the input taxon of the simulated sample, classify it as positive; otherwise, classify it as negative. Feature acquisition: Obtain the alignment sequence vector and weighted alignment sequence vector of the simulated samples, and fit them with the reference matrix and the weighted reference matrix respectively using the non-negative least squares method to obtain the feature of each detected taxon i. and After normalization, we obtain and Then, based on the correspondence between taxon and species, the value of each detected species s is obtained. , , and The and These are all the taxons belonging to this species. and The sum of, the and These are all the taxons belonging to this species. and The sum of the values; and calculate Wrate, Urate, and Wrate according to the following method. and : The Wrate is the weighted sequence number of the species / the number of sequences aligned to the species; The Urate is the number of species-specific sequences divided by the number of sequences aligned to that species. The for ×1420 / number of aligned sequences; The for ×1420 / number of weighted alignment sequences; The , Wrate, Urate and Compositional feature X; Model establishment: Using the characteristics X of each species as input and the result Y of whether the species is included in the simulated sample as output, the XGBoost model is used for training to obtain the species identification model.

6. The method for establishing a species identification model as described in claim 5, characterized in that, Meets one or more of the following conditions: (1) In the step of establishing the reference matrix, the preset range of the full-length consistency of the sequence is ≥99%; (2) In the step of establishing the reference matrix, the method of simulating the generation of detection reads is as follows: obtain the theoretical detection region of the genome, according to the mutation rate of 0.002±0.002, the proportion of insertion and / or deletion in mutation is 0.1±0.1, and according to the sliding window of 40-250bp and the step size of 1, simulate the generation of reads that are positively correlated with the length of the theoretical detection region; (3) In the step of establishing the reference matrix, the threshold for the usearch comparison is 0.97; (4) In the sample simulation step, the reads sampling library is obtained by the following method: obtain the genome sequence of each taxon corresponding to all sequences included in the reference matrix, generate reads by simulating the mutation rate of 0.002±0.002 and the mutation ratio of insertion / deletion of 0.1±0.1, and with a sliding window of 40-250bp and a step size of 1. (5) In the step of establishing the weighted sequence matrix, the weighted sequence number is obtained by the following method: count the number of taxons m that a certain read can be compared with, and include the read with a weight of 1 / m into each taxon that is compared with, and sum the weights of all reads in the positive taxon with a weight of 1 / m, which is the weighted sequence number. (6) In the sample simulation step, 3-8 of the taxon corresponding to the sequence included in the reference matrix are randomly selected. For each taxon, a natural number is randomly drawn from 10-100000. Reads with the corresponding number of reads are drawn from the reads sampling library of the corresponding taxon to form simulated samples with different read gradients. Preferably, the natural number is selected from: 10, 100, 1000, 10000 and 100000. (7) In the feature acquisition step, the absolute abundance vector Calculated using the following method: Constructing a linear relationship: , Where b is the observation vector; S is the reference matrix; and p is the abundance vector to be estimated. Represents unexplained residuals; the observation vector includes an alignment sequence vector and a weighted alignment sequence vector, and the reference matrix includes a reference matrix and a weighted reference matrix; The absolute abundance vector of Taxon is solved using the nonnegative least squares method. and For each taxon i, its absolute abundance is... and ; (8) In the model building step, the model training settings are as follows: ； (9) In the model building step, the training data is divided into a training set and a validation set, and the ratio of the training set to the validation set is 7-9:1-3. The training set is used to train the model, and the validation set is used to optimize and / or evaluate the model performance.

7. The species identification model established by the method according to any one of claims 1-5.

8. A method for species identification, characterized in that, Includes the following steps: Sequencing data of the sample to be identified is used for bioinformatics analysis. The feature X' of the sample to be identified is obtained according to the feature acquisition step of any one of claims 1-6. The feature X' is used as input and substituted into the species identification model of claim 7. The result Y' of the sample to be identified is output, which is the identification result.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the species identification method as described in claim 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the species identification method as described in claim 8.

Images