A Model Construction Method for Predicting Species Attribution of Drug Resistance Genes

Abstract

This application relates to the field of bioinformatics technology, specifically to a method for constructing a model based on gene sequencing read alignment for detecting and identifying drug resistance genes and predicting drug resistance phenotypes. This invention combines important characteristic genes related to drug resistance phenotypes to construct a model that enables direct comparison and detection of target pathogenic bacteria and their carried drug resistance genes based on sequencing read sequences, as well as prediction of drug resistance phenotypes.

Description

[0001] This application is a divisional application of Chinese patent application CN 202111680866.8, filed on December 30, 2021. Technical Field

[0002] This application relates to the field of bioinformatics technology, specifically to a method for identifying drug resistance genes and predicting drug resistance phenotypes based on gene sequencing read alignment detection. Technical Background

[0003] Accurate detection of infectious pathogens and their drug resistance is crucial for guiding precise clinical treatment. Currently, laboratory testing for drug resistance in infectious pathogens is divided into phenotypic and genotypic detection. Regarding phenotypic detection methods, the gold standard clinically is microbial culture combined with drug susceptibility testing. This method provides strong evidence for the diagnosis and treatment of clinical infections, but it also has limitations, such as time-consuming culture (generally 2-4 days), low positive rate of pathogen culture, and susceptibility to various uncertain factors during the culture process. Some difficult-to-culture or rare pathogens may even fail to be cultured successfully. Other detection techniques, such as the Carba NP test, modified carbapenem inactivation assays (including mCIM and eCIM), carbapenemase inhibitor enhancement assays, and time-of-flight mass spectrometry, are for detecting carbapenemases in common Enterobacteriaceae. They are only suitable for detecting carbapenem resistance caused by the production of related enzymes and cannot detect resistance caused by other mechanisms (such as efflux pump mechanisms). Genotypic detection includes enzyme immunochromatography and gene detection techniques. Enzyme immunochromatography and conventional gene detection technologies (such as GeneXpert, Verigene, and Filmarray detection systems) target specific genes and are characterized by rapid and easy interpretation. However, if the gene being tested is different from the target gene, false negative results are likely to occur.

[0004] With the continuous development of sequencing technology and the reduction of sequencing costs, sequencing-based detection of pathogens and their drug resistance genes has gradually become popular. Microbial genome sequencing includes two strategies: whole-genome sequencing (WGS) and metagenomic sequencing (mNGS). For whole-genome sequencing of bacterial strains, studies have shown that models (rule-based models or machine-learning models) can be constructed based on the genome data of all strains of a single population and their corresponding drug resistance phenotypes. These models are then used to analyze the genome data of new strains to predict their drug resistance phenotypes, achieving good results (accuracy rates exceeding 95%). However, a limitation exists: whole-genome sequencing, like other methods, cannot bypass the restriction of pathogen culture enrichment. Metagenomic high-throughput sequencing (mNGS), a novel pathogen detection method developed in the last decade, differs from microbial culture in that it does not require screening to obtain pure cultures of all pathogens in the environment. It can read the base sequences of all pathogen nucleic acids from a small number of samples at once, providing information such as the type of pathogen. mNGS (mundular next-generation sequencing) is characterized by rapid detection, broad pathogen coverage, and unbiasedness. Numerous applied research articles have been published in high-level SCI journals. Furthermore, many of these published articles analyze and discuss the detection of drug resistance genes, demonstrating the potential value of mNGS in predicting the drug susceptibility of pathogens. However, for clinical specimens (bronchoalveolar lavage fluid, blood, cerebrospinal fluid, etc.), we know that they are often heavily contaminated by the host (host nucleic acid accounts for more than 95% of the total nucleic acid extracted from the sample), resulting in low microbial loads. With current conventional sequencing volumes of around 20M reads, the bacterial genome coverage obtained from clinical specimens often does not exceed 1X. Therefore, given this amount of sequencing data, it is necessary to further explore the detection efficacy of mNGS in identifying pathogens and simultaneously accurately detecting drug resistance genes for predicting pathogen resistance. Summary of the Invention

[0005] To address the aforementioned technical challenges, this invention combines key characteristic genes associated with drug resistance phenotypes identified through BGWAS screening to construct a data analysis method and system for directly comparing, detecting, identifying, and predicting drug resistance phenotypes of target pathogenic bacteria and their carried drug resistance genes based on NGS or nanopore sequencing read sequences. For the genomes of all strains included in the training set of the previous BGWAS study, target drug resistance genes were detected by simulating NGS or nonopore sequencing read sequence alignment (read-based) and genome contig sequence alignment (assembly-based). The results of the assembly-based method were used as a reference to validate and optimize the read-based detection process, aiming to achieve accurate genotyping through read-based detection. Then, a custom formula was used to calculate the score, which was used as an indicator to predict antibiotic drug susceptibility. ROC analysis was performed using read sequence simulation tests to determine the optimal cutoff threshold, while simultaneously evaluating the accuracy and performance of the prediction model. Finally, the effectiveness of the analysis system was evaluated using clinical specimens or cultured, pure pathogenic bacterial strains.

[0006] Specifically, this application proposes the following technical solution:

[0007] This application first provides a method for detecting and identifying drug resistance genes and constructing a model for predicting drug resistance phenotypes based on gene sequencing read alignment, comprising the following steps:

[0008] Step 1): Combine, for example, the classification information of drug-resistant genes in the CARD drug resistance database, organize and correct the drug resistance gene family information and the consistency annotation information between gene sequences. Preferably, it also includes organizing and correcting the source species information and / or the mediation mode.

[0009] Step 2): Calculate the weight coefficient of the drug resistance gene family. Based on the weight coefficient of the genes within the family and the sample detection frequency of the corresponding genes in the BGWAS model training set, calculate the gene family weight coefficient.

[0010] Step 3): Detection of drug resistance genes and process correction.

[0011] Furthermore, in step 1),

[0012] The gene families are defined based on the gene family information recorded in the drug resistance database, and are also reviewed and corrected by referring to the drug resistance gene family information recorded in the NCBINDARO database and MEGARes database.

[0013] The source species information is obtained by comparing all drug-resistant reference genes with the NCBI NT database, retaining hits with identity >= 95% and subject coverage >= 95%, to obtain all species annotation information for each drug-resistant gene;

[0014] The mediation method involves querying the reference sequence description information in the alignment to determine whether the drug resistance gene is mediated by a plasmid.

[0015] The gene sequence consistency annotation information is obtained by comparing all drug resistance gene sequences pairwise to obtain the consistency value among all gene sequences.

[0016] Furthermore, in step 2), the weighting coefficient is defined by the following formula:

[0017]

[0018] In the formula, arg_Ni is the number of samples in the BGWAS model training set where the corresponding gene in the target gene family is detected, arg_Wi is the weight coefficient of the corresponding gene in the family, j indicates that there are j key genes in the family, and j+k indicates the total number of genes in the family.

[0019] Furthermore, step 3) of detecting and correcting drug resistance genes includes the following steps:

[0020] a) Simulation of sequencing reads data;

[0021] b) Sequence alignment and annotation statistics;

[0022] c) Screening and filtering of drug resistance genes.

[0023] Furthermore, the a) sequencing reads data simulation is based on strain samples from the BGWAS model training set; preferably:

[0024] Short reads from NGS sequencing of bacterial strain genomes were simulated using ART_Illumina software.

