A pathogen-drug resistance gene synchronous detection method and system for targeted sequencing

Abstract

The present application relates to the technical field of biomedical engineering, and particularly relates to a pathogen-drug resistance gene synchronous detection method and system based on targeted sequencing, which comprises the following steps: constructing a compound knowledge graph about pathogen and drug resistance gene; determining the overlap rate between the pathogen gene sequence and the drug resistance gene sequence of the target pathogen existing in the clinical sample to be detected based on the compound knowledge graph; analyzing the sequence inhibition sensitivity of the pathogen gene sequence and the drug resistance gene sequence of the target pathogen; constructing a collaborative state evolution model in combination with the compound knowledge graph; obtaining the optimal inhibition intensity factor of the drug resistance gene sequence by using the Kalman filtering algorithm; correcting the overlap rate by using the optimal inhibition intensity factor; and finally determining the drug resistance gene sequence and the gene abundance in all pathogen gene sequences. The present application effectively improves the detection accuracy of drug resistance genes.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of biomedical engineering, and particularly relates to a pathogen-drug resistance gene synchronous detection method and system based on targeted sequencing. BACKGROUND

[0002] In the diagnosis of infectious diseases, the identification of drug-resistant pathogens is crucial. In traditional detection methods, such as bacterial culture and biochemical detection, it usually takes several days to obtain the detection results, and in the face of multiple infections or difficult-to-culture pathogens, the detection efficiency and accuracy will be greatly reduced. In recent years, targeted sequencing technology has been introduced into the field of pathogen detection due to its high throughput and high sensitivity, which can simultaneously identify multiple pathogens and their drug resistance genes in a relatively short period of time.

[0003] In the current targeted sequencing technology, pathogens (hosts) and drug resistance genes are usually quantitatively analyzed as independent events, and the common correction methods (such as GC content-based bias correction) are also independently performed on single sequences. Due to the presence of complex PCR inhibitors (such as hematin, polysaccharides, and proteins) in clinical samples, the interference of these inhibitors on different sequences is different (sequence-specific), which can lead to inconsistent coverage depth in sequencing data (for example, the detection abundance of pathogens is high, and the detection abundance of drug resistance genes is low) for a physically connected DNA molecule (for example, a drug resistance plasmid inside a bacterium). The existing technology cannot identify this "decoupling" phenomenon caused by inhibition differences, which can easily lead to inaccurate detection results of pathogens and drug resistance genes. SUMMARY

[0004] In order to solve the above technical problem that the detection result of the drug resistance gene is not accurate due to the influence of the PCR inhibitor, the purpose of the present application is to provide a pathogen-drug resistance gene synchronous detection method and system based on targeted sequencing, and the technical scheme adopted is as follows:

[0005] In a first aspect, the present application provides a pathogen-drug resistance gene synchronous detection method based on targeted sequencing, comprising the following steps:

[0006] Based on the genomic sequence information and drug resistance gene sequence information of different pathogens, a composite knowledge graph about pathogens and drug resistance genes is constructed;

[0007] Based on the genomic sequence information of the different pathogens, the target pathogen present in the to-be-detected clinical sample and the pathogen gene sequence of the target pathogen present in the to-be-detected clinical sample are identified;

[0008] determine a set of drug-resistant gene sequences that the target pathogen is likely to have based on the composite knowledge graph, and compare the pathogen gene sequence with the drug-resistant gene sequences in the set of drug-resistant gene sequences to obtain an overlap rate between the pathogen gene sequence and the drug-resistant gene sequences;

[0009] analyze sequence inhibition sensitivity of the pathogen gene sequence of the target pathogen and the drug-resistant gene sequence, and construct a cooperative state evolution model in combination with the composite knowledge graph, the cooperative state evolution model representing a degree of PCR amplification inhibition by introducing an inhibition intensity factor as a hidden variable, and determining a state prior estimate based on the cooperative state evolution model, and updating a residual error in combination with a state observation value by using a Kalman filtering algorithm, and determining an optimal inhibition intensity factor of the drug-resistant gene sequence by iterative inverse estimation;

[0010] correct the overlap rate by using the optimal inhibition intensity factor, and determine drug-resistant gene sequences and gene abundances thereof in all the pathogen gene sequences based on an overlap rate correction value.

[0011] In combination with the first aspect, in some possible implementation manners, constructing a composite knowledge graph about pathogen and drug-resistant gene includes:

[0012] regarding each pathogen and each drug-resistant gene sequence as a node in the composite knowledge graph;

[0013] based on genomic sequence information and drug-resistant gene sequence information of different pathogens, when a drug-resistant gene sequence exists in a genomic sequence of a pathogen, constructing a directed edge between nodes of the corresponding pathogen and drug-resistant gene sequence, the directed edge being directed from the drug-resistant gene sequence to the pathogen, and determining an edge weight of the directed edge, the edge weight being used to reflect closeness of evolutionary association between the drug-resistant gene sequence and the pathogen.

[0014] In combination with the first aspect, in some possible implementation manners, determining the edge weight of the directed edge includes:

[0015] determining, for the drug-resistant gene sequence and the pathogen corresponding to the directed edge, a copy number of the drug-resistant gene sequence in the genomic sequence of the pathogen;

[0016] determining, for the drug-resistant gene sequence and the pathogen corresponding to the directed edge, an insertion position and a detection frequency of the drug-resistant gene sequence in different strain genomic sequences of the pathogen;

[0017] determining the edge weight corresponding to the directed edge based on the copy number, stability of the insertion position, and the detection frequency.

[0018] In some possible implementation manners of the first aspect, the sequence inhibition sensitivity of the pathogen gene sequence and the drug resistance gene sequence of the target pathogen is analyzed, including:

[0019] determining a core gene sequence in the pathogen gene sequence of the target pathogen;

[0020] determining GC base content and secondary structure minimum free energy of the core gene sequence and the drug resistance gene sequence;

[0021] determining an inhibition sensitivity index of the target pathogen based on the GC base content and the secondary structure minimum free energy of all the core gene sequences, and determining an inhibition sensitivity index of the drug resistance gene sequence based on the GC base content and the secondary structure minimum free energy of the drug resistance gene sequence, the inhibition sensitivity index being used to reflect an ease of the target pathogen or the drug resistance gene sequence being inhibited by an inhibitor.

[0022] In some possible implementation manners of the first aspect, the cooperative state evolution model is constructed, including:

[0023] correcting a baseline inhibition coefficient of a to-be-detected clinical sample based on the inhibition sensitivity index, to obtain an initial inhibition intensity factor of the target pathogen and the drug resistance gene sequence;

[0024] constructing a basic update equation of an inhibition intensity factor corresponding to the target pathogen and the drug resistance gene sequence in a PCR amplification process based on the initial inhibition intensity factor;

[0025] determining each candidate pathogen associated with the drug resistance gene sequence present in the to-be-detected clinical sample based on the complex knowledge graph, and determining a cooperative driving item in the PCR amplification process based on a difference in the inhibition intensity factor between the drug resistance gene sequence and the each candidate pathogen in the PCR amplification process;

[0026] determining a change rate of the inhibition intensity factor corresponding to the drug resistance gene sequence in the PCR amplification process based on the change rate of the basic update equation of the inhibition intensity factor and the cooperative driving item.

[0027] In some possible implementation manners of the first aspect, the cooperative driving item in the PCR amplification process is determined, including:

[0028] determining a difference value of the inhibition intensity factor between the each candidate pathogen and the drug resistance gene sequence in the PCR amplification process, to obtain an inhibition intensity factor difference value;

[0029] The inhibition intensity factor difference value is weighted and accumulated by using the edge weight of the directed edge between each candidate pathogen and the node corresponding to the drug resistance gene sequence in the composite knowledge graph, thereby obtaining a synergistic driving item in the PCR amplification process.

