A bioinformatics method for metagenomic detection of pathogenic microorganisms based on next-generation sequencing

Abstract

The present invention discloses a bioinformatics method for pathogenic microorganism metagenomic detection based on second-generation sequencing technology, covering the steps of creating a species list, downloading genome sequences, building a database, quality control of sequencing data, removing human sequences, obtaining a strain identification list A based on k-mer and LCA algorithms, obtaining a strain identification list B based on bwa-mem algorithms, using list B to correct list A, and strain annotation and screening. Compared with the solutions currently on the market, the method provided by the present invention is significantly different. First, in order to ensure the accuracy and comprehensiveness of the analysis results, a human sequence database, a microbial comparison database, and a microbial annotation database were established. In addition, the present invention combines the current mainstream strain identification method based on k-mer and bwa-mem comparisons, realizes the organic combination of the two algorithms, and shows obvious advantages in the identification of low-abundance species. These improvements bring more reliable and more effective solutions to the detection of pathogenic microorganism metagenomics.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and in particular to a bioinformatics method for detecting pathogenic microorganism metagenomes based on second-generation sequencing. Background Art

[0002] Metagenomics is a method to study the genetic information of all microorganisms in a microbial community. Unlike traditional genomics research, metagenomics aims to analyze the genome of the entire microbial community rather than the genome of a single microorganism. It extracts DNA or RNA directly from environmental samples, performs high-throughput sequencing of all microbial DNA or RNA in them, and then uses bioinformatics technology to analyze and interpret these mixed DNA or RNA sequence data; Illumina second-generation sequencing technology is a large-scale parallel sequencing technology with the characteristics of high throughput, high sequencing accuracy, and short time efficiency. It is the most commonly used high-throughput sequencing technology. Pathogenic microorganisms refer to microorganisms that can cause disease or infection. These microorganisms include bacteria, viruses, fungi, and parasites. They can enter the human body, animals, or plants, destroy the normal physiological functions of the host, and cause the occurrence of diseases. Metagenomic sequencing (mNGS) can identify all microorganisms in the sample and identify new species, which is impossible to achieve using conventional experimental methods such as smears, immunological detection of antigens and antibodies, culture identification, mass spectrometry detection based on metabolites, and detection technology based on multiple PCR. mNGS is widely used in the detection of clinical infectious pathogens such as suspected clinical pathogens, rare pathogens, and critical and severe diseases due to its advantages such as no need for culture, independence of probe sequences, and wide coverage without preference.

[0003] At present, the commonly used bioinformatics methods for species identification in the field of mNGS include Kraken2 and BWA. Kraken2 is a k-mer-based rapid species identification tool that identifies the microbial composition in a sample by matching the sequence with the k-mer in the known genome database. In contrast, BWA (Burrows-WheelerAligner) is a commonly used sequence alignment tool, which is usually used to align sequencing data with the reference genome to determine the source of the sequence. However, both methods have obvious disadvantages. The identification results of Kraken2 depend on the database built, which may be limited by the coverage and accuracy of the database, and the accuracy of the identification of homologous or similar species has always been insufficient; BWA is a fast and efficient sequence alignment tool, but BWA itself does not provide species identification function, so in the process of species identification, additional analysis steps and tools are often required to interpret the alignment results and perform species identification. Based on this, developing a method that combines Kraken2 and BWA is of great significance to improve the analysis efficiency and accuracy of mNGS data. Summary of the invention

[0004] The present invention provides a bioinformatics method for pathogenic microorganism metagenomic detection based on second-generation sequencing, which specifically comprises the following steps:

[0005] 1. Create a species list. The present invention uses the Taxonomy classification system in the NCBI database to create a species list including bacteria, fungi, viruses and parasites.