[0025] The ReadSim software was used to simulate Nanopore sequencing reads of bacterial strain genomes.

[0026] More preferably, the simulation simulates gradient data volumes of 0.05X, 0.1X, 0.2X, 0.3X, 0.4X, 0.5X, 0.6X, 0.7X, 0.8X, 0.9X, 1X, 2X, 3X, 5X, 10X, and 30X.

[0027] Furthermore, the sequence alignment and annotation statistics in section b) include:

[0028] The simulated read sequences were compared with the drug resistance gene library to filter out low-quality hits; final gene annotation was performed, and the number of specific reads, multiple alignment reads, and reads belonging to the drug resistance gene family were counted in the sample to calculate the coverage of the detected genes.

[0029] Preferably, the best hit is selected as the final gene annotation of each read sequence using the best alignment and LCA algorithm. If there are multiple hits with the same value, i.e., multiple alignments, the LCA algorithm is used to annotate the read sequence. That is, for a single read sequence, if the genotype cannot be annotated due to multiple alignments, the annotation is moved to a higher level to the gene family level.

[0030] More preferably,

[0031] For NGS sequencing data, the blastn software was used to compare the simulated read sequences with the drug resistance gene library, filtering out and retaining hits with an identity greater than 90%.

[0032] For Nanopore sequencing data, minimap2 software was used to compare with a drug resistance reference gene library, filtering out hits with an identity below 0.7 or a subject coverage below 0.4.

[0033] Furthermore, the screening and filtering of drug resistance genes in step c) is as follows: screening and filtering drug resistance genes that have been aligned with the reads sequence in step b).

[0034] Preferably, the filtering criteria include one or more of the following:

[0035] A) Evaluate the impact of sequence consistency among different genotypes in the drug resistance gene reference library on read-based accurate detection of drug resistance gene genotyping: By selecting drug resistance genes with different maximum sequence consistency in the database, short read sequences were simulated for read-based detection and analysis. The number of specific reads detected for the target gene and the number of reads that matched the target gene were counted. For the strategy of identifying genotypes by detecting specific read sequences, 95% identity was selected as the threshold standard for whether the target gene can be easily accurately genotyped.

[0036] B) Screening based on drug resistance genotyping: For target genes with high similarity to other genes in the database, the gene with the highest number of reads that can be compared to the target gene is considered a true positive. For those not ranked first, genes with 100% coverage calculated based on precise read comparisons are retained as true positives. For target genes with low similarity to other genes in the database, the number of specific reads detected is used to determine whether the result is a true positive. Gene families for which no specific reads were detected are directly filtered out.

[0037] C) Evaluate the accuracy of read-based detection for identifying drug resistance gene genotypes: Using the assembly-based detection results of drug resistance genes during the BGWAS model training process as a reference, statistically analyze the accuracy, sensitivity, and specificity of read-based detection of drug resistance genes or gene families for important genes or gene families screened by the BGWAS model.

[0038] Further, step 4): Define and calculate the Score value for positive and negative judgment indicators, and determine the reporting rules and cutoff threshold based on ROC analysis;

[0039] The score value is calculated as follows:

[0040]

[0041] In the formula, arg_Wi represents the weight coefficient of the corresponding genotype, and genefamily_Wi represents the weight coefficient of the corresponding gene family. When a genotype is detected and the genotype weight coefficient is >0, the genotype weight coefficient is used for calculation. When a genotype is detected but the genotype weight coefficient is 0 or there is no weight coefficient, the gene family weight coefficient is used for calculation.

[0042] Furthermore, the sequencing reads are first-generation, second-generation, or third-generation sequencing reads, preferably NGS or Nanopore sequencing reads; more preferably NGS or Nanopore metagenomic sequencing reads.

[0043] This invention also provides a model construction method for predicting the species attribution of drug resistance genes, characterized in that the method includes the following steps:

[0044] Step 1): Alignment of the genome sequence of the target pathogen and calculation of the number of detected sequences, genome coverage, and coverage depth;

[0045] Step 2): Based on the detection results of drug resistance genes in the BGWAS model training set specimens, count the copy number of drug resistance genes carried by the target pathogen species;

[0046] Step 3): Based on the assumed gene-species attribution relationship, calculate the copy number of the drug resistance gene and determine the species attribution.

[0047] Furthermore, step 1) specifically involves:

[0048] Commonly known clinical pathogens were selected as target pathogens. The reference genome of the target pathogen was searched and downloaded from the NCBI genome database and used as a reference sequence library for the identification of the target pathogen species.

[0049] Each sequenced read was compared with the above reference sequence library, and the number of detected sequences, genome coverage, and coverage depth of the target pathogen species were calculated. The total number of detected pathogen sequences, genome coverage, and coverage depth were then statistically obtained.

[0050] Furthermore, step 2) specifically involves:

[0051] Based on the assembly-based drug resistance gene detection results of the training set samples during BGWAS model training, the detection distribution and copy change range of drug resistance genes and drug resistance gene families of the target pathogen species were statistically obtained.

[0052] Furthermore, step 3) specifically involves:

[0053] When assuming a gene-species correspondence for drug resistance, the main criteria are: a) whether the species annotation of the reference gene in the database includes the target species; if so, the hypothesis of the gene-species attribution is accepted. b) if a is not met, check the mediation mode annotation of the reference gene in the database to see if it includes plasmid-mediated mediation; if so, the hypothesis of the gene-species attribution is accepted. c) if a and b are not met, infer the species origin based on the species annotation of ARG-like reads, and calculate the copy number of the drug resistance gene using the following formula:

[0054]

[0055] If the calculated copy number of the drug resistance gene falls within the normal copy number range of the target gene family obtained from the BGWAS model training set, then the assumed gene-species attribution relationship is accepted; otherwise, it is rejected.

[0056] This invention also provides a method for detecting drug resistance in metagenomic sequencing data, comprising the following steps:

[0057] 1) Perform quality control and remove human nucleic acid sequences from the sample sequencing data;

[0058] 2) Detect and identify the drug-resistant genes contained in the sample: Based on the above detection and identification method, conduct drug-resistant gene alignment and annotation statistics on the sample sequence to detect and identify the drug-resistant genes contained in the sample;

[0059] 3) Predict the species attribution of the detected drug-resistant genes in the sample: Identify the target pathogenic bacteria contained in the sample according to the above species attribution prediction method, and predict the species attribution of the detected drug-resistant genes;

[0060] 4) For the target pathogenic bacteria, according to the detected drug-resistant gene carriage situation, calculate the score value of the target species - antibiotic drug according to the above score calculation method, and compare it with the cutoff value: When score >= cutoff, it is predicted as R; when score < cutoff, if the detected pathogen genome coverage is higher than the minimum genome coverage or data volume required for model stability, it is predicted as S, otherwise it is reported as unknown.

[0061] The present invention also provides a model for detecting and identifying drug-resistant genes based on gene sequencing reads alignment, including the following modules:

[0062] Module 1): Used to combine the classification information of drug-resistant genes in the CARD drug-resistant database, sort and correct the information of drug-resistant gene families and the consistency annotation information between gene sequences. Preferably, it also includes sorting and correcting the source species information and / or the mediation method;

[0063] Module 2): Used to calculate the weight coefficient of the drug-resistant gene family, and calculate the weight coefficient of the gene family based on the weight coefficient of the member genes within the family and the sample detection frequency of the corresponding genes in the BGWAS model training set;