[0030] In combination with the first aspect, in some possible implementation manners, the optimal inhibition intensity factor of the drug resistance gene sequence is determined by iterative inverse estimation based on the state prior estimation determined by the synergistic state evolution model and combined with the state observation value, including:

[0031] The inhibition factor observation value and the abundance observation value of the target pathogen corresponding to each PCR amplification cycle in the Kalman filtering process are obtained, and the state observation value corresponding to each PCR amplification cycle in the Kalman filtering process is formed;

[0032] The inhibition factor prediction value of the target pathogen and the drug resistance gene sequence corresponding to each PCR amplification cycle in the Kalman filtering process is determined based on the initial inhibition intensity factor, the basic update equation of the inhibition intensity factor, and the change rate of the inhibition intensity factor;

[0033] The PCR exponential growth model of the target pathogen is corrected based on the inhibition factor prediction value of the target pathogen, and the abundance prediction value of the target pathogen corresponding to each PCR amplification cycle in the Kalman filtering process is determined;

[0034] The state prior estimation corresponding to each PCR amplification cycle in the Kalman filtering process is determined based on the inhibition factor prediction value and the abundance prediction value of the target pathogen;

[0035] The observation residual is determined based on the state prior estimation and the state observation value, and the inhibition factor prediction value of the drug resistance gene sequence is taken as a state update object, and the Kalman filtering cycle iteration is performed with the PCR amplification cycle to perform residual update, and finally the optimal inhibition intensity factor corresponding to the drug resistance gene sequence is obtained.

[0036] In combination with the first aspect, in some possible implementation manners, the overlap rate is corrected by using the optimal inhibition intensity factor, including:

[0037] The ratio of the overlap rate to the optimal inhibition intensity factor is determined;

[0038] The product of the ratio and a probe capture efficiency coefficient is determined as an overlap rate correction value.

[0039] In combination with the first aspect, in some possible implementations, based on the overlap rate correction value, the drug-resistant gene sequences and their gene abundances in all the pathogen gene sequences are determined, including:

[0040] The pathogen gene sequences with the overlap rate correction value less than or equal to a set overlap rate threshold value are removed to obtain remaining pathogen gene sequences after removal;

[0041] The remaining pathogen gene sequences after removal are subjected to repeat sequence deduplication, and reads from the same original deoxyribonucleic acid template are combined to form a UMI family, sequence consistency of the pathogen gene sequences in the UMI family is calculated, and the UMI family with the sequence consistency greater than a set consistency threshold value is retained;

[0042] The pathogen gene sequences containing chimeric sequences are removed in the retained UMI family, and finally the drug-resistant gene sequences and their gene abundances in all the pathogen gene sequences are obtained.

[0043] In a second aspect, the present application further provides a pathogen-drug resistance gene synchronous detection system for targeted sequencing, including a memory and a processor. The memory is used to store executable computer program code, and the processor is used to call and run the executable computer program code from the memory, so that the system executes the pathogen-drug resistance gene synchronous detection method for targeted sequencing in the first aspect or any one of the possible implementations of the first aspect.

[0044] This invention offers the following advantages: First, by constructing a composite knowledge graph of pathogens and drug-resistant genes to reflect the correlation between them, readable structured prior knowledge is obtained. Second, the presence of target pathogens and their original gene sequences in the clinical samples to be tested is identified to quickly pinpoint the core target. The pathogen gene sequence is then compared with drug-resistant gene sequences in the set of possible drug-resistant gene sequences of the target pathogen determined by the composite knowledge graph, yielding the overlap rate between the pathogen and drug-resistant gene sequences, thus providing initial input for subsequent fine-tuning. Next, the sequence inhibition sensitivity of the pathogen gene sequence and drug-resistant gene sequence of the target pathogen is analyzed to quantify the degree to which different genes are "naturally" susceptible to inhibition. This, combined with the composite knowledge graph, further enhances the ability to... By utilizing the biological principles of co-amplification of pathogens and drug-resistant genes, a co-state evolution model was constructed. This model introduces an inhibition strength factor as a latent variable to characterize the degree of PCR amplification suppression. Based on the state prior estimate determined by the co-state evolution model, and combined with the state observation values ​​from sequencing data, the Kalman filter algorithm was used for residual update. Through iterative inverse estimation, the optimal inhibition strength factor for the drug-resistant gene sequence was determined, thereby decoupling the interference term "inhibition factor" from the messy sequencing data. Finally, the optimal inhibition strength factor was used to correct the overlap rate between pathogen gene sequences and drug-resistant gene sequences, restoring the true biological signal, and thus determining the drug-resistant gene sequences and their abundance in all pathogen gene sequences, effectively improving the detection accuracy. Attached Figure Description

[0045] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0046] Figure 1 This is a flowchart illustrating the steps of a targeted sequencing method for simultaneous detection of pathogen-drug resistance genes according to an embodiment of the present invention.

[0047] Figure 2 This is a flowchart illustrating the steps involved in constructing a composite knowledge graph about pathogens and drug resistance genes, as described in an embodiment of the present invention.

[0048] Figure 3 This is a flowchart illustrating the steps of analyzing the sequence inhibition sensitivity of pathogen gene sequences and drug resistance gene sequences of a target pathogen according to an embodiment of the present invention.

[0049] Figure 4A step flow chart for constructing a collaborative state evolution model of an embodiment of the present application;

[0050] Figure 5 A step flow chart for determining an optimal inhibitory strength factor of a drug resistance gene sequence of an embodiment of the present application;

[0051] Figure 6 A step flow chart for determining a drug resistance gene sequence and its gene abundance in all pathogen gene sequences of an embodiment of the present application;

[0052] Figure 7 A structural schematic diagram of a pathogen-drug resistance gene synchronous detection system of an embodiment of the present application. DETAILED DESCRIPTION

[0053] To clearly illustrate the technical features of the present application, the present application will be described in detail below with reference to specific embodiments and in conjunction with the accompanying drawings.

[0054] Embodiments of the present application will be described in more detail by referring to the accompanying drawings. Although certain embodiments of the present application are shown in the drawings, it is understood that the present application can be implemented in various forms and should not be construed as being limited to the embodiments set forth herein, but rather these embodiments are provided so as to more completely and thoroughly understand the present application. It is understood that the drawings and embodiments of the present application are for exemplary purposes only and are not intended to limit the scope of the present application.

[0055] It should be understood that each of the steps recited in the method embodiments of the present application can be performed in different orders and / or in parallel. In addition, the method embodiments can include additional steps and / or omit performing the steps shown. The scope of the present application is not limited in this respect.

[0056] The term "comprising" and variations thereof as used herein are open-ended, that is, "comprising but not limited to". The term "based on" is "based, at least in part, on". The term "one embodiment" means "at least one embodiment"; the term "another embodiment" means "at least one additional embodiment"; the term "some embodiments" means "at least some embodiments". Related definitions will be given in the description below.

[0057] It should be noted that the concepts of "first", "second", etc. mentioned in the present application are only used to distinguish different devices, modules or units, and are not used to limit the order or interdependence of the functions performed by these devices, modules or units.

[0058] Although the operations or steps are described in a particular order in the figures, it should be understood that the operations or steps can be performed in other sequences or in a serial or parallel manner without departing from the scope of the present application. In some embodiments, some of the operations or steps can be performed concurrently.