[0006] 2. Download of genome sequences. The data that needs to be downloaded in the present invention include the genome sequences of all species in the list. The genome sequences used to construct the database come from various public databases, including Refseq, GenBank, FDA-ARGOS and PATRIC for bacteria; FungiDB, Refseq, GenBank for fungi; RVDB, Refseq, GenBank, GISAID for viruses; and Refseq, GenBank, and WormBaseParaSite for parasites.

[0007] 3. Database construction. The databases constructed by the present invention include a microbial comparison database, a microbial reference genome database, a microbial annotation database, and a human sequence database. The present invention uses kraken2 software to construct a microbial comparison database; the genome sequences marked as reference in the NCBI database in the list are summarized together to form a microbial reference genome database; the microbial annotation database includes the genus, biological classification, Latin name, pathogenicity, clinical attention, and species epidemiological information of each microorganism in the list; the human sequence database includes the hg19 version of the human reference genome and transcriptome, the Yanhuang No. 1 genome, and T2T-CHM13.

[0008] 4. Quality control of the next-generation sequencing data of pathogen metagenomic (mNGS) clinical samples and negative control samples, including removal of adapter sequences, removal of low-quality bases at the beginning and end, removal of sequences containing N bases in reads, removal of sequences with a single base accounting for more than 80% or two bases accounting for more than 90%, and removal of repeated sequences and low-complexity sequences.

[0009] 5. Remove human sequences, compare the quality-controlled data with the human sequence database, remove the matched sequences, and the remaining ones are high-quality non-human sequences.

[0010] 6. According to the microbial comparison database, kraken2 software was used, based on k-mer and LCA algorithms and combined with self-written Python scripts to obtain the bacterial species identification list A. List A includes indicators such as the number of bacterial species reads, the number of unique reads of bacterial species, the relative abundance of bacterial species, the proportion of mer aligned to the species level, the proportion of unique reads, and the foldchange value; According to the microbial reference genome database, based on the bwa-mem algorithm and combined with self-written Python scripts, the bacterial species identification list B was obtained. List B includes indicators such as the number of aligned reads, the number of aligned unique reads, the average MAPQ value of the aligned reads, the average match rate of the aligned reads, the genome coverage, and the foldchange value.

[0011] 7. Use strain list B to correct list A to get the final strain list. The correction principle is as follows:

[0012] (1) Extract the common bacterial species in List A and List B to form Set S1. The number of reads detected for the common strains is based on the number of reads of the species with more reads in the two methods. The bacterial species in List A but not in List B form Set S2. The bacterial species not in List A but in List B form Set S3.

[0013] (2) Take out the bacterial species in list A in set S1, and calculate the score of each bacterial species according to the five indicators in the A-S1 indicator set: relative abundance of bacterial species, proportion of mer aligned to species level, proportion of unique reads, foldchange value, and proportion of bacterial species reads in all reads aligned to microorganisms. And assign a corresponding coefficient to each indicator. When the foldchange value is positive infinity, the foldchange value is assigned to 1. When the foldchange value is greater than or equal to 5 and not positive infinity, it is assigned to 0.5. When the foldchange value is less than 5, it is assigned to 0.2. Calculate the ScoreA of each bacterial species according to the following formula:

[0014]

[0015] a i is the corresponding coefficient of each indicator in the A-S1 indicator set, X i is the score of each indicator of the A-S1 indicator set, and m is the number of indicators. The A-S1 indicator set includes species relative abundance, species-level mer ratio, unique read ratio, foldchange value, and the ratio of species reads to all reads aligned to microorganisms;

[0016] (3) Take the bacterial species in list B in set S1, and calculate the score of each bacterial species based on the six indicators in the B-S1 indicator set: the number of aligned reads, the number of aligned unique reads, the average MAPQ value of aligned reads, the average match rate of aligned reads, the genome coverage, and the foldchange value; similarly, when the foldchange value is positive infinity, the foldchange value is assigned to 1, when the foldchange value is greater than or equal to 5 and not positive infinity, it is assigned to 0.5, and when the foldchange value is less than 5, it is assigned to 0.2; assign a corresponding coefficient to each indicator, and calculate the ScoreB of each bacterial species according to the following formula:

[0017]

[0018] Where bi is the corresponding coefficient of each indicator in the B-S1 indicator set, Yi is the score of each indicator in the B-S1 indicator set, and n is the number of indicators. The B-S1 indicator set includes the number of aligned reads, the number of aligned unique reads, the average MAPQ value of aligned reads, the average match rate of aligned reads, the genome coverage, and the foldchange value;

[0019] (4) Sort the ScoreA values ​​of all bacterial species in List A in Set S1 from small to large, take the first quartile as Threshold 1, calculate ScoraA of all bacterial species in S2, and those with ScoreA values ​​above Threshold 1 are credible results; Sort the ScoreB values ​​of all bacterial species in List B in Set S1 from small to large, take the first quartile as Threshold 2; calculate ScoreB of all bacterial species in S3, and those with ScoreB values ​​above Threshold 2 are credible results;

[0020] (5) Summarize the reliable results in sets S2 and S3 with those in set S1 to obtain the final strain identification list.

[0021] 8. Based on the self-written Python script, use the microbial annotation database to annotate each microorganism in the strain identification list. The annotation information includes each microorganism's genus, biological classification, pathogenicity, clinical concern, and species epidemiological information.

[0022] 9. Screen the final strain identification list after annotation, and screen out the final potential pathogenic microorganisms based on the specific values ​​of each indicator in the strain list, microbial pathogenicity, clinical complaints, patient symptoms and other information.

[0023] Compared with the prior art, the present invention has the following beneficial effects:

[0024] (1) The present invention is based on metagenomic second-generation sequencing technology, which can simultaneously detect all microorganisms in the sample, including bacteria, viruses, fungi, etc., and can discover microorganisms that may be missed by traditional methods, thereby improving the comprehensiveness and sensitivity of the detection; at the same time, there is no need to culture microorganisms, and the detection and analysis of samples can be completed in a shorter time, greatly shortening the detection cycle.

[0025] (2) The bioinformatics method provided by the present invention combines the advantages of Kraken2 and BWA, solving the problem that Kraken2 has low classification accuracy for similar or homologous species and BWA can only compare but not identify species. At the same time, it improves the detection accuracy of low-abundance species, realizes rapid and accurate detection of pathogenic microorganisms, and provides important technical support for disease diagnosis and epidemiological research. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 Flowchart of a bioinformatics method for metagenomic detection of pathogenic microorganisms based on next-generation sequencing.

[0027] Figure 2 Flowchart for determining threshold 1 according to set S1A.

[0028] Figure 3 Flow chart for determining threshold 2 according to set S1B.

[0029] Figure 4 Box plot of average accuracy of identification of Escherichia coli of 2000 simulated samples by the present invention and Kraken2.

[0030] Figure 5 Box plot of the average accuracy of identification of Shigella flexneri of 2000 simulated samples by the present invention and Kraken2.

[0031] Figure 6 Box plot of the average number of detected reads for identification of Escherichia coli in 2000 simulated samples by the present invention and kraken2.

[0032] Figure 7 Box plot of the average number of detected reads for identification of Shigella flexneri of 2000 simulated samples by the present invention and kraken2. DETAILED DESCRIPTION

[0033] In order to more clearly illustrate the contents disclosed by the present invention, including the purpose, technical solutions and advantages, we will use drawings and detailed descriptions to illustrate. The purpose of providing these embodiments is to make the understanding of the present invention more thorough and comprehensive.

[0034] Unless otherwise defined, all technical and scientific terms used herein are generally understood by those skilled in the art in the art of the present invention. The terms used in the specification of the present invention herein are only for the purpose of describing specific embodiments and are not intended to limit the present invention.

[0035] Example 1

[0036] like Figure 1 As shown, a bioinformatics method for pathogenic microorganism metagenomic detection based on second-generation sequencing provided by a specific embodiment of the present invention comprises the following steps:

[0037] 1. Create a species list, including bacteria, fungi, viruses and parasites. According to the classification of Taxonomy in the NCBI database, sort it down from kingdom, phylum, class, order, family, genus and species step by step, and finally organize it into a complete species list. Each line of the list represents a genome sequence information of the species, including genome assembly_accession number, species Latin name, taxid, species_taxid, assembly level, sequence release time, species classification, genome size and database source, among which taxid is the genome sequence number annotated in the NCBI database, and species_taxid is the species number.

[0038] 2. Obtain the genome sequence in the species list, enter through the link provided in the table below, find the corresponding genome and download it.

[0039] Table 1. Databases used in the present invention and their download links

[0040] database Download Link RefSeq https: / / ftp.ncbi.nlm.nih.gov / genomes / refseq / GenBank https: / / ftp.ncbi.nlm.nih.gov / genomes / genbank / PATRIC https: / / www.patricbrc.org / FDA-ARGOS https: / / www.ncbi.nlm.nih.gov / bioproject / 231221 FungiDB https: / / fungidb.org / fungidb / app / downloads RVDB https: / / rvdb.dbi.udel.edu / GISAID https: / / gisaid.org / WormBaseParaSit https: / / parasite.wormbase.org / ftp.html

[0041] 3. Database construction. The databases constructed by the present invention include a microbial comparison database, a microbial reference genome database, a microbial annotation database, and a human sequence database. The genome sequences in the list are subjected to a series of filtering operations, including the removal of plasmids, fuzzy names, sequences less than 1000bp in length, and low repetitive sequences; the comparison database is constructed using the kraken2 software. First, the command kraken2-build--download-taxonomy--threads 96--dbkraken2_custom_db is used to download the classification information of the NCBI database Taxonomy, and the classification information is stored in the kraken2_custom_db folder; then the genome sequence is added to the database through find genome-name'*.fna'-print0|xargs-0-I{}-n1kraken2-build--add-to-library{}--db kraken2_custom_db; finally, kraken2-build--build--threads 96--db The kraken2_custom_db command constructs a microbial comparison database, and kraken2_custom_db is the name of the created microbial comparison library; the microbial reference genome database is composed of the genome sequences of the species in the list that are marked as reference in the NCBI database; the microbial annotation database includes the genus, biological classification, Latin name, pathogenicity, clinical attention, and species epidemiological information of each microorganism in the list; in addition, the present invention constructs a comprehensive human sequence database that is more suitable for the Chinese race. In addition to the conventional hg19 version of the human reference genome and transcriptome, it also adds the Yanhuang No. 1 genome sequence specifically for the Chinese race and the gapless version of the human genome sequence T2T-CHM13.

[0042] 4. In the present invention, fastp v0.23.2 is used to perform quality control on the original data of pathogen metagenomic (mNGS) sequencing, and the parameter is set to --dedup--length_required 45-q 20-n 1. First, the repeated sequences are removed, and then the reads containing N bases are filtered out, and the read fragments with a quality value lower than 20 are cut off. In addition, the adapter sequence and the low-complexity sequence are cut off, and the sequence with a read length of less than 45bp after quality control is filtered out; in addition, for the entire sample, the sample is qualified only if the three conditions of Q30 base data volume greater than 80%, adapter contamination ratio less than 1% and the effective data volume after quality control must be greater than 70% are met at the same time.

[0043] 5. In the present invention, bowtie v2-2.4.4 is used to align the data after quality control to the human sequence database. The human sequence database consists of the hg19 version of the reference genome, transcriptome data (refMrna), the world's first Chinese standard genome sequence Yanhuang No. 1, and the gapless version of the human genome sequence T2T-CHM13 published in March 2022. The reads that cannot be aligned to the human sequence database are saved through the --un parameter to obtain high-quality non-human sequences.