[0064] Module 3): Used for the detection of drug-resistant genes and process correction.

[0065] The present invention also provides a prediction model for the species attribution of drug-resistant genes. The method modules are as follows:

[0066] Module 1): Used for the alignment of the genome sequence of the target pathogenic bacteria species and the calculation of the detected sequence number, genome coverage and coverage depth;

[0067] Module 2): Used to count the copy number of drug-resistant genes carried by the target pathogenic bacteria species based on the drug-resistant gene detection results of the specimens in the BGWAS model training set;

[0068] Module 3): Used for the calculation of the copy number of drug-resistant genes and the judgment of species attribution based on the assumed gene - species attribution relationship.

[0069] The further limitations of each module in the above model are the same as the limitations of each step in any of the above methods.

[0070] The present invention also provides an apparatus comprising: at least one memory for storing a program; and at least one processor for loading the program to perform the method as described in any of the preceding claims.

[0071] The present invention also provides a storage medium storing processor-executable instructions, which, when executed by a processor, are used to implement the method as described in any of the preceding claims.

[0072] This invention also provides the following:

[0073] Application of genes AAC(3)-IIe, AAC(3)-IV, AAC(3)-IId, rmtC, armA, rmtF, rmtB, AAC(6')-33 and ANT(2”)-Ia as non-core drug resistance genes in the auxiliary drug susceptibility prediction of Klebsiella pneumoniae;

[0074] The drug sensitivity prediction includes drug resistance prediction and sensitivity prediction, preferably sensitivity prediction;

[0075] More preferably, the drug sensitivity test is specific to gentamicin.

[0076] Application of detection reagents for non-core drug resistance genes AAC(3)-IIe, AAC(3)-IV, AAC(3)-IId, rmtC, armA, rmtF, rmtB, AAC(6')-33 and ANT(2”)-Ia in the preparation of Klebsiella pneumoniae auxiliary drug susceptibility prediction kit;

[0077] The drug sensitivity prediction includes drug resistance prediction and sensitivity prediction, preferably sensitivity prediction;

[0078] More preferably, the drug sensitivity test is specific to gentamicin.

[0079] Further preferred methods involve simultaneously detecting the AAC(3)-IIe, AAC(3)-IV, AAC(3)-IId, rmtC, armA, rmtF, rmtB, AAC(6')-33, and ANT(2”)-Ia genes, which are the main mediators of drug resistance and occur frequently with high weights. If all test results are negative, the drug is considered sensitive.

[0080] Application of genes AAC(3)-IV, AAC(3)-IId, AAC(6')-Ib', AAC(6')-Ib-cr, AAC(6')-Ib-Hangzhou, AAC(6')-Ib4, mphE, ANT(2”)-Ia and aadA24 as non-core drug resistance genes in the auxiliary drug susceptibility prediction of Klebsiella pneumoniae;

[0081] The drug sensitivity prediction includes drug resistance prediction and sensitivity prediction, preferably sensitivity prediction;

[0082] More preferably, the drug sensitivity test is specific to tobramycin.

[0083] Application of detection reagents for non-core drug resistance genes AAC(3)-IV, AAC(3)-IId, AAC(6')-Ib', AAC(6')-Ib-cr, AAC(6')-Ib-Hangzhou, AAC(6')-Ib4, mphE, ANT(2”)-Ia and aadA24 in the preparation of Klebsiella pneumoniae auxiliary drug susceptibility prediction kit;

[0084] The drug sensitivity prediction includes drug resistance prediction and sensitivity prediction, with sensitivity prediction being preferred.

[0085] More preferably, the drug sensitivity test is specific to tobramycin.

[0086] Further preferred methods involve simultaneously detecting the AAAC(3)-IV, AAC(3)-IId, AAC(6')-Ib', AAC(6')-Ib-cr, AAC(6')-Ib-Hangzhou, AAC(6')-Ib4, mphE, ANT(2”)-Ia, and aadA24 genes, which are the main mediators of drug resistance and occur frequently with high weights. If all the test results are negative, the drug is considered sensitive.

[0087] Application of genes CTX-M-55, CTX-M-11, CTX-M-15, SHV-155, SHV-5, SHV-11, SHV-12, SHV-76, SHV-30, SHV-53, SHV-124, SHV-182, DHA-1, KPC-3, and KPC-2 as non-core drug resistance genes in the auxiliary drug susceptibility prediction of Klebsiella pneumoniae.

[0088] The drug sensitivity prediction includes drug resistance prediction and sensitivity prediction, preferably sensitivity prediction;

[0089] More preferably, the drug sensitivity test is for ceftazidime.

[0090] Application of detection reagents for non-core drug resistance genes CTX-M-55, CTX-M-11, CTX-M-15, SHV-155, SHV-5, SHV-11, SHV-12, SHV-76, SHV-30, SHV-53, SHV-124, SHV-182, DHA-1, KPC-3, and KPC-2 in the preparation of an auxiliary antimicrobial susceptibility prediction kit for Klebsiella pneumoniae;

[0091] The drug sensitivity prediction is a drug resistance prediction;

[0092] Preferably, the drug sensitivity test is for ceftazidime.

[0093] More preferably, by detecting the CTX-M-55, CTX-M-11, CTX-M-15, SHV-155, SHV-5, SHV-11, SHV-12, SHV-76, SHV-30, SHV-53, SHV-124, SHV-182, DHA-1, KPC-3, and KPC-2 genes, which are the main mediators of drug resistance and occur frequently with high weights, drug resistance is presumed if all the test results are positive.

[0094] Application of genes dfrA12, dfrA15, dfrA17, dfrA19, dfrA30, dfrA8, dfrA5, dfrA15b, dfrA14, dfr22, dfrA27 and dfrA1 as non-core drug resistance genes in the auxiliary drug susceptibility prediction of Klebsiella pneumoniae;

[0095] The drug sensitivity prediction is a drug resistance prediction;

[0096] Preferably, the drug sensitivity test is for the compound sulfamethoxazole drug.

[0097] Application of detection reagents for non-core drug resistance genes dfrA12, dfrA15, dfrA17, dfrA19, dfrA30, dfrA8, dfrA5, dfrA15b, dfrA14, dfr22, dfrA2 and dfrA1 in the preparation of an auxiliary antimicrobial susceptibility prediction kit for Klebsiella pneumoniae;

[0098] The drug sensitivity prediction is a drug resistance prediction;

[0099] Preferably, the drug sensitivity test is for the compound sulfamethoxazole drug.

[0100] More preferably, by simultaneously detecting the dfrA12, dfrA15, dfrA17, dfrA19, dfrA30, dfrA8, dfrA5, dfrA15b, dfrA14, dfr22, dfrA27 and dfrA1 genes that mainly mediate drug resistance and occur frequently with high weights, if the test results are all positive, drug resistance is presumed.

[0101] The beneficial technical effects of this application are:

[0102] 1) This invention is based on a nucleic acid molecular detection method for drug resistance, which bypasses the limitations of traditional culture and directly performs metagenomic sequencing on clinical specimens to identify target pathogens and their drug resistance gene carrying status. Furthermore, it predicts the drug sensitivity results of antibiotics based on the presence or absence of drug resistance genes. This invention is also applicable to pure strain specimens.

[0103] 2) This invention directly aligns and detects drug resistance genes based on NGS or non-Apore sequencing reads. Compared to genome contig-based alignment, this bypasses the genome assembly step and has higher detection sensitivity. Specifically, a read-based drug resistance gene alignment and detection method and a corresponding database are constructed, and the results of asembly-based drug resistance gene detection are used as a reference to compare read-...

[0104] The performance of the read-based process was validated and evaluated to ensure the accuracy of read-based detection of drug resistance genes.

[0105] 3) This invention addresses read-based drug resistance gene detection by constructing a specific reference database for drug resistance genes, particularly focusing on the organization of drug resistance gene genotyping and classification, enabling the use of an LCA (lowest level of autologous) annotation strategy for query sequences. Specifically, it includes all gene sequences from the CARD drug resistance public library as reference genes. Referring to the multi-level gene annotation method of the MEGARes database and the family information of drug resistance genes recorded in the NCBI NDARO database, each reference gene is annotated at six levels, thus enabling LCA (lowest level of autologous ...

[0106] The common ancestors annotation strategy, i.e., obtaining the number of detection-specific reads at each level. Using OXA-181...

[0107] For example, the labeling for a gene at level 6 is: OXA-181__1(L1_geneST),OXA-

[0108] 181(L2_genetype),OXA-48subfamily(L3_subgroup),OXA family(L4_Group),

[0109] Class_D_betalactamases(L5_Mechanism),betalactams(L6_Class).

[0110] 4) This invention employs a combined strategy of two rules for the precise detection of read-based drug resistance genotyping, effectively improving the genotyping capability and accuracy of the detection process. First, for target genes with high similarity to other genes in the database (e.g., consistency exceeding 95%), the gene with the highest ARG-like read count (i.e., the total number of reads that can be compared to the target gene, or the number of specific reads plus the number of multiple alignment reads) is considered a true positive. For genes not ranked first, those with 100% coverage calculated based on precise alignment reads are retained as true positives. Second, for target genes with low similarity to other genes in the database (consistency below 95%), the determination of a true positive result is primarily based on the number of detected specific reads. For some genes without detected specific reads...

[0111] The family of reads will be considered false positives and filtered out directly.

[0112] 5) This invention defines a method for calculating the weight coefficient of the corresponding drug resistance gene family based on the weight coefficient of drug resistance gene typing.

[0113] (or formula) In cases where genotyping may be inaccurate, it replaces the use of gene family weights in calculating and predicting drug sensitivity results, effectively avoiding false positives that may be caused by inaccurate genotyping.

[0114] 6) This invention defines a method for predicting drug sensitivity results, namely, defining a formula for calculating the positive / negative interpretation index Score, and simultaneously combining gradient simulation tests with different sequencing data volumes to evaluate the performance of the prediction model and determine the cutoff.

[0115] The threshold is ultimately used to effectively predict antibiotic susceptibility results.

[0116] 7) This invention addresses metagenomic sequencing (mixed microbial community sequencing) of clinical specimens, defining a method (or formula) for calculating the copy number of drug-resistant genes based on the detected sequences of target pathogenic bacteria and drug-resistant genes. Based on the calculated copy number, the possible pathogenic species origin of the drug-resistant gene is predicted and assessed. Specifically, when inferring the correspondence between drug-resistant genes and pathogenic species, the drug-resistant gene is first assumed to originate from a target species based on the actual detected drug-resistant gene and pathogenic species information. Then, the copy number of the drug-resistant gene is calculated, and it is checked whether the calculated copy number falls within the normal range. If it is normal, the attribution is considered acceptable; otherwise, it is rejected. When assuming a gene-species correspondence for drug resistance, the main criteria are: a) whether the species annotation of the reference gene in the database includes the target species; if so, the hypothesis of gene-species attribution is accepted. b) if a is not met, the mediation information of the reference gene annotated in the database is checked to see if plasmid-mediated mediation is included; if so, the hypothesis of gene-species attribution is accepted. c) if a and b are not met, the species origin is inferred based on the species annotation of ARG-like reads.