[0059] Meanwhile, it can be understood that the data involved in the technical solutions of the present application (including but not limited to the data itself, the acquisition or use of the data) should comply with the requirements of the corresponding laws, regulations and relevant provisions. Unless otherwise defined, all technical and scientific terms used in the present application have the same meaning as understood by those skilled in the art to which the present application belongs. Meanwhile, in all division operations and logarithmic operations involved in the present application, in order to prevent the calculation from collapsing or producing invalid values due to the denominator being zero or the input being zero, a protection mechanism is adopted. The implementation means of the protection mechanism can be reasonably set according to the actual situation, for example, when the denominator item of the division operation or the real number item of the logarithmic function is zero, a protection parameter with the same dimension or dimensionless as the denominator item or the real number item can be added, the value of the protection parameter can be a very small value greater than zero, so as to ensure the robustness and implementability of the algorithm under extreme working conditions. In addition, the normalization function mentioned in the present application adopts maximum and minimum value normalization unless otherwise specified, in order to normalize the normalization result to the [0, 1] interval or other continuous intervals. Among them, the maximum value and the minimum value used in the maximum and minimum value normalization can be obtained according to the actual situation, for example, when a plurality of values can be obtained during implementation and the relationship between different values needs to be compared, the maximum value and the minimum value can be obtained by counting a plurality of values, and when only a single value can be obtained during implementation, the maximum value and the minimum value can be obtained based on a large amount of historical experimental data or prior data obtained in advance.

[0060] Next, a pathogen-drug resistance gene synchronous detection method and system for targeted sequencing provided by an embodiment of the present application will be described in detail with reference to the accompanying drawings.

[0061] Figure 1 The basic flowchart of a pathogen-drug resistance gene synchronous detection method for targeted sequencing provided by an embodiment of the present application is shown in FIG. 1, which specifically includes the following steps: Figure 1

[0062] Step S100: Based on the genomic sequence information and drug resistance gene sequence information of different pathogens, a composite knowledge graph about pathogens and drug resistance genes is constructed.

[0063] ​The genomic sequence information of the disclosed pathogen and the disclosed drug-resistant gene sequence information are collected from a public database, and the drug-resistant gene sequence information is compared with the genomic sequence information of the pathogen to construct a complex knowledge graph about the pathogen and the drug-resistant gene. The complex knowledge graph reflects the pathogen and the drug-resistant gene as a "drug resistance ecological map", reflects the correlation between the pathogen and the drug-resistant gene, for example, which drug-resistant genes are carried by a pathogen, and the co-occurrence relationship between different drug-resistant genes, etc., providing readable structured prior knowledge for subsequent identification of the relationship between the pathogen and the drug-resistant gene.

[0064] Step S200: Based on the genomic sequence information of different pathogens, the target pathogen present in the clinical sample to be detected and the pathogen gene sequence of the target pathogen present in the clinical sample to be detected are identified.

[0065] The pathogen present in the clinical sample to be detected is determined as the target pathogen, and all DNA sequences of the sequencing data of the clinical sample to be detected are sequenced and obtained. According to the genomic sequence information of the disclosed pathogen, the pathogen gene sequence of the target pathogen present in the clinical sample to be detected is determined.

[0066] Step S300: Based on the complex knowledge graph, a set of drug-resistant gene sequences that the target pathogen may exist is determined, and the pathogen gene sequence is compared with the drug-resistant gene sequence in the set of drug-resistant gene sequences to obtain the overlap rate between the pathogen gene sequence and the drug-resistant gene sequence.

[0067] In order to identify whether there is a drug-resistant gene in the pathogen gene sequence, according to the complex knowledge graph obtained above, a set of fixed drug-resistant genes that may exist in the target pathogen is determined, that is, all drug-resistant gene sequence nodes that exist directed edges with the target pathogen node in the complex knowledge graph are determined, and the set of drug-resistant gene sequences corresponding to these drug-resistant gene sequence nodes is used as the set of drug-resistant gene sequences. Using the sequence comparison tool, the pathogen gene sequence of the target pathogen is compared with any drug-resistant gene sequence in the above drug-resistant gene set, and the overlap rate between the gene sequences is calculated, which can be used as a basis for verifying whether the drug-resistant gene sequence exists in the pathogen gene sequence.

[0068] It should be understood that a gene sequence is usually represented by a sequence of nucleotides in a straight line. For example, ATGCGTAC is a simple gene fragment, and the above overlap rate refers to the proportion of the number of matched nucleotides in the length of the drug-resistant gene sequence after the two gene fragments are matched by sequence alignment (matching the same nucleotides). The core is to measure the "homology" of the sequence (the higher the similarity, the stronger the sequence homology). For example, fragment 1: ATGCGTAC (length 8); fragment 2: ATGCGTAG (length 8); the number of matched nucleotides = 7, the total length of the alignment = 8, and the overlap rate is 7 / 8. For fragments of different lengths, fragment X: ATGCGTAC (length 8); fragment Y: CGTAC (length 5); the number of matches = 5, the total length of the drug-resistant gene sequence = 5, and the overlap rate (similarity) = 5 / 5 = 1.

[0069] Step S400: Analyze the sequence inhibition sensitivity of the pathogen gene sequence and the drug-resistant gene sequence of the target pathogen, and construct a collaborative state evolution model combined with the composite knowledge graph. The collaborative state evolution model introduces an inhibition intensity factor as a hidden variable to represent the degree of PCR amplification inhibition. The Kalman filtering algorithm is used to determine the state priori estimation based on the collaborative state evolution model, and the residual error is updated combined with the state observation value. Through iterative inverse estimation, the optimal inhibition intensity factor of the drug-resistant gene sequence is determined.

[0070] For the overlap rate between the pathogen gene sequence and the drug-resistant gene sequence in the above determined set of pathogen gene sequences and drug-resistant gene sequences, if the overlap rate is high, it can be verified that there is a certain drug-resistant gene sequence fragment in the pathogen gene sequence, and the abundance of these verified drug-resistant genes can be calculated to obtain the detection result. However, since the clinical samples (especially blood, sputum, and pus) may contain hemoglobin, polysaccharides, proteins, and other PCR inhibitors, these inhibitors have different effects on the amplification efficiency of different gene sequences, resulting in systematic deviation in the consistency of pathogen genome and drug-resistant gene detection between different samples. Specifically, if multiple samples contain the same pathogen and its carried drug-resistant genes, theoretically these genes should be detected in all samples (i.e., complete overlap), but in actual sequencing data, there is an inconsistent phenomenon of partial sample detection and partial sample missed detection, resulting in a lower gene overlap rate than expected. This is the result of the combined effect of differential interference of inhibitors on amplification and sequencing randomness.

[0071] The embodiment considers decoupling and compensating the overlap ratio error by using the state tracking mechanism of Kalman filtering to eliminate the inhibitor interference. The core advantage of Kalman filtering is to recursively estimate the optimal state value from the noisy observation by using the dynamic evolution law of the system and the statistical characteristics of the observation data. However, since the Kalman filtering assumes that the state evolution is a Markov process, that is, the system state needs to be "independent evolution", such as the change of A is not directly bound to B, while the pathogen-drug resistance gene system has a special constraint structure, the copy number of drug resistance gene is not independent evolution, which is completely constrained by the number of pathogen carrying it, that is, the more pathogen, the more drug resistance gene. This binding relationship cannot be handled by standard Kalman filtering, so it needs to be improved and adjusted.

[0072] Therefore, the embodiment analyzes the sequence inhibition sensitivity of the pathogen gene sequence and the drug resistance gene sequence of the target pathogen, combines the composite knowledge graph, uses the biological characteristics of the synergistic amplification of the pathogen and the drug resistance gene, introduces the inhibition intensity factor as a hidden variable to represent the degree of PCR amplification inhibition, constructs a synergistic state evolution model corresponding to the pathogen and the drug resistance gene sequence, and then uses the Kalman filtering algorithm to determine the state priori estimation based on the synergistic state evolution model, and combines the state observation value to update the residual error. Through iterative reverse estimation, the optimal inhibition intensity factor of the drug resistance gene sequence is determined. Subsequently, by using the optimal inhibition intensity factor to correct the overlap ratio between the pathogen gene sequence and the drug resistance gene sequence, the accuracy of the detection of the drug resistance gene sequence in the pathogen gene sequence can be effectively improved.

[0073] Step S500: using the optimal inhibition intensity factor to correct the overlap ratio, and determining the drug resistance gene sequence and its gene abundance in all pathogen gene sequences based on the overlap ratio correction value.