[0044] 6. According to the microbial comparison database kraken2_custom_db, kraken2 software was used, and the command kraken2--db kraken2_custom_db sample.remove-hg19-chm13-yh-tr.fq--threads 10--report sample_kraken2_h.report--output sample_kraken2_h.output--use-names was used in combination with the self-written Python script to obtain the bacterial species identification list A. The following table is an example of list A (part of the rows and columns are intercepted):

[0045] Table 2. Example of extracting some rows and columns from List A

[0046]

[0047]

[0048] List A includes the number of strain reads, the number of strain unique reads, the relative abundance of strains, the proportion of mer aligned to the species level, the foldchange value and other indicators; according to the microbial reference genome database, based on the bwa-mem algorithm, the bwa operation command: bwa mem-t 10{bwaindex}sample.remove-hg19-chm13-yh-tr.fq>sample.sam, and then samtools is used to sort the sam file, create indexes and other operations and combine with the self-written Python script to obtain the strain identification list B, which includes the number of aligned reads, the number of aligned unique reads, the average MAPQ value of aligned reads, the average match rate of aligned reads, the genome coverage and the foldchange value indicators. The following table is an example file of List B (partial rows and columns are intercepted):

[0049] Table 3. Example of extracting some rows and columns from List B

[0050]

[0051] 7. The present invention uses a self-written Python script to calibrate list A using strain list B according to the established calibration principles. The calibration principles are as follows:

[0052] (1) Extract the bacterial species common to List A and List B to form Set S1. The detection read number of the common strains is based on the bacterial species with the larger read number in the two methods. The bacterial species in List A but not in List B form Set S2. The bacterial species not in List A but in List B form Set S3.

[0053] (2) Take the bacterial species in list A in set S1 as S1A, such as Figure 2 As shown, the score of each bacterial species is calculated according to five indicators, namely, relative abundance of bacterial species, proportion of mer aligned to species level, proportion of unique reads, foldchange value, and proportion of bacterial species reads in all reads aligned to microorganisms, and a corresponding coefficient is assigned to each indicator, wherein the proportion of bacterial species reads in all reads aligned to microorganisms is also called the absolute proportion of reads; in addition, considering that the foldchange value ranges from 0 to positive infinity, the present invention assigns a foldchange value when calculating the score of each bacterial species. When the foldchange value is positive infinity, the foldchange value is assigned to 1, when the foldchange value is greater than or equal to 5 and not positive infinity, the foldchange value is assigned to 0.5, and when the foldchange value is less than 5, the foldchange value is assigned to 0.2; the ScoreA of each bacterial species is calculated according to the following formula:

[0054]

[0055] a i is the coefficient corresponding to each indicator, X i is the score of each indicator, and m is the number of indicators.

[0056] (3) Take the bacterial species in list B in set S1 as S1B, such as Figure 3 As shown in the figure, the score of each bacterial species is calculated based on six indicators: the number of aligned reads, the number of aligned unique reads, the average MAPQ value of aligned reads, the average match rate of aligned reads, the genome coverage and the foldchange value. Similarly, the foldchange value is assigned a value. When the foldchange value is positive infinity, the foldchange value is assigned a value of 1. When the foldchange value is greater than or equal to 5 and not positive infinity, the foldchange value is assigned a value of 0.5. When the foldchange value is less than 5, the foldchange value is assigned a value of 0.2. Each indicator is assigned a corresponding coefficient, and the ScoreB of each bacterial species is calculated according to the following formula:

[0057]

[0058] Among them, bi is the corresponding coefficient of each indicator, Yi is the score of each indicator, and n is the number of indicators.

[0059] (4) Sort the ScoreA values ​​of all bacterial species in List A in Set S1 from small to large, take the first quartile as threshold 1, calculate the ScoraA of all bacterial species in S2, and the results with ScoreA values ​​above threshold 1 are credible results; sort the ScoreB values ​​of all bacterial species in List B in Set S1 from small to large, take the first quartile as threshold 2, calculate the ScoreB of all bacterial species in S3, and the results with ScoreB values ​​above threshold 2 are credible results.