[0117] 8) Taking the detection of antibiotic resistance in Klebsiella pneumoniae as an example, this invention can accurately identify Klebsiella pneumoniae species and the drug resistance genes they carry, and can effectively predict the resistance of carbapenems (imipenem, meropenem) and aminoglycosides.

[0118] The drug susceptibility results for gentamicin and tobramycin, as well as the prediction of resistance to ceftazidime and trimethoprim-sulfamethoxazole, show an accuracy rate exceeding 90%. Clinical specimen sampling verification shows that the accuracy rate for predicting carbapenem susceptibility can reach 100%, with over 80% of samples providing clear susceptibility predictions. This invention can assist in the clinical detection and diagnosis of drug-resistant bacteria. Attached Figure Description

[0119] Figure 1 Technical roadmap of the present invention;

[0120] Figure 2 The figure shows the results of a drug resistance gene detection test using simulated 100X target gene read data for target genes with different identities from other genes in the database. In the figure, ARG-like represents the proportion of the number of detected target gene or non-target gene sequences to the number of detected sequences of the target gene family, and Specific represents the number of detected target gene-specific sequences.

[0121] Figure 3 A graph showing the accuracy performance of read-based genotyping or gene family typing, using assembly-based drug resistance gene detection results as a reference.

[0122] Figure 4 Technical flowchart for Score calculation and reporting rules based on read-based drug resistance gene detection Figure 5 Performance (AUC) curves of six antibiotic resistance prediction models under different sequencing data volumes Figure 6 Performance (AUC) values ​​of six antibiotic drug prediction models were simulated for training and validation sets with 30X genomic data. Detailed Implementation

[0123] The embodiments of this application will be described in detail below with reference to examples. However, those skilled in the art will understand that the following examples are for illustrative purposes only and should not be considered as limiting the scope of this application. Unless otherwise specified in the examples, conventional conditions or conditions recommended by the manufacturer shall apply. Reagents or instruments whose manufacturers are not specified are all conventional products that can be purchased on the market.

[0124] Definitions of some terms

[0125] Unless otherwise defined below, all technical and scientific terms used in the specific embodiments of this application are intended to have the same meaning as commonly understood by those skilled in the art. While it is believed that the following terms will be well understood by those skilled in the art, the following definitions are set forth to better explain this application.

[0126] As used in this application, the terms “comprising,” “including,” “having,” “containing,” or “involving” are inclusive or open-ended and do not exclude other unlisted elements or method steps. The term “consisting of” is considered a preferred embodiment of the term “comprising.” If a group is defined below as comprising at least a certain number of embodiments, this should also be understood to disclose a group that preferably consists only of those embodiments.

[0127] When referring to a singular noun, the indefinite or definite article used, such as "a" or "a kind of," "the," includes the plural form of the noun.

[0128] The term "approximately" in this application refers to an accuracy range that, as would be understood by those skilled in the art, still guarantees the technical effects of the discussed features. This term typically indicates a deviation from the indicated value of ±10%, preferably ±5%.

[0129] Furthermore, the terms first, second, third, (a), (b), (c), and similar terms used in the specification and claims are for distinguishing similar elements and are not necessary for the order of description or chronological sequence. It should be understood that such terms are interchangeable in appropriate contexts, and the embodiments described herein can be implemented in a different order than that described or illustrated herein.