[0074] The above optimal inhibition intensity factor of the drug resistance gene sequence output by the Kalman filtering is used as the core to correct the overlap ratio between the pathogen gene sequence and the drug resistance gene sequence to eliminate the measurement interference caused by the inhibitor. Further, based on the overlap ratio correction value, the drug resistance gene sequence and its gene abundance in all pathogen gene sequences are finally accurately determined.

[0075] Based on the above technical scheme, the complex knowledge graph about the pathogen and the drug-resistant gene is constructed to reflect the correlation between the pathogen and the drug-resistant gene, so as to obtain readable structured prior knowledge; the target pathogen and the pathogen gene sequence existing in the target pathogen are identified in the to-be-detected clinical sample to lock the core target, and the pathogen gene sequence is compared with the drug-resistant gene sequence in the drug-resistant gene sequence set that may exist in the target pathogen based on the complex knowledge graph, so as to obtain the overlap rate between the pathogen gene sequence and the drug-resistant gene sequence, thereby providing original input for subsequent fine correction; the sequence inhibition sensitivity of the pathogen gene sequence and the drug-resistant gene sequence of the target pathogen is analyzed, so as to quantify the degree to which different genes are naturally easy to be inhibited, and a cooperative state evolution model is constructed by combining the complex knowledge graph and using the biological law of cooperative amplification of the pathogen and the drug-resistant gene, the cooperative state evolution model introduces an inhibition intensity factor as a hidden variable to represent the degree of PCR amplification inhibition; the state prior estimation determined based on the cooperative state evolution model is combined with the state observation value of the sequencing observation data, the Kalman filtering algorithm is used for residual error updating, the optimal inhibition intensity factor of the drug-resistant gene sequence is determined through iterative reverse estimation, so that the interference term of the inhibition factor is decoupled from the chaotic sequencing data; the optimal inhibition intensity factor is used to correct the overlap rate between the pathogen gene sequence and the drug-resistant gene sequence, and the real biological signal is restored, and then the drug-resistant gene sequence and the gene abundance in all pathogen gene sequences are determined, so as to effectively improve the detection accuracy.

[0076] In a possible implementation manner, as shown in FIG. 1, the step S100 of constructing the complex knowledge graph about the pathogen and the drug-resistant gene includes: Figure 2

[0077] Step S101: Each pathogen and each drug-resistant gene sequence is taken as a node in the complex knowledge graph.

[0078] Step S102: Based on the genomic sequence information and the drug-resistant gene sequence information of different pathogens, when the drug-resistant gene sequence exists in the genomic sequence of the pathogen, a directed edge is constructed between the nodes of the corresponding pathogen and the drug-resistant gene sequence, and the edge weight of the directed edge is determined, and the directed edge is pointed from the drug-resistant gene sequence to the pathogen.

[0079] ​Specifically, for the construction process of the composite knowledge graph, each pathogen species is taken as a node, and each drug-resistant gene sequence is also taken as a node. When a certain drug-resistant gene sequence has been confirmed to exist in the genome of a certain pathogen, a directed edge is established between the corresponding pathogen node and the drug-resistant gene sequence node, and the direction of the directed edge is “drug-resistant gene sequence node→pathogen node”. The directed edge represents the subordinate relationship between the drug-resistant gene sequence and the pathogen, that is, the drug-resistant gene sequence belongs to one of the genomes of the pathogen. At the same time, the edge weight of the directed edge is determined, which is used to reflect the closeness of the evolutionary association between the drug-resistant gene and the pathogen.

[0080] Further, in a possible implementation, the edge weight of the directed edge can be determined in the following manner:

[0081] First, for the drug-resistant gene sequence and the pathogen corresponding to the directed edge, the copy number of the drug-resistant gene sequence in the genome sequence of the pathogen is determined. The copy number (Copy Number) refers to how many copies of the drug-resistant gene are contained within a single pathogen (such as a bacterial cell).

[0082] Second, for the drug-resistant gene sequence and the pathogen corresponding to the directed edge, the insertion position and the detection frequency of the drug-resistant gene sequence in the different strain genome sequences of the pathogen are determined. That is, for the drug-resistant gene sequence node and the pathogen node corresponding to each directed edge, the multiple strain genome sequences of the pathogen species are counted, the insertion position of the drug-resistant gene sequence in different strains is compared (which can be represented by the sequence number of the inserted gene sequence), and when the insertion position is more stable, it means that the drug-resistant gene sequence is more conservative in the genome position of the corresponding pathogen. At the same time, the proportion of the drug-resistant gene sequence in the sequenced strains of the pathogen is counted to obtain the detection frequency of the drug-resistant gene sequence in the pathogen strain. The detection frequency refers to the proportion of strains in which the drug-resistant gene sequence is detected in all publicly available and sequenced strain genomes of the pathogen.

[0083] Finally, based on the copy number, the stability of the insertion position, and the detection frequency, the edge weight corresponding to the directed edge is determined. That is, for the drug-resistant gene sequence node and the pathogen node corresponding to each directed edge, the copy number is normalized to the range [0, 1] using the maximum-minimum value normalization function to obtain the copy number normalized value, and is denoted as C; the standard deviation of the position coordinates of all insertion positions is calculated, which is used to reflect the conservation degree of the drug-resistant gene position. The smaller the standard deviation, the more conservative it is, and vice versa. The drug-resistant gene may be obtained through horizontal gene transfer, and the position variability is large. The standard deviation is normalized to the range [0, 1] using the maximum-minimum value normalization function to obtain the standard deviation normalized value, and is denoted as D. The detection frequency is normalized to the range [0, 1] using the maximum-minimum value normalization function to obtain the detection frequency normalized value, and is denoted as P. Further, based on the copy number normalized value C, the standard deviation normalized value D, and the detection frequency normalized value P, the edge weight corresponding to the directed edge is calculated ; wherein , and respectively represent the reference weights of the copy number, the standard deviation, and the detection frequency, which are determined according to the influence degree of the copy number, the standard deviation, and the detection frequency on the relationship between the drug-resistant gene sequence and the pathogen. The higher the influence degree, the higher the reference weight, and . For example, for the drug resistance monitoring scene, the position conservation degree reflects whether the gene is fixed in the genome. The smaller the standard deviation (high conservation degree), the more likely it is that the drug-resistant gene is a characteristic of the pathogen, rather than a temporary horizontal transfer. Increasing its weight helps to filter out more stable association relationships and reduce interference caused by plasmid exchange. The standard deviation can be set to , , . The edge weight corresponding to the directed edge is larger, indicating that the association between the drug-resistant gene and the pathogen is closer and more stable.

[0084] In addition, since there is also a co-occurrence relationship between drug-resistant gene sequences, when multiple drug-resistant gene sequences often appear simultaneously on the same mobile genetic element, an undirected edge is also established between the nodes corresponding to these drug-resistant gene sequences, and the co-occurrence frequency of these drug-resistant gene sequences in the known pathogen genome is taken as the edge weight. Thus, a complex knowledge graph is constructed, which includes pathogen nodes, drug-resistant gene sequence nodes, and describes the association relationship between them (directed graph) and the co-occurrence relationship between drug-resistant genes (undirected graph).

[0085] Based on the above technical scheme, by converting the biological pathogen-drug resistance gene attribution relationship into a quantitative topological constraint weight, differentiated synergistic evolution stiffness is provided for subsequent algorithms, effectively solving the ambiguity of tracing the source of drug resistance genes in mixed infections, and preventing over-correction of unstable associated genes, thereby significantly improving the accuracy of complex sample detection.

[0086] In one possible implementation, as shown in FIG. 4, the sequence inhibition sensitivity of the pathogen gene sequence and the drug resistance gene sequence of the target pathogen in step S400 includes: Figure 3

[0087] Step S401: Determine the core gene sequence in the pathogen gene sequence of the target pathogen.