[0060] (5) Summarize the reliable results in sets S2 and S3 with those in set S1 to obtain the final strain identification list.

[0061] 8. Write a Python script. The present invention uses the microbial annotation database to perform detailed annotations on all microorganisms in the final strain identification obtained in step 7. The annotation information includes but is not limited to the genus classification, biological classification, pathogenicity assessment, clinical concern level, and related information of the species in the field of epidemiology.

[0062] 9. Taking into account the specific values ​​of each indicator in the strain list, the pathogenic characteristics of the microorganisms, and the clinical manifestations and symptoms of the patient, potential pathogenic microorganisms are screened out from the annotated final strain list and sent to the patient in the form of a report.

[0063] Example 2

[0064] In order to evaluate the actual effect of the bioinformatics method for pathogenic microorganism metagenomic detection based on second-generation sequencing provided in the above Example 1, the identification ability of Escherichia coli and Shigella flexneri was compared from both qualitative and quantitative perspectives using the method of the present invention and Kraken2.

[0065] 1. Create a simulated data set

[0066] (1) 500ul negative simulated plasma samples were used for Illumina high-throughput sequencing to obtain negative control sample data, with a single-end read length of 50bp and a data volume of 20M.

[0067] (2) The reference genome sequences of Escherichia coli and Shigella flexneri were downloaded from the NCBI database. The reference genome sequence of Escherichia coli is GCF_000005845.2_ASM584v2_genomic.fna, and the reference genome sequence of Shigella flexneri is GCF_000006925.2_ASM692v2_genomic.fna.

[0068] (3) The reference genomes of Escherichia coli and Shigella flexneri were randomly broken into 50 bp sequences, and 100 sequences were taken from each, for a total of 200 sequences, with no duplication among the 200 sequences.

[0069] (4) Repeat step (3) 2000 times, add the extracted sequences to the negative control samples respectively to form 2000 simulated samples, and then divide the simulated samples into 10 groups, with 200 samples in each group.

[0070] 2. Comparison of analysis results of 2000 simulated samples

[0071] (1) Qualitative analysis results are as follows Figure 4 and Figure 5 As shown, the average accuracy of the present invention for identifying Escherichia coli and Shigella flexneri is higher than that of kraken2; the average accuracy of the present invention and kraken2 for identifying Escherichia coli is 99.99% and 99.83% respectively; the average accuracy of the present invention and kraken2 for identifying Shigella flexneri is 99.94% and 99.77% respectively.

[0072] (2) Quantitative analysis results are as follows Figure 6 and Figure 7 As shown, the detection read numbers of Escherichia coli and Shigella flexneri of the present invention are higher than those of kraken2, indicating that the species detection sensitivity of the present invention is higher than that of kraken2; and the average detection read numbers of Escherichia coli and Shigella flexneri of the present invention are 88 and 83, respectively, and the ratio is close to 1:1, which is consistent with the actual situation; while the average detection read numbers of Escherichia coli and Shigella flexneri of kraken2 are 61 and 28, respectively, and the detection ratio is inconsistent with the actual situation, indicating that kraken2 cannot accurately distinguish between the two bacteria, and a part of the reads of Shigella flexneri are allocated to Escherichia coli, resulting in the detection read number of Escherichia coli being significantly higher than that of Shigella flexneri.

Claims

