CN112614544B - Optimization method of Kraken2 software output results and method of identifying species types in samples - Google Patents

Abstract

The invention provides a Kraken software output result optimization method, which comprises the following steps: matching sub-reads of each read in the sequencing result with species sequences in a known database, obtaining the number of kmers of the sub-reads matched with each species in each read, selecting the maximum value of the number of kmers in each read, and recording the maximum value as kmermax numbers; comparing the kmermax number with a first threshold, and removing the reads corresponding to the kmermax number when the kmermax number is smaller than or equal to the first threshold so as to filter all reads. The method can accurately optimize Kraken software output results, avoid false positive phenomenon, and can be applied to identifying species types in samples.

Description

Optimization method of Kraken2 software output results and identification of species types in samples

Technical Field

The present invention relates to the biological field. Specifically, the present invention relates to a method for optimizing the output results of Kraken2 software and a method for identifying species types in a sample.

Background Art

Metagenome refers to the sum of all biological genetic materials in a specific environment. By taking it as the research object, we can obtain the biological composition of the sample and the relationship between organisms and between organisms and the environment through sequencing analysis, functional gene screening, etc.

Metagenome sequencing, abbreviated as mNGS, stands for metagenomics next generation sequencing, which is a technology that performs mixed sequencing of the genomes of all organisms in the environment without separation.

Pathogenic microorganisms refer to microorganisms that can cause disease in humans or animals, including parasites, fungi, bacteria and viruses.

mNGS can identify all microorganisms in the sample and can identify new species, which is impossible to achieve using conventional experimental methods such as smears, biochemical identification, culture identification, and multiplex PCR-based detection technology. mNGS is widely used in clinical suspected pathogens, rare pathogens, critical illnesses, and other clinical infectious pathogen detection due to its advantages such as no need for culture, no reliance on probe sequences, and wide coverage without bias.

mNGS detects the genomes of all organisms in a sample. Before testing, the composition of microorganisms in a sample is mostly unknown. During testing, the reads obtained by sequencing need to be compared with known databases to confirm their expression pathways or species composition, etc.

BLAST, a well-known tool for sequence alignment, can obtain good alignment results for longer sequences, but the tool is not very effective when applied to short sequence data obtained by high-throughput sequencing. In addition, since mNGS often needs to match large databases when determining species, the alignment speed of BLAST limits its application in time-sensitive scenarios such as clinical testing.

Kraken2 software is a commonly used mNGS species identification tool. It constructs unique kmers for each species when building a library, and counts the kmer status of each read during alignment. Kraken2 has an extremely fast running speed, but because it does not filter the results by default, the false positive rate in the original results is very high, which can reach more than 85%.

Summary of the invention

The present invention aims to solve at least one of the technical problems existing in the prior art to a certain extent. To this end, the present invention proposes a method for optimizing the output results of Kraken2 software and a method for identifying species types in samples, which can effectively optimize the output results of Kraken2 software and accurately know the species types, so that it can be applied to the identification of species types in samples and reduce the false positive rate.

The present invention proposes a method for optimizing the output results of Kraken2 software. According to an embodiment of the present invention, the method includes: matching the sub-read segments of each read segment in the sequencing result with the species sequence in the known database, obtaining the number of kmers in each read segment that matches the sub-read segment of each species, selecting the maximum value of each kmer number in each read segment, and recording it as the kmermax number; comparing the kmermax number with a first threshold value to filter all read segments and remove read segments that do not meet the requirements.

In the optimization method according to an embodiment of the present invention, for each read in the output parameter output file of the Kraken2 software, the maximum number of kmers matched to the species is counted, and the maximum value is compared with the first threshold to filter out reads that do not meet the requirements, thereby facilitating the filtering of all reads, improving the accuracy of the results, and avoiding false positives.

According to an embodiment of the present invention, the optimization method of the output result of the above-mentioned Kraken2 software may also have the following additional technical features:

According to an embodiment of the present invention, the first threshold is 15 to 30. Thus, the reads corresponding to kmermax numbers less than or equal to 15 to 30 can be removed to avoid interference, which may lead to inaccurate detection results and false positives. It should be noted that the first threshold described in the present invention refers to any specific value between 15 and 30, such as 20.

According to an embodiment of the present invention, the method further includes: taking the remaining reads after the filtering as candidate reads, and the species corresponding to the kmermax number in each candidate read as a candidate species; for each candidate species, selecting the sum of the kmer numbers of all reads matching the candidate species, recorded as the kmersum number; comparing the kmersum number with a second threshold, and when the kmersum number is less than the second threshold, removing the candidate species corresponding to the kmersum number, so as to filter the candidate species, remove species that do not meet the requirements, and the remaining species are target species.

As mentioned above, for each species, by counting the kmermax number and comparing it with the first threshold, the reads that do not meet the requirements are filtered out. Further, the candidate species matched by the remaining reads are taken as the research object, and the sum of the kmer numbers of all reads matched by each species (i.e., the kmersum number) is counted, and the kmersum number is compared with the second threshold, so as to further filter the candidate species, remove the species that do not meet the requirements, and thus obtain the target substance. In this way, the accuracy of the results can be improved and the false positive phenomenon can be reduced.

According to an embodiment of the present invention, the second threshold is 900 to 1100. Thus, the reads corresponding to the kmersum numbers less than 900 to 1100 and the candidate species matched therewith can be removed to avoid their interference, which may lead to inaccurate detection results and false positive phenomena. It should be noted that the second threshold described in the present invention refers to any specific value between 900 and 1100, such as 1000.

According to an embodiment of the present invention, the method further includes: arranging the kmersum numbers corresponding to different candidate species in descending order, and when the previous kmersum number divided by the next kmersum number is greater than or equal to a third threshold and the next kmersum number is less than a second threshold, the candidate species corresponding to the previous kmersum number and the kmersum number before it is the target species.

By arranging the kmersum numbers corresponding to different candidate species in descending order, when the third threshold value decreases for the first time and the kmersum number is less than the second threshold value, the candidate species corresponding to the previous kmersum number and the kmersum number before it are the target species. In this way, the accuracy of the detection results can be further guaranteed and false positives can be avoided.

According to an embodiment of the present invention, the third threshold is 3 to 5. Thus, the accuracy of the detection result can be further guaranteed and false positive phenomenon can be avoided.

According to an embodiment of the present invention, the number of sequencing results is 10K-30M. It should be noted that when the sequencing results contain host data, the host data needs to be deleted in advance, and the above 10K-30M refers to the number of sequencing results after removing the host data.

For ease of understanding, the following example explains the optimization method in detail:

Example 1: Figure 1 shows a schematic diagram of the output results of a Kraken2 output parameter.

Example 2: Assume that there are a total of 10 reads in the sequencing results, recorded as read1, read2, ..., read10.

For read1, it is divided into three sub-segments, which are recorded as read1-1, read1-2, and read1-3. Read1-1, read1-2, and read1-3 are compared with the species database information in the Kraken2 software. It is found that the number of kmers matched by read1-1 and species S1 is 5; the number of kmers matched by read1-2 and species S1 is 15, and the number of kmers matched by read1-2 and species S2 is 20; the number of kmers matched by read1-2 and species S3 is 80. Therefore, the number of species kmers matched in read1 is: the number of kmers of species S1 is 20, the number of kmers of species S2 is 20, and the number of kmers matched by species S3 is 80, among which the largest number of kmers (i.e., kmermax number) is 80.

Similarly, the kmermax numbers of read2, read3, ..., read10 and their matching species are shown in the following table. Since the kmermax number of read3 is less than 10, the read result is removed. The candidate species are S1, S2, S3, S4, and S5.

The read counts corresponding to species S1 to S5 were counted, and the sum of the kmer counts of each species were calculated and arranged in descending order, as shown in the following table. It can be seen that from S4 to S5, the kmersum number dropped by more than 3 times for the first time, and 600 was less than 1000. Therefore, the data with kmersum numbers of 600 and below and the corresponding candidate species were removed. Finally, the species in the sample were determined to be S3, S1, and S4.