[0088] Specifically, for all pathogen gene sequences (pathogens as complete genomes, containing tens to thousands of gene sequences) obtained by sequencing in the to-be-detected clinical sample, gene sequences with a length of 500bp-5000bp (belonging to the mainstream length range of clinical PCR amplification target genes, too short sequence characteristics are not representative, and too long is easy to cause low amplification efficiency) are retained; for these retained gene sequences, gene sequences with a sequencing coverage of ≥90% are further retained to avoid GC content and free energy calculation bias caused by incomplete sequences; then a plurality of sequence clusters are obtained through sequence similarity clustering (CD-HIT algorithm, similarity threshold 95%), only one representative gene sequence is retained in the same clustering cluster to avoid feature repetition caused by highly homologous genes; and finally, 10 or so core gene sequences are selected from the representative gene sequences in order from high to low according to the copy amount.

[0089] Step S402: Determine the GC base content and the minimum free energy of the secondary structure of the gene sequence for the core gene sequence and the drug resistance gene sequence.

[0090] Specifically, different DNA sequences have different sensitivities to inhibitors, and sequences with high GC content and stable secondary structure are more easily blocked by inhibitors to prevent polymerase extension. Therefore, the GC base content and the minimum free energy of the secondary structure of each core gene sequence are calculated, and the GC base content and the minimum free energy of the secondary structure of the drug resistance gene sequence are also calculated. Sequences with more GC bases are more easily disturbed, and structures with small free energy are more stable and are more easily targeted by inhibitors.

[0091] ​Step S403: determining the inhibition sensitivity index of the target pathogen based on the GC base content and the minimum free energy of the secondary structure of all core gene sequences, and determining the inhibition sensitivity index of the drug-resistant gene sequence based on the GC base content and the minimum free energy of the secondary structure of the drug-resistant gene sequence, which reflects the ease of inhibition of the target pathogen or the drug-resistant gene sequence by the inhibitor.

[0092] Specifically, for each core gene sequence, the absolute value of the minimum free energy of the secondary structure and the GC base content corresponding to the sequence are normalized to the range of [0, 1] respectively by using the maximum-minimum value normalization function, and the average of the two normalized values is calculated as the inhibition sensitivity index of each core gene sequence. The average of the inhibition sensitivity indexes of all core gene sequences is calculated, and the average is taken as the inhibition sensitivity index of the corresponding target pathogen. At the same time, for the drug-resistant gene sequence, the absolute value of the minimum free energy of the secondary structure and the GC base content corresponding to the sequence are also normalized to the range of [0, 1] respectively by using the maximum-minimum value normalization function, and the average of the two normalized values is calculated as the inhibition sensitivity index of the drug-resistant gene sequence.

[0093] Based on the above technical solution, by determining the core gene sequences in the pathogen gene sequences of the target pathogen and identifying the GC base content and the minimum free energy of the secondary structure in the core gene sequences and the drug-resistant gene sequences, the inhibition sensitivity of the target pathogen and the drug-resistant gene sequence can be accurately quantified.

[0094] In one possible implementation, as shown in Figure 4 , the step S400 of constructing the cooperative state evolution model includes:

[0095] Step S411: based on the inhibition sensitivity index, the baseline inhibition coefficient of the clinical sample to be detected is corrected to obtain the initial inhibition intensity factor of the target pathogen and the drug-resistant gene sequence.

[0096] Specifically, the initial value of the inhibition intensity factor of the sequence with high inhibition sensitivity (rich in GC and stable secondary structure) should be smaller but the decay rate should be faster, and the initial value of the sequence with low inhibition sensitivity should be larger but the decay should be slower. The baseline inhibition coefficient of different clinical samples (such as blood, sputum) is obtained , which can be obtained by calculating the average inhibition level of historical sample data, such as 0.8 for blood samples and 0.6 for sputum / pus.

[0097] According to the inhibition sensitivity index of the target pathogen and the baseline inhibition coefficient of the clinical sample to be detected to which it belongs, the initial inhibition intensity factor of the target pathogen is calculated ; wherein:​ represents the cycle time consumption; represents the balance coefficient representing the weight for balancing the baseline inhibition and the sequence-specific inhibition, which is set empirically ; represents the inhibition sensitivity index of the target pathogen.

[0098] Similarly, according to the inhibition sensitivity index of the drug resistance gene sequence and the baseline inhibition coefficient of the clinical sample to be detected to which the drug resistance gene sequence belongs, the initial inhibition intensity factor of the drug resistance gene sequence is calculated ; wherein: represents the inhibition sensitivity index of the drug resistance gene sequence.

[0099] Step S412: Based on the initial inhibition intensity factor, a basic update equation of the inhibition intensity factor of the target pathogen and the drug resistance gene sequence in the PCR amplification process is constructed.

[0100] The inhibition intensity factor is a hidden variable in the state space, which dynamically changes with the PCR amplification cycle: a smaller value (less than 1 indicates inhibition, which will cause the observed value to be lower than the true value) is taken when the inhibitor concentration is high at the beginning of amplification, and approaches 1 (indicating no inhibition) when the inhibitor is diluted as the amplification proceeds; each target pathogen and each drug resistance gene sequence has its specific inhibition intensity factor.

[0101] Specifically, according to the initial inhibition intensity factor of the target pathogen and the drug resistance gene sequence, a basic update equation of the inhibition intensity factor of the target pathogen and the drug resistance gene sequence in the PCR amplification process is constructed:

[0102]

[0103] In the formula: represents the update value of the inhibition intensity factor of the target pathogen as the PCR amplification proceeds to cycle number t; represents the update value of the inhibition intensity factor of the drug resistance gene sequence as the PCR amplification proceeds to cycle number t; represents the initial inhibition intensity factor of the target pathogen ; represents the initial inhibition intensity factor of the drug resistance gene sequence ; represents the decay coefficient of the clinical sample to be detected, which is set according to the type of inhibitor in the clinical sample to be detected, such as 0.05 / cycle. When multiple inhibitors coexist, the most difficult to eliminate (slowest decay) will determine the inhibition tailing time of the entire reaction, so the minimum value of the decay coefficients of all inhibitors is selected as the decay coefficient of the clinical sample to be detected.

[0104] In the above formula, As the cycle approaches 0, the first term or (the initial inhibition effect) gradually disappears, and the second term or (the non-inhibition state compensation term) gradually approaches or , showing the actual law that the inhibition effect decays with the cycle.

[0105] Step S413: Based on the composite knowledge graph, determine each candidate pathogen associated with the drug-resistant gene sequence present in the to-be-detected clinical sample, and based on the difference in the inhibition intensity factor between the drug-resistant gene sequence and each candidate pathogen in the PCR amplification process, determine the synergistic driving term in the PCR amplification process.

[0106] In the composite knowledge graph, extract the directed edge of the pathogen-drug-resistant gene, thereby determining all pathogen types associated with the drug-resistant gene sequence (i.e., there is a directed edge between the pathogen node and the drug-resistant gene sequence node), and obtaining all pathogen types present in the to-be-detected clinical sample, which constitute the candidate pathogen of the drug-resistant gene sequence. It should be understood that when the number of candidate pathogens is large, only the top K (such as Top 5) candidate pathogens in the abundance (Reads Count) can be selected to participate in the subsequent synergistic evolution calculation.

[0107] Consider the standard Kalman filter to assume that the state evolves independently in advance, but the pathogen and the drug-resistant gene it carries (the same DNA molecule) should be subject to the same inhibition, and if the difference in the inhibition intensity factor between the two is large, it does not conform to logic. Therefore, in this embodiment, the inhibition intensity factor of the drug-resistant gene sequence and the pathogen is analyzed by analyzing the difference in the inhibition intensity factor between the drug-resistant gene sequence and each candidate pathogen in the PCR amplification process, defining the synergistic driving term of the synergistic evolution of the two, which is used to replace the "independent state transition equation" of the standard Kalman filter, so that the change of the inhibition intensity factor of the drug-resistant gene is simultaneously driven by the sequence characteristics of the drug-resistant gene and the associated pathogen factor, ensuring the synergy between the two.