[0130] The "drug resistance" mentioned in this application, also known as antibiotic resistance, refers to the tolerance of microorganisms, parasites, and tumor cells to the effects of drugs. Once drug resistance develops, the effectiveness of the drug is significantly reduced. This application preferably refers to the resistance of bacteria in vivo to antibiotics.

[0131] The “drug resistance phenotype” mentioned in this application generally refers to the drug resistance characteristics presented by an organism, which is called the drug resistance phenotype, and the drug resistance gene it possesses is called the drug resistance genotype.

[0132] The "non-core genes" described in this application refer to genes that exist only in some strains of a particular bacterial population, as opposed to core genes, which are present in all strains. The drug resistance genes detected by the method in this application are primarily targeted at these non-core genes.

[0133] The “important characteristic genes” mentioned in this application refer to the aforementioned non-core drug resistance genes, namely, drug resistance characteristics or drug resistance genes that are significantly associated with the drug resistance phenotype of a certain antibiotic.

[0134] The “read-based” mentioned in this invention refers to read sequence alignment: directly aligning the sequenced reads with a drug resistance gene library to detect and analyze drug resistance genes.

[0135] The “assembly-based” mentioned in this invention refers to genome contig sequence alignment; the sequencing reads are assembled into species genomes to obtain contigs, and then drug resistance gene libraries are aligned based on the contig sequences to detect and analyze drug resistance genes.

[0136] The “BGWAS” or “BGWAS model” mentioned in this invention refers to bacterial genome-wide association analysis, which involves analyzing the association between bacterial genome data and drug resistance phenotype data to screen for important drug resistance characteristics or drug resistance genes that are significantly associated with drug resistance phenotypes.

[0137] Correspondingly, the “BGWAS model training set” refers to the data of all bacterial strains used in conducting bacterial genome-wide association analysis, i.e., the model training set.

[0138] For details regarding "BGWAS" or "BGWAS model," please refer to the applicant's earlier patent CN202111400540.5. This model specifically includes the following modules:

[0139] Module 1) is used to acquire the genomic data of the target bacterial strain and collect the corresponding drug susceptibility test results.

[0140] Module 2) is used for alignment and annotation of antibiotic resistance databases based on contig sequences of bacterial genomes;

[0141] Module 3) is used to perform genotype and drug resistance phenotype data association analysis for the target drug, screen important characteristic genes related to drug resistance, and calculate the weight coefficient of important characteristic genes; preferably, the important characteristic genes are non-core drug resistance genes.

[0142] Module 4) ROC analysis evaluates the performance of models that predict drug susceptibility results based on screened important genes.

[0143] The ROC analysis is as follows: Based on the matrix of important gene weight coefficients obtained in step 3), the Score value is defined and calculated, and used as the indicator for positive and negative interpretations. An ROC curve is plotted, and the cut-off value is determined. The model performance is then validated and evaluated using a validation set. Among them ar / _W i This represents the weighting coefficient value for the detected gene.

[0144] Furthermore, in step 1), the number of strain genomes is >= 100, the strain sources cover various subtypes, and the ratio of drug-resistant to susceptible strains is balanced; in some preferred embodiments, the acquisition is done by searching and downloading published target genome sequences from public databases, or by sequencing and assembling bacterial strains identified through current clinical culture; in some more preferred embodiments, the search and download from public databases is as follows: collecting information on bacterial strains with drug susceptibility test results from the NCBI NDARO database and PATRIC database platform, organizing phenotypic data, and downloading genome data in batches from the NCBI genome database according to the genome assembly ID number or from the PATRIC database according to the PATRIC ID. Furthermore, the alignment annotation in step 2) involves: aligning the contig sequence with the CARD drug resistance gene reference sequence library, filtering out hits with low identity and coverage (preferably, first filtering hits with identity less than 90% or reference gene coverage less than 90%), then selecting the best hit from each contig alignment as the final alignment result for that contig region, and adding annotation information for the drug resistance gene. Furthermore, the association analysis in step 3) uses a Lasso regression model for association analysis. Furthermore, the association analysis method of the Lasso regression model in step 3) specifically involves: using the gene detection distribution matrix and the antibiotic susceptibility test result data matrix as input, performing association analysis of genotype and drug resistance phenotype data using the glmnet package, and performing k-fold (preferably k = 5-15) cross-validation to screen out important characteristic genes related to the drug resistance phenotype, and calculating the weight coefficients of the important characteristic genes; furthermore, the important characteristic genes are specifically selected based on the model CV error rate and AUC change curves under different numbers of characteristic genes, choosing the gene corresponding to the lowest CV error rate and the relatively stable model AUC value at this point as the important characteristic gene. Further, step 3) may include manual recall, which involves manually recalling genes with a high PPV (preferably PPV >= 0.8) associated with the drug resistance phenotype, and calculating the weight coefficients of the recalled genes based on the weight coefficient values ​​of the important genes obtained above. Furthermore, the bacteria described in this application include, but are not limited to, Escherichia coli, Klebsiella pneumoniae, Acinetobacter baumannii, Pseudomonas aeruginosa, Enterobacter cloacae complex, Staphylococcus aureus, Enterococcus faecalis, Enterococcus faecium, Streptococcus pneumoniae, Streptococcus pyogenes, Haemophilus influenzae, and Staphylococcus epidermidis; Klebsiella pneumoniae is preferred.Furthermore, the drug resistance phenotypes described in this application include, but are not limited to, phenotypes of resistance to carbapenems, cephalosporins, penicillins, β-lactam antibiotic inhibitors, aminoglycosides, sulfonamides, tetracyclines, quinolones, glycopeptides, oxazolidinones, and polymyxins; preferably, the drug resistance phenotype is a phenotype of resistance to carbapenems.

[0145] The present application will now be described in conjunction with specific embodiments.

[0146] Example 1: Establishment of the Method of the Invention

[0147] Figure 1 The technical roadmap of this invention is shown below, with each step described as follows:

[0148] 1) Based on the important genes screened by BGWAS, a drug resistance gene detection and phenotypic prediction process based on sequencing read alignment was constructed, and the results were tested and verified using simulated data to determine the positive cutoff value.

[0149] 1.1 Based on the classification information of drug-resistant genes in the CARD drug resistance database (V3.1.0), the annotation information, including the gene family, possible source species, mediation mode, and sequence consistency, was reorganized and corrected. Gene family definitions were based on the gene family information recorded in the CARD database, while also referencing and correcting information from the NCBI NDARO and MEGARes databases. For example, the OXA family can be further subdivided into subfamilies such as OXA-48 family and OXA-51 family. Different subfamilies may contribute differently to the development of resistance to different antibiotics; therefore, specific OXA genotypes need to be defined at the subfamily level, such as OXA-181 and OXA-232 belonging to the OXA-48 family. Secondly, all drug-resistant reference genes were aligned with the NCBI NT database, retaining hits with identity >= 95% and subject coverage >= 95% to obtain species-wide annotation information for each drug-resistant gene. The keyword "plasmid" was then searched in the reference sequence descriptions to determine if the drug-resistant gene was plasmid-mediated. Finally, BLASTN software was used to perform pairwise alignments of all drug-resistant gene sequences to obtain the consistency values ​​among all gene sequences.