Species kmersum number S3 2000 S1 1200 S4 600 S5 100 S2 80

In another aspect of the present invention, the present invention proposes a method for identifying species types in a sample. According to an embodiment of the present invention, the method includes: (1) performing metagenomic sequencing on the sample to obtain sequencing results; (2) analyzing the sequencing results using Kraken2 software, and optimizing the results using the optimization method of the output results of the Kraken2 software described above, so as to determine the species types in the sample. Thus, the method according to an embodiment of the present invention can accurately identify the species types in the sample with a low false positive rate.

Additional aspects and advantages of the present invention will be given in part in the following description and in part will be obvious from the following description, or will be learned through practice of the present invention.

BRIEF DESCRIPTION OF THE DRAWINGS

The above and/or additional aspects and advantages of the present invention will become apparent and easily understood from the description of the embodiments in conjunction with the following drawings, in which:

Figure 1 is a schematic diagram of the output results of Kraken2 output parameters;

FIG2 is a schematic diagram of the relationship between different data volume sensitivities and kmermax (kmer maximum value) in one embodiment;

FIG3 is a schematic diagram of the relationship between relative specificity of different data amounts and kmermax (kmer maximum value) in one embodiment;

FIG. 4 is a schematic diagram of the relationship between the accuracy of different data amounts and kmermax (kmer maximum value) in one embodiment.

DETAILED DESCRIPTION

The scheme of the present invention will be explained below in conjunction with the embodiments. It will be appreciated by those skilled in the art that the following embodiments are only used to illustrate the present invention and should not be considered as limiting the scope of the present invention. Where specific techniques or conditions are not indicated in the embodiments, the techniques or conditions described in the literature in this area or the product specifications are used. The reagents or instruments used are not indicated by the manufacturer and are all conventional products that can be obtained commercially.

Example 1: Analysis of simulated sequencing data of 40 species

First, 40 representative species (Table 1) genomes were randomly selected from the database, covering eukaryotic, bacterial, and viral species, including some closely related species of the same genus. Then, the sequencing data simulation tool ART was used to generate double-end 75bp sequencing data and run Kraken2 to identify these species. Finally, the results of Kraken2 were screened using the two parameters kmermax and kmersum to calculate the sensitivity, relative specificity, and accuracy under different data volumes.

Sensitivity is the ratio of the number of true positive species X detected to the theoretical number of true positive species (40), that is, X/40.

Since there is no true negative value in the original data, in order to evaluate the results, we use the original results of Kraken2 as the basis, and the false positive species in the results are used as the theoretical true negative species, and their number is N. When calculating the relative specificity, these true negative species are used as references. For the number of detected species D obtained in the experimental group S, its true negative TN is N-(D-X), and the relative specificity is TN/N.

The accuracy rate is the ratio of true positives X to the number of detected species D, that is, X/D.

Table 1 Genomes of 40 representative species tested in simulation

The test data and results are shown in Tables 2 and 3. It can be seen that the use of method 3 can improve the accuracy of detection and reduce the incidence of false positives. At the same time, Figures 2 to 4 show the relationship between kmermax and detection sensitivity, relative specificity and accuracy under different data volumes. The corresponding number of reads for the data volume is shown in Table 4. It can be seen that different kmermax will significantly affect the detection sensitivity, relative specificity and accuracy. Considering comprehensively, kmermax less than or equal to 20 is selected as the standard for screening and removing inappropriate data. By removing data with kmermax less than or equal to 20, the accuracy of the results can be effectively improved, with high sensitivity, relative specificity and accuracy.

Table 2 Species with kmermax>0 (no screening) and kmermax>20 (removing kmermax<=20) (data volume: 1)

Table 3 40 species screening results (data volume: 1)

Table 4. Number of read segments corresponding to the data volume in Figures 2 to 4