[0108] Further, in a possible implementation, determining the synergistic driving term in the PCR amplification process includes: determining the difference value of the corresponding inhibition intensity factor of each candidate pathogen and the drug-resistant gene sequence in the PCR amplification process, to obtain the inhibition intensity factor difference value; using the edge weight of the directed edge between the corresponding nodes of each candidate pathogen and the drug-resistant gene sequence in the composite knowledge graph, to weight and accumulate the inhibition intensity factor difference value, thereby obtaining the synergistic driving term in the PCR amplification process.

[0109] Step S414: Determine the change rate of the inhibition intensity factor corresponding to the drug resistance gene sequence in the PCR amplification process based on the change rate of the basic update equation of the inhibition intensity factor and the synergistic driving term.

[0110] Specifically, the change rate of the inhibition intensity factor corresponding to the drug resistance gene sequence in the PCR amplification process is equal to the addition value of the exponential decay term of its own inhibition intensity factor and the synergistic driving term, at which time the change rate of the inhibition intensity factor of the drug resistance gene sequence is:

[0111]

[0112] In the formula: represents the change rate of the inhibition intensity factor of the drug resistance gene sequence ; represents the updated value of the inhibition intensity factor of the drug resistance gene sequence as the PCR amplification proceeds to cycle number t; represents the edge weight of the directed edge between the node of the drug resistance gene sequence and the node of the candidate pathogen ; represents the updated value of the inhibition intensity factor of the candidate pathogen as the PCR amplification proceeds to cycle number t, which can be determined according to the above basic update equation of the inhibition intensity factor corresponding to the target pathogen in the PCR amplification process; represents the updated value of the inhibition intensity factor of the drug resistance gene sequence as the PCR amplification proceeds to cycle number t; represents the total number of candidate pathogens.

[0113] In the above formula, the first term is the derivative equation of the dynamic update equation of the inhibition intensity factor of the drug resistance gene sequence, and the first term is the basic update rate of the drug resistance gene inhibition intensity; the second term, the synergistic driving term, is the weighted sum of the difference between the inhibition intensity factors of all candidate pathogens and the inhibition intensity factor of the drug resistance gene sequence, i.e. wherein constitutes a negative feedback, when the synergistic driving term is negative, causing to decrease, the growth rate slows down or decreases, and vice versa the growth rate accelerates, ultimately realizing the synchronization of the inhibition intensity factor of each candidate pathogen and the inhibition intensity factor of the drug resistance gene sequence .

[0114] Based on the above technical scheme, the synergistic driving term is constructed in the PCR amplification process to realize the synergistic evolution mechanism, which not only retains the attenuation law of the inhibition effect of the drug-resistant gene sequence itself, but also converts the topological constraint of the knowledge graph into a dynamic coupling, so that the inhibition intensity factors of the pathogen-drug-resistant gene pair derived from the same DNA molecule automatically tend to be consistent.

[0115] Further, the Kalman filtering algorithm is used to determine the optimal inhibition intensity factor of the drug-resistant gene sequence based on the state prior estimation determined by the synergistic state evolution model and the residual error update combined with the state observation value through iterative reverse estimation.

[0116] In a possible implementation manner, as shown in Figure 5 the step S400 of using the Kalman filtering algorithm to determine the optimal inhibition intensity factor of the drug-resistant gene sequence based on the state prior estimation determined by the synergistic state evolution model and the residual error update combined with the state observation value through iterative reverse estimation includes the following steps.

[0117] Step S421: Obtain the inhibition factor observation value and the abundance observation value of the target pathogen corresponding to each PCR amplification cycle in the Kalman filtering process, and form the state observation value corresponding to each PCR amplification cycle in the Kalman filtering process.

[0118] Specifically, in each PCR amplification cycle, two types of information need to be tracked at the same time:

[0119] The first type is the target detection quantity: the measured abundance of the target pathogen, that is, the measured copy number. The abundance of the drug-resistant gene and the abundance of the pathogen are bound. The more pathogen genes, the more drug-resistant genes. In the case of uncertain drug-resistant genes, the abundance of the drug-resistant gene is indirectly observed according to the abundance of the pathogen amplification.

[0120] Among them, the theoretical amplification abundance (copy number) of the target pathogen under the experimental environment, that is, in the ideal environment without inhibitory amplification substances, when the target pathogen undergoes different PCR amplification cycles can be obtained by traditional bacterial culture and biochemical detection methods (which takes several days and is pre-prepared experimental data).

[0121] The second type is the interference correction quantity: the real-time inhibition intensity factor of the target pathogen, which is used to quantify the degree of inhibition of the target pathogen in the amplification process.

[0122] Wherein, whether the amplification process is inhibited can be determined by calculating the difference between the theoretical amplification abundance of the target pathogen and the measured abundance obtained above, the theoretical amplification abundance refers to the number that should be reached when the target pathogen undergoes different PCR amplification cycles under ideal conditions without inhibition (efficiency 100%), at this time the real-time inhibition intensity factor of the target pathogen = difference / theoretical amplification abundance, complete non-inhibition is recorded as zero, and complete non-amplification is recorded as 1.

[0123] The amplification abundance of the target pathogen undergoing different PCR amplification cycles and the real-time inhibition intensity factor of the target pathogen obtained above are taken as the observation values of the current cycle state corresponding to different PCR amplification cycles.

[0124] It should be understood that, since the amplification process cannot be interrupted, in a specific example, a small amount of discrete real observation values can be obtained in advance by a limited number of discrete termination experiments (a conventional experimental method in the art: artificially terminate the PCR reaction at different cycle numbers t, and obtain the abundance of the cycle by sequencing), and then the observation values of each cycle are generated by interpolation fitting. In another specific example, based on the abundance in the final sequencing result and the theoretical amplification efficiency E, the observation values of each cycle are obtained by the formula , t is the cycle number, is the abundance at cycle number t, is the initial template amount, and the virtual observation sequence of the reconstructed abundance is obtained, and then the observation values of each cycle are obtained.

[0125] At the same time, two types of basic noise are introduced: process noise (describing the random fluctuations of PCR amplification) and observation noise (reflecting the inherent error of the sequencing platform), both of which conform to Gaussian distribution. The values of these two types of noise can be determined by statistical analysis of historical experimental data, which is the basic composition of Kalman filtering, and will not be repeated here.

[0126] Step S422: Based on the initial inhibition intensity factor, the basic update equation of the inhibition intensity factor, and the change rate of the inhibition intensity factor, the prediction value of the inhibition factor of the target pathogen and the drug-resistant gene sequence corresponding to each PCR amplification cycle in the Kalman filtering process is determined.

[0127] Specifically, the state prediction is performed based on the co-evolution, and the state prediction includes the inhibition factor prediction and the abundance prediction. In the inhibition factor prediction, the initial inhibition strength factor of the target pathogen in the constructed co-state evolution model is used, and the basic update equation of the corresponding inhibition strength factor of the target pathogen in the PCR amplification process is combined to determine the inhibition factor prediction value of the target pathogen corresponding to each PCR amplification cycle in the Kalman filtering process; meanwhile, the initial inhibition strength factor of the drug resistance gene sequence is used, and the basic update equation of the corresponding inhibition strength factor of the drug resistance gene sequence in the PCR amplification process and the change rate of the inhibition strength factor are combined to determine the inhibition factor prediction value of the drug resistance gene sequence corresponding to each PCR amplification cycle in the Kalman filtering process, so that the inhibition factor change of the drug resistance gene is simultaneously driven by the sequence characteristics and the associated pathogen, and the sequence inhibition effect from the same DNA molecule source is ensured to be consistent.