[0150] 1.2 Calculation of Weight Coefficients for Drug Resistance Gene Families. Based on the weight coefficients of genes within a family and the sample detection frequencies of the corresponding genes in the BGWAS model training set (see the applicant's prior patent CN202111400540.5 for details), the weight coefficients of the gene family are calculated. The calculation formula is as follows:

[0151]

[0152] In the formula, arg_N i arg_W represents the number of samples in the training set (i.e., the samples used for model training in the previous BGWAS analysis) where the corresponding gene within the target gene family was detected. i This represents the weighting coefficient of the corresponding gene within the family. j indicates that there are j key genes in the family, and j+k represents the total number of genes in the family.

[0153] 1.3 Simulating NGS or ONT sequencing reads for drug resistance gene detection and workflow correction

[0154] 1.3.1 Reads Data Simulation

[0155] Based on the training set of bacterial strains used in the previous BGWAS model, short reads of bacterial genome sequencing (Ilumina SE75) were simulated using ART_Illumina software (Version 2.5.8) (parameter settings: -ss NS50-l75-f 5-nf 0-rs 1). Nanopore sequencing reads were simulated using ReadSim software (Version 1.6) with the following parameters: --rev_strd on-tech nanopore--read_mu 3000--read_dist normal. Considering that the detection depth of pathogenic bacterial genomes is generally no more than 1X under routine sequencing of 20M reads in clinical specimens, gradient data volumes of 0.05X, 0.1X, 0.2X, 0.3X, 0.4X, 0.5X, 0.6X, 0.7X, 0.8X, 0.9X, 1X, 2X, 3X, 5X, 10X, and 30X were simulated.

[0156] 1.3.2 Sequence Alignment and Annotation Statistics

[0157] The simulated Illumina reads were compared with the drug resistance gene database using BLASTN software (parameter settings: -evalue 1e-5-outfmt 6). First, the reads were filtered to retain only those with an identity greater than 90%. Then, the best hit was selected as the final hit for each read sequence. If multiple hits with the highest score had the same value (i.e., multiple alignments), the LCA algorithm was used to annotate the read sequence (i.e., for a single read sequence, if multiple alignments do not annotate the genotype, the annotation is moved to a higher level to the gene family level). The number of specific reads, multiple alignment reads, and family-specific reads of drug resistance genes were then counted in the sample, and the coverage index of each detected gene was calculated.

[0158] For the simulated nanopore reads sequence data, minimap2 software (version 2.17) was used for alignment with the drug resistance reference gene library. The parameters were set as follows: -cx map-ont-L--secondary=no. Hit sequences with an identity below 0.7 or subject coverage below 0.4 were then filtered out. Finally, the best alignment or LCA algorithm was used to perform gene annotation on the read sequences, and the number of specific reads, multiple alignment reads, and reads belonging to the drug resistance gene family were counted to calculate the coverage of the detected genes.

[0159] 1.3.3 Screening and filtering of drug resistance genes with aligned read sequences

[0160] A) Evaluate the impact of sequence identity among different genotyping genes in the drug resistance gene reference library on read-based accurate detection of drug resistance gene genotyping. By selecting drug resistance genes with different maximum sequence identity in the database, a read-based detection process was performed, simulating short read sequences. The number of specific reads and ARG-like reads (i.e., the total number of reads aligned to the target gene) were counted. Finally, regarding the strategy of identifying genotypes based on the detection of specific read sequences, 95% identity can be selected as the threshold standard for whether accurate genotyping of the target gene can be easily achieved (e.g., ...). Figure 2 ).

[0161] B) The identification of drug resistance genotypes employs a strategy combining two rules: First, for target genes with high similarity to other genes in the database (e.g., consistency exceeding 95%), the gene with the highest ARG-like read count (i.e., the total number of reads that can be compared to the target gene, or the number of specific reads plus the number of multiple alignment reads) is considered a true positive. For genes not ranked first, those with 100% coverage calculated based on precise alignment reads are retained as true positives. Second, for target genes with low similarity to other genes in the database (e.g., consistency below 95%), the primary criterion for a true positive result is the number of detected specific reads. Gene families without detected specific reads are considered false positives and directly filtered out.

[0162] C) Evaluate the accuracy of read-based detection for identifying drug resistance gene genotypes. Using the assembly-based detection results of drug resistance genes during BGWAS model training as a reference, statistically analyze the accuracy, sensitivity, and specificity of read-based detection of drug resistance genes or gene families selected by the BGWAS model, such as... Figure 3 The results show that the read-based detection process performs well in detecting most drug resistance genes or gene families.

[0163] 1.4 Define and calculate the positive / negative criterion Score, and perform ROC analysis in conjunction with simulation tests under different sequencing data volumes to determine the reporting rules and cutoff threshold.

[0164] 1.4.1 Based on the important genes (and gene families) and their weight coefficients obtained through screening, and according to the detection results of drug resistance genes in the samples, the Score index value is defined and calculated. The calculation formula is as follows:

[0165]

[0166] In the formula, arg_W i Genefamily_W represents the weighting coefficient for the corresponding genotype. i This represents the weight coefficient of the corresponding gene family. When a genotype (genetype, i.e., gt) is detected and the genotype weight coefficient is >0, the genotype weight coefficient is used for calculation; when a genotype is detected but the genotype weight coefficient is 0 or there is no weight coefficient, the gene family weight coefficient is used for calculation (if the gene family has no weight coefficient, it is recorded as 0).

[0167] Specifically, according to Figure 4The rules, based on read-based drug resistance gene genotyping results and gene weight coefficient matrix, and considering the positive consistency rate (PPV) between each gene and drug susceptibility results in the BGWAS model training set, are used to calculate the score. Figure 4 In Chinese: gt represents ARG type, i.e., drug resistance gene typing, and gf represents ARG family, i.e., drug resistance gene family.

[0168] 1.4.2 For simulation tests based on training set strains with different amounts of data, ROC analysis was performed based on the drug resistance gene detection results and the actual drug susceptibility results of the training set strains to evaluate the performance (AUC value) of the drug susceptibility prediction model under different sequencing data amounts and to determine the cutoff threshold.

[0169] Specifically, considering the influence of model performance and sequencing data volume, thresholds are set separately for reporting "drug resistance" and "susceptibility". The threshold for reporting "drug resistance" is set as the Score value corresponding to the maximum Youden index when the target strain genome sequencing data volume is sufficient and the model performance is stable (e.g., 30X), denoted as R_cutoff. The threshold for reporting "susceptibility" uses two sets of threshold standards. One is the R_cutoff determined under the condition of sufficient sequencing data volume (30X data volume, stable model), denoted as S_cutoff2, where the NPV must be above 0.9; otherwise, "susceptibility" cannot be reported. The second is the minimum sequencing data volume (denoted as gf_LOD1) corresponding to stable model performance, found directly based on the score value calculated using gene family weight coefficients. Under this data volume, the maximum Score value that satisfies an NPV exceeding 0.9 (this Score value must be less than or equal to S_cutoff2) is used as the threshold for reporting "susceptibility," denoted as S_cutoff1. Secondly, for various simulated sequencing data volumes between gf_LOD1 and 30X, we examined the feasibility of using S_cutoff2 as the reporting "sensitivity" threshold standard, i.e., finding the minimum data volume corresponding to NPV>0.9, denoted as gf_LOD2. Therefore, we can ultimately determine two sets of reporting "sensitivity" threshold standards: when the sequencing data volume is above gf_LOD2, S_cutoff2 is used as the reporting "sensitivity" threshold; when the sequencing data volume is between gf_LOD1 and gf_LOD2, S_cutoff1 is used as the reporting "sensitivity" threshold.