Data volume Number of read segments 0.01 15560 0.02 31148 0.05 77854 0.1 155720 0.2 311408 0.5 778572 1 1557092 2 3114186

Example 2: Analysis of sequencing data of ZymoBIOMICS ^™ MICROBIAL Community Standard (Catalog No. D6300)

ZymoBIOMICS ^TM MICROBIAL Community Standard (Table 5) was added to triple distilled water for library construction and sequencing. After quality control and removal of human origin for the double-end 75bp sequencing data, Kraken2 was run for species identification and screening was performed using the kmermax and kmersum parameters. Finally, the sensitivity, relative specificity and accuracy were calculated under different data amounts.

Table 5 ZymoBIOMICS ^TM MICROBIAL Community Standard

Sensitivity is the ratio of the number of true positive species detected (X) to the theoretical number of true positive species (10), that is, X/10.

Since there is no true negative value in the original data, in order to evaluate the results, we use the original results of Kraken2 as the basis, and the false positive species in the results are used as the theoretical true negative species, and their number is N. When calculating the relative specificity, these true negative species are used as references. For the number of detected species D obtained in the experimental group S, its true negative TN is N-(D-X), and the relative specificity is TN/N.

The accuracy rate is the ratio of the number of true positive species X to the number of detected species D, that is, X/D.

The test data and results are shown in Tables 6 and 7.

Table 6 Species with kmermax>0 (no screening) and kmermax>20 (removing kmermax<=20) (data volume: 197472 reads)

Table 7 ZymoBIOMICS ^™ MICROBIAL Community Standard test results (data volume: 197472 reads)

In the screening results of Table 7, the relative specificity and accuracy of the condition group "remove kmermax <= 20 or kmersum < 1000 and the first occurrence of a 3-fold drop" are both 100%. One false negative is caused by the high similarity between Bacillus subtilis and the species of its genus. During the screening, only a small number of reads matching this species were found, resulting in a small kmersum. When the influence of the similarity of species at the genus level was considered and the kmersum was recalculated, Bacillus subtilis could be detected normally, and the sensitivity reached 100%, while the relative specificity and accuracy were not affected.

In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, without contradiction.

Although the embodiments of the present invention have been shown and described above, it is to be understood that the above embodiments are exemplary and are not to be construed as limitations of the present invention. A person skilled in the art may change, modify, replace and vary the above embodiments within the scope of the present invention.

Claims (6)

1. A method for optimizing the output results of Kraken2 software, characterized by comprising: Match the sub-reads of each read in the sequencing results with the species sequences in the known database, obtain the number of kmers in each read that matches the sub-reads of each species, and select the maximum value of each kmer number in each read, recorded as kmermax number; Comparing the kmermax number with a first threshold, and when the kmermax number is less than or equal to the first threshold, removing the reads corresponding to the kmermax number, so as to filter all the reads; The remaining reads after the filtering are used as candidate reads, and the species corresponding to the kmermax number in each candidate read is used as a candidate species; For each candidate species, the sum of the kmer numbers of all reads matching the candidate species is selected and recorded as the kmersum number; The kmersum numbers corresponding to the different candidate species are arranged in descending order, and when the previous kmersum number divided by the next kmersum number is greater than or equal to the third threshold and the next kmersum number is less than the second threshold, the candidate species corresponding to the previous kmersum number and the kmersum number before it are the target species. 2. The method according to claim 1, characterized in that the first threshold is 15~30. 3. The method according to claim 1, characterized in that the second threshold is 900~1100. 4. The method according to claim 1, characterized in that the third threshold is 3-5. 5. The method according to claim 1, characterized in that the number of sequencing results is 10K-30M. 6. A method for identifying species types in a sample, comprising: (1) Perform metagenomic sequencing on the sample to obtain sequencing results; (2) Analyzing the sequencing results using Kraken2 software, and optimizing the results using the optimization method of the Kraken2 software output according to any one of claims 1 to 5, so as to determine the species type in the sample.