[0128] Step S423: Based on the inhibition factor prediction value of the target pathogen, the PCR exponential growth model of the target pathogen is corrected to determine the abundance prediction value of the target pathogen corresponding to each PCR amplification cycle in the Kalman filtering process.

[0129] Specifically, in the abundance prediction in the state prediction, the existing PCR exponential growth model is used, where is the initial template amount, E is the amplification efficiency (generally assumed to be 1), t is the cycle number, is the abundance at the cycle number t. In the existing PCR exponential growth model , the inhibition factor prediction value X2 of the target pathogen is introduced to regulate the amplification efficiency, and the stronger the inhibition, the slower the abundance growth, so that the corrected PCR exponential growth model is obtained. Using the corrected model, the predicted abundance of the target pathogen at different PCR amplification cycles can be determined.

[0130] Step S424: Based on the inhibition factor prediction value and the abundance prediction value of the target pathogen, the state prior estimate corresponding to each PCR amplification cycle in the Kalman filtering process is determined.

[0131] Specifically, the inhibition factor prediction value and the abundance prediction value of the target pathogen and the drug resistance gene sequence are integrated together to form a vector, and the vector is taken as the state prior estimate of the current cycle.

[0132] Step S425: Based on the state prior estimate and the state observation value, the observation residual is determined, the inhibition factor prediction value of the drug resistance gene sequence is taken as the state update object, and the Kalman filtering cycle iteration is performed with the PCR amplification cycle to update the residual, and finally the optimal inhibition strength factor corresponding to the drug resistance gene sequence is obtained.

[0133] Specifically, according to the state prior estimation obtained above, the observation value and the predicted value of the inhibitory factor of the target pathogen in the state observation value and the observation value and the predicted value of the abundance, the observation residual can be obtained for indirectly evaluating the inhibitory factor prediction value of the drug-resistant gene sequence, and the inhibitory factor prediction value of the drug-resistant gene sequence is taken as the state value, that is, the state update object, and the two are fused by using the Kalman gain to correct the residual error, if the prediction uncertainty is large, the observation data is more trusted; if the observation noise is large, the prediction result is more retained, and finally the optimal value of the inhibitory factor of the drug-resistant gene sequence in the current PCR amplification cycle is updated, and then the real abundance can be inferred, and this process completes one compensation cycle of “inhibitory interference stripping-real abundance restoration”, and the state estimation accuracy is gradually improved with the iteration of PCR amplification. It should be understood that, due to the binding relationship between the pathogen and the drug-resistant gene, the essence of the observation residual is the indirect projection of the state error, which does not directly modify the state, but through the three steps of “bias decoding-weight distribution-incremental correction”, the transmission from the observation bias to the state update is completed, that is, the update of the inhibitory factor of the drug-resistant gene. PCR amplification usually performs 35-40 cycles, and the filtering iteration ends with the termination of amplification, and finally the optimal inhibitory strength factor of the drug-resistant gene sequence is output.

[0134] Based on the above technical solution, by using the inhibitory factor prediction value and the abundance prediction value of the target pathogen and the drug-resistant gene sequence, the state prior estimation of different PCR amplification cycles is constructed, and by using the inhibitory factor observation value and the abundance observation value of the target pathogen and the drug-resistant gene sequence, the state observation value of different PCR amplification cycles is constructed, so as to realize the observation update process by using Kalman filtering, and finally realize the accurate estimation of the state and obtain the optimal inhibitory strength factor of the drug-resistant gene sequence.

[0135] In a possible implementation, step S500 corrects the overlap rate by using the optimal inhibitory strength factor, including: determining the ratio of the overlap rate to the optimal inhibitory strength factor; and determining the product of the ratio and the probe capture efficiency coefficient as the overlap rate correction value.

[0136] Specifically, the optimal inhibitory strength factor of the drug-resistant gene sequence is used to correct the overlap rate between the pathogen gene sequence and the drug-resistant gene sequence, and the overlap rate correction value is determined by the following formula:

[0137]

[0138] In the formula: represents the overlap rate correction value between the pathogen gene sequence and the drug-resistant gene sequence; represents the overlap rate between the pathogen gene sequence and the drug-resistant gene sequence; represents the optimal inhibitory strength factor of the drug-resistant gene sequence; represents the probe capture efficiency coefficient, which can be determined by standard sequencing data calibration, and the value is the ratio of the known concentration of the standard to the detected concentration of the sequencing). It should be understood that, since the drug resistance gene will be inhibited to a certain extent by the inhibitor, the value of is usually not 0, when the value of is 0, it means that the drug resistance gene will not be inhibited by the inhibitor, at this time, the overlap rate correction value is obtained directly by the formula .

[0139] Based on the above technical scheme, by taking the optimal inhibition intensity factor of the drug resistance gene sequence as the core, and combining the probe binding efficiency characteristics of targeted sequencing, a correction model is constructed, so that the double interference of the inhibitor and the probe preference can be eliminated, and the real amplification level underestimated by the inhibitor can be effectively restored.

[0140] In a possible implementation, as shown in Figure 6 step S500, based on the overlap rate correction value, the drug resistance gene sequence and its gene abundance in all pathogen gene sequences are determined, including:

[0141] Step S501: removing the pathogen gene sequences with an overlap rate correction value less than or equal to a set overlap rate threshold value, to obtain the remaining pathogen gene sequences after removal.

[0142] Specifically, based on the sequence alignment result, the reads that completely match the target gene region are retained, that is, the pathogen gene sequences with an overlap rate correction value greater than 95% are retained, and the interference sequences homologous to the human genome or non-target pathogen are removed, thereby obtaining the remaining pathogen gene sequences after removal.

[0143] Step S502: removing duplicate sequences from the remaining pathogen gene sequences after removal, and merging the reads from the same original deoxyribonucleic acid template to form a UMI family, calculating the sequence consistency of the pathogen gene sequences in the UMI family, and retaining the UMI family with a sequence consistency greater than a set consistency threshold value.

[0144] Specifically, for the remaining pathogen gene sequences after removal, duplicate sequences are removed by molecular tags (UMI), and reads from the same original deoxyribonucleic acid template are merged to form a UMI family. Then, the proportion of identical bases of all pathogen gene sequences in each UMI family is calculated to obtain the sequence consistency of each UMI family, and the family with a sequence consistency greater than 95% is retained to exclude amplification errors.

[0145] Step S503: removing the pathogen gene sequences containing chimeric sequences in the retained UMI family, and finally obtaining the drug resistance gene sequence and its gene abundance in all pathogen gene sequences.

[0146] Specifically, in the reserved UMI family, the read alignment breakpoint is detected, the read containing the chimeric sequence (transgene fragment fusion) is removed, the verified existing drug-resistant gene is obtained, the gene abundance of these "verified identity" genes is recalculated, and finally the detection result is output.

[0147] Based on the above technical solution, by sequentially performing small overlap rate removal, repeat sequence deduplication, small sequence consistency removal and chimeric sequence containing removal and other noise removal operations on the pathogen gene sequence, the drug-resistant gene sequence and its gene abundance in all pathogen gene sequences are finally accurately determined, and the detection accuracy is effectively improved.

[0148] Based on the same inventive concept, the embodiments of the present application also provide a pathogen-drug gene synchronous detection system for targeted sequencing, as shown in Figure 7 The system includes a memory, a processor, and computer program code stored in the memory and running on the processor, wherein when the processor executes the computer program code, the system can execute any of the above-described pathogen-drug gene synchronous detection methods for targeted sequencing.

[0149] The embodiments of the present application can divide the functional modules of the system according to the above method examples, for example, each functional module can be corresponding, or two or more functions can be integrated in one processing module, and the above integrated module can be realized in the form of hardware. It should be noted that the division of modules in this embodiment is illustrative, and is only a logical function division, and actual implementation can have another division method.

[0150] It should be noted that: the above-described embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application, and should be included in the protection scope of the present application.