[0170] The specific reporting rules for drug susceptibility prediction for a specific antibiotic are as follows: When the score value is greater than R_cutoff, it is reported as "possibly resistant"; when the genome coverage of the detected pathogen is greater than or equal to gf_LOD2 and the score value is less than S_cutoff2, it is reported as "possibly sensitive"; or when the genome coverage of the detected pathogen is between gf_LOD1 and gf_LOD2 and the score value is less than S_cutoff1, it is reported as "possibly sensitive"; when the genome coverage of the detected pathogen is less than gf_LOD1 and the score is less than S_cutoff1, it is reported as " / ", i.e., unknown.

[0171] (ii) Species Attribution Prediction of Drug Resistance Genes

[0172] 2.1 Detection of target pathogen species and prediction of species attribution of drug resistance genes

[0173] 2.1.1 Target pathogen genome sequence alignment and calculation of the number of detected sequences, genome coverage, and coverage depth.

[0174] Commonly known clinical pathogens (Klebsiella pneumoniae, Escherichia coli, Acinetobacter baumannii, Pseudomonas aeruginosa, Enterobacter cloacae, etc.) were selected as target pathogens. The reference genomes of the target pathogens were searched and downloaded from the NCBI Genome Database and used as a reference sequence library for the identification of target pathogen species.

[0175] Using minimap2 software (v2.17), Illumina reads were aligned with the target pathogen reference genome sequence library set above (alignment parameter: -x sr-a--secondary=no-L). The number of detected sequences, genome coverage, and coverage depth of the aligned target pathogen species were then calculated. For alignment of nanopore sequencing reads, the parameter was set to -x map-ont-a--secondary=no-L.

[0176] Next, the total number of reads, genome coverage, and coverage depth of the detected pathogenic bacteria were statistically analyzed.

[0177] 2.1.2 Based on the detection results of drug resistance genes in the BGWAS model training set specimens, the copy number of drug resistance genes carried by the target pathogen species was statistically analyzed.

[0178] Based on the assembly-based drug resistance gene detection results of the training set samples during BGWAS model training, the detection distribution and copy change range of drug resistance genes and drug resistance gene families of the target pathogen species were statistically obtained.

[0179] 2.1.3 Based on the assumed gene-species attribution relationship, calculate the copy number of drug resistance genes and determine species attribution.

[0180] When assuming a gene-species correspondence for drug resistance, the main criteria are: a) whether the species annotation of the reference gene in the database includes the target species; if so, the hypothesis of the gene-species attribution is accepted. b) if a is not met, the mediation mode annotation of the reference gene in the database is checked to see if it includes plasmid-mediated modes; if so, the hypothesis of the gene-species attribution is accepted. c) if a and b are not met, the species origin is inferred based on the species annotation of ARG-like reads. Then, the copy number of the drug resistance gene is calculated using the following formula:

[0181]

[0182] If the calculated copy number of the drug resistance gene falls within the normal copy number range of the target gene family obtained from the BGWAS model training set, then the assumed gene-species attribution relationship is accepted; otherwise, it is rejected.

[0183] (iii) Metagenomic sequencing of clinical specimens for drug resistance gene detection and validation

[0184] Clinical specimens that underwent culture and drug susceptibility testing were collected and transported to the medical laboratory for pretreatment, nucleic acid extraction, library construction, and sequencing (Illumina CN500 SE75 or Nanopore sequencing). The reads obtained from sequencing were then compared and analyzed with the target species genome and drug resistance gene database to identify the pathogenic bacteria in the samples and the drug resistance genes they carry.

[0185] 3.1 Perform quality control filtering and remove human nucleic acid sequences from the raw sequencing low-quality sequences.

[0186] The sequence data obtained from the Illumina platform are processed as follows:

[0187] a. Use bcl2fastq (v2.20.0.422) to process the sequencing data, converting the BCL format data into fastq format sequence data;

[0188] b. Use the fastp (v0.19.5) software to filter the obtained raw fastq sequence data (parameter setting: -q15 -u 40 -l read_length*0.67) to remove low-quality and short sequences; at the same time, use the komplexity (v0.3.6) software to calculate the sequence information complexity (parameter setting: -Ft 0.4) and filter out sequences with low complexity.

[0189] c. Align the clean sequences obtained by quality control filtering with the human reference genome sequence (human_38) using the bowtie2 (v2.3.4.3) software (parameter settings: --mm --very-sensitive -k 1) to filter out human-derived sequences.

[0190] 3.2 Detect and identify the drug-resistant genes contained in the sample and their species attribution

[0191] According to step 1.3.2, use the blastn software to perform drug-resistant gene alignment and annotation statistics on the sample sequences, and according to step 2.1.1, use the minimap2 software to perform genome alignment of the target species and calculate the detected genome coverage. Then, according to step 2.1.3, evaluate and predict the species attribution of the detected drug-resistant genes.

[0192] 3.3 For the pathogenic bacteria of concern and detected, according to the detected situation of their drug-resistant genes, calculate the score value of the target species - antibiotic drug according to the score calculation method defined in 1.4.1, and compare it with the cutoff value: when score >= cutoff, it is predicted as R; when score < cutoff, if the detected pathogenic genome coverage is higher than the minimum genome coverage or data volume required for model stability determined based on step 1.4.2 (such as > 40%), it is predicted as S, otherwise report as " / " (indicating unknown)

[0193] 3.4 Finally, compare the above drug sensitivity prediction results with the actual drug sensitivity test results of the clinical specimens collected at the same time, count the accuracy rate of the prediction results and the proportion of the number of samples with effective reports, and evaluate the performance of the drug resistance detection process.

[0194] Example 2: Perform drug-resistant gene detection and phenotypic prediction analysis on clinical specimens for Klebsiella pneumoniae

[0195] 1. Based on the Klebsiella pneumoniae strain genome, BGWAS screened out important drug-resistant genes related to antibiotics and their corresponding gene families, and calculated the weight coefficients of the genes or gene families

[0196] Collect and download the genome data of Klebsiella pneumoniae strains with drug sensitivity test result information from the NCBI NDARO database and the PATRIC database, and finally obtain 3072 strain samples (2410 and 662 cases for the training set and validation set) after screening and filtering. Then, based on the machine learning method, screen out the important genes related to antibiotic drug resistance and their weight coefficient matrix, and calculate the weight coefficients of the families to which these important genes belong according to the formula in technical solution 1.2. The results are shown in the following table:

[0197]

[0198]

[0199]

[0200]

[0201] 2. Based on the training set of 2410 Klebsiella pneumoniae strain genomes, the read-based drug resistance gene detection process was tested and validated using 75bp short reads (NGS sequencing platform) following steps 1.3 of the technical scheme. Different data gradients were simulated: 0.05X, 0.1X, 0.2X, 0.3X, 0.4X, 0.5X, 0.6X, 0.7X, 0.8X, 0.9X, 1X, 2X, 3X, 5X, 10X, and 30X. Drug resistance gene detection was then performed to obtain the detection results for each simulated sample. Score values ​​were calculated according to steps 1.4.1 of the technical scheme, and ROC curve analysis was conducted to obtain the AUC values ​​of each antibiotic under different data volumes. The AUC value variation curves of the model performance are then plotted as follows. Figure 5 When the data volume is 30X, the model performance has stabilized. The AUC values ​​for each antibiotic model at this point are as follows: Figure 6 Following the steps in technical solution 1.4.2, the reporting rules and cutoff thresholds for each antibiotic were finally determined as shown in the table below.