1. A bioinformatics method for pathogenic microorganism metagenomic detection based on second-generation sequencing, characterized in that: The following steps are involved: S-1. Create a list of species, including bacteria, fungi, viruses, and parasites; S-2. Genome sequence download; S-3. Database construction; S-4. Quality control of next-generation sequencing data; S-5. Remove human sequences according to the human sequence database; S-6. Based on the microbial comparison database, a bacterial species identification list A is obtained based on the k-mer and LCA algorithms; S-7. Bacterial species identification list B obtained based on the bwa-mem algorithm according to the microbial reference genome database; S-8. Formulate correction principles and use List B to correct List A to obtain a more accurate list of bacterial species identification; S-9. Microorganisms whose species identification list is annotated according to the microbial annotation database; S-10. Screening of potential pathogenic bacteria; The correction principle of using list B to correct list A in step S-8 is as follows: (1) Extract the common bacterial species in List A and List B to form Set S1. The number of reads detected for the common strains is based on the number of reads of the species with more reads in the two methods. The bacterial species in List A but not in List B form Set S2. The bacterial species not in List A but in List B form Set S3. (2) Take out the bacterial species in list A in set S1, and calculate the ScoreA of each bacterial species based on the five indicators of relative abundance of bacterial species in set S1A, the proportion of mer aligned to the species level, the proportion of unique reads, the foldchange value, and the proportion of bacterial species reads in all reads aligned to microorganisms, and assign a corresponding coefficient to each indicator; The foldchange value assignment operation is: when the foldchange value is positive infinity, the foldchange value is assigned to 1; when the foldchange value is greater than or equal to 5 and not positive infinity, the foldchange value is assigned to 0.5; when the foldchange value is less than 5, the foldchange value is assigned to 0.2; The ScoreA of each bacterial species was calculated according to the following formula: Among them, a i is the corresponding coefficient of each indicator in set S1A, X i is the score of each indicator in set S1A, m is the number of indicators, and set S1A includes species relative abundance, species-level mer proportion, unique read proportion, foldchange value, and the proportion of species reads in all reads aligned to microorganisms; The corresponding coefficient assignments for each indicator are: species relative abundance coefficient is 0.25, coefficient of mer proportion aligned to species level is 0.2, coefficient of unique sequence proportion is 0.2, coefficient of foldchange value is 0.25, coefficient of bacterial species read proportion in all reads aligned to microorganisms is 0.1; (3) Take the bacterial species in list B in set S1, and calculate the ScoreB of each bacterial species based on the six indicators of the number of aligned reads, the number of aligned unique reads, the average MAPQ value of aligned reads, the average match rate of aligned reads, the genome coverage, and the foldchange value in set S1B; similarly, when the foldchange value is positive infinity, the foldchange value is assigned to 1, when the foldchange value is greater than or equal to 5 and not positive infinity, it is assigned to 0.5, and when the foldchange value is less than 5, it is assigned to 0.2; assign a corresponding coefficient to each indicator, and calculate the ScoreB of each bacterial species according to the following formula: Among them, b i is the corresponding coefficient of each indicator in set S1B, Y i is the score of each indicator in set S1B, n is the number of indicators, set S1B includes the number of aligned reads, the number of aligned unique reads, the average MAPQ value of aligned reads, the average match rate of aligned reads, genome coverage and foldchange value; The corresponding coefficients of each indicator in the set S1B are: the coefficient of the number of aligned reads is 0.25, the number of aligned unique reads is 0.25, the average MAPQ value of aligned reads is 0.1, the average match rate of reads is 0.1, the coefficient of genome coverage is 0.1, and the foldchange value is 0.2; (4) Sort the ScoreA values ​​of all bacterial species in List A in Set S1 from small to large, take the first quartile as Threshold 1, calculate the ScoreA of all bacterial species in S2, and those with ScoreA values ​​above Threshold 1 are credible results; Sort the ScoreB values ​​of all bacterial species in List B in Set S1 from small to large, take the first quartile as Threshold 2; calculate the ScoreB of all bacterial species in S3, and those with ScoreB values ​​above Threshold 2 are credible results; (5) Summarize the reliable results in sets S2 and S3 with those in set S1 to obtain the final strain identification list.

2. According to the bioinformatics method for pathogenic microorganism metagenomic detection based on second-generation sequencing as described in claim 1, it is characterized in that: The microbial reference genome database in step S-7 is composed of genome sequences marked as reference in the NCBI database in the species list created in S-1.

Images