Claims

1. A method for simultaneous detection of pathogen-drug resistance genes by targeted sequencing, characterized in that, The method comprises the following steps: constructing a complex knowledge graph about pathogens and drug-resistant genes based on genomic sequence information and drug-resistant gene sequence information of different pathogens; identifying a target pathogen present in a to-be-detected clinical sample and a pathogen gene sequence of the target pathogen present in the to-be-detected clinical sample based on the genomic sequence information of the different pathogens; determining a set of drug-resistant gene sequences that the target pathogen may have based on the complex knowledge graph, and comparing the pathogen gene sequence with the drug-resistant gene sequences in the set to obtain an overlap rate between the pathogen gene sequence and the drug-resistant gene sequences; analyzing sequence inhibition sensitivity of the pathogen gene sequence and the drug-resistant gene sequence of the target pathogen, and constructing a collaborative state evolution model in combination with the complex knowledge graph, wherein the collaborative state evolution model characterizes a degree of PCR amplification inhibition by introducing an inhibition intensity factor as a hidden variable, and determines an optimal inhibition intensity factor of the drug-resistant gene sequence by iterative inverse estimation based on state prior estimation determined by the collaborative state evolution model and in combination with state observation values for residual update; correcting the overlap rate by using the optimal inhibition intensity factor, and determining drug-resistant gene sequences and gene abundances in all the pathogen gene sequences based on an overlap rate correction value.

2. The method according to claim 1, wherein, The method for constructing a complex knowledge graph about pathogens and drug-resistant genes comprises: regarding each pathogen and each drug-resistant gene sequence as a node in the complex knowledge graph; based on genomic sequence information and drug-resistant gene sequence information of different pathogens, when a drug-resistant gene sequence exists in a genomic sequence of a pathogen, constructing a directed edge between nodes corresponding to the pathogen and the drug-resistant gene sequence, and determining an edge weight of the directed edge, wherein the directed edge is directed from the drug-resistant gene sequence to the pathogen, and the edge weight is used to reflect closeness of evolutionary association between the drug-resistant gene sequence and the pathogen.

3. The method according to claim 2, wherein, Determining the edge weight of the directed edge comprises: determining a copy number of the drug-resistant gene sequence in the genomic sequence of the pathogen for the drug-resistant gene sequence and the pathogen corresponding to the directed edge; determining an insertion position and a detection frequency of the drug-resistant gene sequence in different strain genomic sequences of the pathogen for the drug-resistant gene sequence and the pathogen corresponding to the directed edge; determining the edge weight corresponding to the directed edge based on the copy number, stability of the insertion position, and the detection frequency.

4. The method according to claim 1, wherein, Analyzing sequence inhibition sensitivity of the pathogen gene sequence and the drug-resistant gene sequence of the target pathogen comprises: determining a core gene sequence in the pathogen gene sequence of the target pathogen; determining GC base content and secondary structure minimum free energy of gene sequences for the core gene sequence and the drug-resistant gene sequence; Determine the inhibition susceptibility index of the target pathogen based on the GC base content and the minimum free energy of the secondary structure of all the core gene sequences, and determine the inhibition susceptibility index of the drug-resistant gene sequence based on the GC base content and the minimum free energy of the secondary structure of the drug-resistant gene sequence, which reflects the ease of inhibition of the target pathogen or the drug-resistant gene sequence by an inhibitor.

5. The method according to claim 4, wherein the pathogen-drug resistance gene is selected from the group consisting of the genes listed in Table 1. Construct a synergistic state evolution model, including: Based on the inhibition susceptibility index, correct the baseline inhibition coefficient of the to-be-detected clinical sample to obtain the initial inhibition intensity factor of the target pathogen and the drug-resistant gene sequence; Based on the initial inhibition intensity factor, construct the basic update equation of the corresponding inhibition intensity factor of the target pathogen and the drug-resistant gene sequence in the PCR amplification process; Based on the complex knowledge graph, determine each candidate pathogen associated with the drug-resistant gene sequence present in the to-be-detected clinical sample, and determine the synergistic driving term in the PCR amplification process based on the difference in inhibition intensity factor between the drug-resistant gene sequence and each candidate pathogen in the PCR amplification process. Determine the change rate of the inhibition intensity factor of the drug-resistant gene sequence in the PCR amplification process based on the change rate of the basic update equation of the inhibition intensity factor and the synergistic driving term.

6. The method according to claim 5, wherein the pathogen-drug resistance gene is selected from the group consisting of the genes listed in Table 1. Determine the synergistic driving term in the PCR amplification process, including: Determine the difference value of the inhibition intensity factor of the candidate pathogen and the drug-resistant gene sequence in the PCR amplification process to obtain the inhibition intensity factor difference value; Use the edge weight of the directed edge between the nodes corresponding to the candidate pathogen and the drug-resistant gene sequence in the complex knowledge graph to weight and accumulate the inhibition intensity factor difference value, thereby obtaining the synergistic driving term in the PCR amplification process.

7. The method according to claim 5, wherein the pathogen-drug resistance gene is selected from the group consisting of the genes listed in Table 1. Using Kalman filtering algorithm, based on the state priori estimation determined by the synergistic state evolution model, and combining with the state observation value to update the residual error, through iterative inverse estimation to determine the optimal inhibition intensity factor of the drug-resistant gene sequence, including: Obtain the inhibition factor observation value and abundance observation value of the target pathogen corresponding to each PCR amplification cycle in the Kalman filtering process, and form the state observation value corresponding to each PCR amplification cycle in the Kalman filtering process; Based on the initial inhibition intensity factor, the basic update equation of the inhibition intensity factor, and the change rate of the inhibition intensity factor, determine the inhibition factor prediction value of the target pathogen and the drug-resistant gene sequence corresponding to each PCR amplification cycle in the Kalman filtering process; Based on the inhibition factor prediction value of the target pathogen, correct the PCR exponential growth model of the target pathogen to determine the abundance prediction value of the target pathogen corresponding to each PCR amplification cycle in the Kalman filtering process; Based on the inhibition factor prediction value and the abundance prediction value of the target pathogen, determine the state priori estimation corresponding to each PCR amplification cycle in the Kalman filtering process; Based on the state prior estimation and the state observation value, an observation residual is determined, and a predicted value of an inhibition factor of the drug resistance gene sequence is taken as a state update object, and a Kalman filtering cycle iteration is performed along with a PCR amplification cycle to update the residual, and finally an optimal inhibition intensity factor corresponding to the drug resistance gene sequence is obtained.

8. The method according to claim 1, wherein, Using the optimal inhibition intensity factor, the overlap rate is corrected, including: determining a ratio of the overlap rate to the optimal inhibition intensity factor; determining a product of the ratio and a probe capture efficiency coefficient as an overlap rate correction value.

9. The method according to claim 1, wherein the method is characterized by, Based on the overlap rate correction value, the drug resistance gene sequence and its gene abundance in all the pathogen gene sequences are determined, including: eliminating the pathogen gene sequence with an overlap rate correction value less than or equal to a set overlap rate threshold value to obtain a remaining pathogen gene sequence after elimination; performing repeat sequence deduplication on the remaining pathogen gene sequence after elimination, and merging reads from the same original deoxyribonucleic acid template to form a UMI family, calculating sequence consistency of the pathogen gene sequence in the UMI family, and retaining the UMI family with a sequence consistency greater than a set consistency threshold value; eliminating the pathogen gene sequence containing a chimeric sequence in the retained UMI family, and finally obtaining the drug resistance gene sequence and its gene abundance in all the pathogen gene sequences.

10. A simultaneous detection system of pathogen-drug resistance gene targeted sequencing, characterized in that, The computer program product comprises a memory, a processor, and executable computer program code stored in the memory and executable on the processor, and the processor executes the computer program code to perform the pathogen-drug resistance gene synchronous detection method of targeted sequencing according to any one of claims 1-9.

Images