[0202] mNGS antimicrobial susceptibility prediction cutoff thresholds for six antibiotics in Klebsiella pneumoniae.

[0203]

[0204] Note: " / " indicates that "sensitive" cannot be reported.

[0205] 3. A total of 48 clinical specimens containing Klebsiella pneumoniae were collected and cultured, and stored at -80℃. Nucleic acid was then extracted from these 48 specimens to construct metagenomic second-generation (intercalation fragment length 200-400bp) libraries, which were then sequenced using Illumina NextSeq CN500 SE75. The resulting data were used to identify the pathogen and its drug resistance genes, and score values ​​and drug susceptibility predictions were calculated. Finally, the detection and identification results of the target pathogen and its drug resistance genes, as well as the drug susceptibility prediction results, were obtained for each specimen. Some sample results are shown in the table below:

[0206]

[0207]

[0208]

[0209]

[0210] Note: ND indicates not detected, and " / " indicates unknown.

[0211] The statistics show that the accuracy rate of drug sensitivity prediction and the proportion of reportable samples are as follows:

[0212]

[0213] The results show that the present invention can effectively and accurately identify pathogenic bacteria and their drug resistance genes in clinical specimens, as well as predict the drug susceptibility of antibiotics, and can be used to assist in the clinical detection and diagnosis of drug-resistant bacteria.

[0214] Further, the above results show that AAC(3)-IIe, AAC(3)-IV, AAC(3)-IId, rmtC, armA, rmtF, rmtB, AAC(6')-33, and ANT(2”)-Ia are important resistance genes for gentamicin. Considering factors such as gene weight, gene family weight, gene occurrence frequency, and all possible mechanisms of resistance, in practice, the high-frequency and high-weight genes that mainly mediate resistance, AAC(3)-IIe, AAC(3)-IV, AAC(3)-IId, rmtC, armA, rmtF, rmtB, AAC(6')-33, and ANT(2”)-Ia, can be detected simultaneously. If all the test results are negative, it can be inferred that the drug is sensitive, that is, the purpose of drug sensitivity detection, especially sensitive detection, is achieved.

[0215] Further, the above results indicate that AAC(3)-IV, AAC(3)-IId, AAC(6')-Ib', AAC(6')-Ib-cr, AAC(6')-Ib-Hangzhou, AAC(6')-Ib4, mphE, ANT(2”)-Ia, and aadA24 are important resistance genes for tobramycin. Considering factors such as gene weight, gene family weight, gene occurrence frequency, and all possible mechanisms of resistance development, in practice... The AAAC(3)-IV, AAC(3)-IId, AAC(6')-Ib', AAC(6')-Ib-cr, AAC(6')-Ib-Hangzhou, AAC(6')-Ib4, mphE, ANT(2”)-Ia, and aadA24 genes, which are the main mediators of drug resistance, can be detected simultaneously. If all the test results are negative, it can be inferred that the drug is sensitive, thus achieving the purpose of drug sensitivity detection, especially sensitivity detection.

[0216] Further analysis of the above results reveals that, against ceftazidime, CTX-M-55, CTX-M-11, CTX-M-15, SHV-155, SHV-5, SHV-11, SHV-12, SHV-76, SHV-30, SHV-53, SHV-124, SHV-182, DHA-1, KPC-3, and KPC-2 are important resistance genes. This is further supported by considering gene weights, gene family weights, gene frequency, and all possible factors contributing to resistance. Mechanisms and other factors can be considered in practice. This can be addressed by detecting high-frequency, high-weight genes that primarily mediate drug resistance, such as CTX-M-55, CTX-M-11, CTX-M-15, SHV-155, SHV-5, SHV-11, SHV-12, SHV-76, SHV-30, SHV-53, SHV-124, SHV-182, DHA-1, KPC-3, and / or KPC-2. If all test results are positive, drug resistance can be inferred. However, according to the 30X genome simulated reads test results, a suitable threshold cannot be found where the calculated score is less than the threshold, resulting in a high NPV (e.g., >0.9), thus failing to achieve the goal of sensitive detection.

[0217] Further analysis of the above results reveals that, for compound sulfamethoxazole, dfrA12, dfrA15, dfrA17, dfrA19, dfrA30, dfrA8, dfrA5, dfrA15b, dfrA14, dfr22, dfrA27, and dfrA1 are important resistance genes. Considering factors such as gene weight, gene family weight, gene frequency, and all possible mechanisms of resistance, in practice, simultaneous detection of the high-frequency and high-weight dfrA12, dfrA15, dfrA17, dfrA19, dfrA30, dfrA8, dfrA5, dfrA15b, dfrA14, dfr22, dfrA27, and / or dfrA1 genes that primarily mediate resistance can be performed. If all test results are positive, resistance can be inferred. This achieves drug susceptibility testing, especially sensitive detection. Meanwhile, according to the 30X genome simulated reads test results, there is no suitable threshold that results in a high NPV (e.g., >0.9) when the calculated Score is less than the threshold, thus failing to achieve the purpose of sensitive detection.

[0218] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of this application.

Claims

1. A method for constructing a model to predict the species attribution of drug-resistant genes, characterized in that, The method includes the following steps: Step 1): Alignment of the target pathogen's genome sequence and calculation of the number of detected sequences, genome coverage, and coverage depth; Step 2): Based on the detection results of drug resistance genes in the BGWAS model training set specimens, count the copy number of drug resistance genes carried by the target pathogen species; Step 3): Based on the assumed gene-species attribution relationship, calculate the copy number of the drug resistance gene and determine the species attribution: When assuming a gene-species correspondence for drug resistance, the following criteria are used: a) Does the species annotation of the reference gene in the database include the target species? If so, the hypothesis of the gene-species attribution is accepted. b) If a is not met, the mediation mode annotation of the reference gene in the database is checked to see if it includes plasmid-mediated modes. If so, the hypothesis of the gene-species attribution is accepted. c) If a and b are not met, the species origin is inferred based on the species annotation of ARG-like reads, and the copy number of the drug resistance gene is calculated using the following formula: If the calculated copy number of the drug resistance gene falls within the normal copy number range of the target gene family obtained from the BGWAS model training set, then the assumed gene-species attribution relationship is accepted; otherwise, it is rejected.

2. The model construction method according to claim 1, characterized in that, Step 1) specifically refers to: Commonly known clinical pathogens were selected as target pathogens. The reference genome of the target pathogen was searched and downloaded from the NCBI genome database and used as a reference sequence library for the identification of the target pathogen species. Each sequenced read was compared with the above reference sequence library, and the number of detected sequences, genome coverage, and coverage depth of the target pathogen species were calculated. The total number of detected pathogen sequences, genome coverage, and coverage depth were then statistically obtained.

3. The model construction method according to claim 1, characterized in that, Step 2) specifically refers to: Based on the assembly-based drug resistance gene detection results of the training set samples during BGWAS model training, the detection distribution and copy change range of drug resistance genes and drug resistance gene families of the target pathogen species were statistically obtained.

4. A model construction device for predicting the species attribution of drug-resistant genes, characterized in that, include: At least one memory for storing programs; At least one processor is configured to load the program to execute the model building method as described in any one of claims 1-3.

Images