Method for identifying large insertions in target genomic regions and uses thereof

Abstract

The application discloses a method for identifying large fragment sequence insertion of a target genomic region and application thereof, and belongs to the technical field of bioinformatics. In view of the problem that large fragment insertion variation is difficult to be accurately recognized in clinical metagenomic sequencing due to short sequencing read length, insufficient coverage and other factors, the application proposes to construct a reference sequence which can represent the insertion variation by means of manual construction, and to realize efficient identification of the insertion event of the target genomic region by combining a short read-based fast alignment process. The method overcomes the dependence of existing structural variation detection tools on high sequencing depth and long read length, has the advantages of fast identification speed, high sensitivity and high accuracy, and is suitable for rapid screening of large fragment insertion related to drug resistance mechanism in clinical samples. Meanwhile, the method can be popularized for insertion variation analysis of other pathogen drug resistance related genes or genomic regions, and has a good clinical application prospect.

Description

Technical Field

[0001] This invention belongs to the field of biomedical technology, and in particular relates to a method for identifying large sequence insertions in a target genomic region and its application. Background Technology

[0002] Clinical metagenomic sequencing (mNGS), as an emerging method for detecting infectious pathogens, boasts advantages such as high sensitivity, broad pathogen coverage, and rapid response. It also offers potential value for further analysis of drug resistance genes and indicative of bacterial resistance. It is well known that bacterial antibiotic resistance is typically encoded genetically, and resistance can develop through various mechanisms, including overexpression or duplication of existing genes, point mutations, or acquisition of novel genes via horizontal gene transfer (HGT). Accordingly, bacterial resistance to antibiotics can be predicted by analyzing genomic characteristics such as the presence or variation of resistance genes. Researchers have developed models based on clinical mNGS data to directly predict antibiotic resistance in Acinetobacter baumannii, Klebsiella pneumoniae, and Pseudomonas aeruginosa, including carbapenems (imipenem / meropenem), cephalosporins (cefotaxime / ceftriaxone / ceftazidime), and fluoroquinolones (levofloxacin / ciprofloxacin), with a concordance rate exceeding 90% with conventional culture susceptibility testing results. However, resistance to some other antibiotics is mainly caused by large-segment insertion variations in the genome. For example, an IS sequence insertion in the mgrB gene of Klebsiella pneumoniae mediates the upregulation of PhoP / PhoQ expression, leading to colistin resistance; an IS sequence insertion in the acrR gene of Escherichia coli mediates the upregulation of the efflux pump AcrAB, leading to quinolone and β-lactam resistance. Due to the short length of mNGS sequencing reads (e.g., 50 bp or 75 bp) and the low genome sequencing depth, existing conventional bioinformatics alignment software such as DELLY, Manta, BWA, and KMA often require a large amount of sequencing data, making it difficult to identify large-segment (e.g., >50 bp) genomic insertion variations. Summary of the Invention

[0003] In view of the above limitations, in order to achieve rapid identification of large-fragment sequence insertion variants in the genome of pathogens based on clinical mNGS data, this invention artificially constructs a reference sequence that can represent large-fragment insertion variants and a corresponding read-based rapid alignment identification process. The effectiveness of the method is evaluated using the colistin-associated drug resistance gene mgrB of Klebsiella pneumoniae as an example.

[0004] This invention is implemented as follows: a method for identifying large sequence insertions in a target genomic region, using short sequencing read alignment, includes the following steps:

[0005] S1. Collect the genome and phenotypic information of the target pathogen strain, and perform quality control to remove genomes with poor assembly quality;

[0006] S2. Select the wild-type target pathogen genome region as the reference sequence, compare the genomes of all strains collected in S1 with the reference sequence, find the insertion site and the IS sequence at the insertion site; cut a short sequence of 10-100 bp from the end of the IS sequence and manually insert it into the corresponding insertion site in the reference sequence to construct a representative sequence that can indicate the insertion of large fragment sequences in the target genome region.

[0007] S3. Establish a genotyping and identification process for representative sequences of target genomic regions based on the alignment of genome contigs and WGS short sequencing reads, and evaluate and determine the optimal length of the short sequences inserted to construct representative sequences of target genomic regions and the accuracy of the identification process.

[0008] Furthermore, in S1, genome and phenotypic information of the target pathogen strain are collected from the NCBI NDARO, BV-BRC public database or local hospitals, and genomes with poor assembly quality, such as those with more than 500 genome contigs, genome length exceeding 1.5 times the average length or less than 0.5 times the average length, are removed.

[0009] Furthermore, in S2, the alignment results with the reference sequence are divided into the following five modes:

[0010] M0: Only one alignment exists, and it covers the entire region of the reference sequence. In this case, no IS sequence is inserted.

[0011] M1-1: There is only one alignment, which covers the right region of the reference sequence. The alignment breakpoint is the insertion site of the IS sequence, and the unaligned query contig sequence to the left of the alignment breakpoint is the IS sequence.

[0012] M1-2: There is only one alignment, which covers the left region of the reference sequence. The alignment breakpoint is the insertion site of the IS sequence, and the unaligned query contig sequence to the right of the alignment breakpoint is the IS sequence.

[0013] M2-1: There are two alignments for the same contig, covering the left L1 and right R2 regions of the reference sequence, respectively. If there is an overlap between the two alignments, there is one insertion site. Further, for the overlapping region, if the number of SNPs and InDels in the L1 alignment is higher than that in the R2 alignment, the insertion site is located at the alignment breakpoint on the left side of the overlapping region; otherwise, the insertion site is located at the alignment breakpoint on the right side of the overlapping region. If there is no overlap between the two alignments, there are two insertion sites. Further, the insertion sites are located at the alignment breakpoint on the right side of the L1 alignment and the alignment breakpoint on the left side of the R2 alignment. The unaligned sequence between L1 and R2 in the query is the IS sequence.

[0014] M2-2: Two different contigs each have one alignment, resulting in two alignments, covering the left L1 and right R2 regions of the reference sequence, respectively. If there is an overlap between the two alignments, there is one insertion site. Further, regarding the overlap region, if the number of SNPs and InDels in the L1 alignment is higher than that in the R2 alignment, the insertion site is located at the alignment breakpoint on the left side of the overlap region; otherwise, the insertion site is located at the alignment breakpoint on the right side of the overlap region. If there is no overlap between the two alignments, there are two insertion sites. Further, the insertion sites are located at the alignment breakpoint on the right side of the L1 alignment and the alignment breakpoint on the left side of the R2 alignment. The unaligned region on the right side of the query contig corresponding to the L1 alignment and the unaligned region on the left side of the query contig corresponding to the R2 alignment constitute the IS sequence.

[0015] Furthermore, in S3, the assembly-based and read-based alignment identification processes were used to complete the identification and analysis of the target genomic region sequence types of all strain genomes obtained in S1. Using the alignment breakpoint analysis results obtained in S2 as a reference standard, the identification accuracy of the assembly-based and read-based processes and the optimal short sequence insertion length when constructing representative sequences were compared and confirmed.

[0016] Furthermore, the assembly-based alignment and identification process involves using Blastn software to directly align the contig sequences of all S1 strain genomes with the wild-type reference sequence of the target genome region and artificially constructed representative sequences. The Hit with the highest score is selected as the final alignment result, which is then identified as the final target gene or region sequence type. Simultaneously, the m0 alignment format file is parsed to obtain SNPs and InDels information on the target gene or region.

[0017] Furthermore, the read-based alignment and identification process is as follows: First, based on all genomes collected by S1, the ART software is used to simulate the 10X WGS reads sequence of each strain; then, the KMA software is used to align the simulated 10X WGS reads sequence of the strain with the wild-type reference sequence of the target genome region and the artificially constructed representative sequence, and the final results are identified and annotated.

[0018] Furthermore, the identification and annotation process is as follows: First, for each query read, only the best alignment result is retained. Then, the number of read sequences aligned with each reference sequence is counted and sorted from high to low according to the number of aligned reads. If the number of reads in the top 2 is less than the number of reads in the top 1, then the sequence type corresponding to the top 1 reference sequence is finally identified. If the number of reads in the top 2 is equal to the number of reads in the top 1, then the identification resolution is recorded as family level, that is, at this time, the specific target genomic region sequence type cannot be identified.

[0019] A computer-readable storage medium includes a stored computer program, wherein, when the computer program is executed, it controls the device on which the computer-readable storage medium is located to perform the method described above.

[0020] A computer device includes a memory, a processor, and a program stored in the memory and executable thereon, the program being executed by the processor to implement the steps of the method described above.

[0021] The advantages and technical effects of this invention are as follows:

[0022] 1. The reference sequence library and corresponding short read-based alignment identification process constructed in this invention can indicate large fragment sequence insertions, which solves the problem that clinical specimen mNGS cannot directly and quickly identify IS sequence insertion variations in target genomic regions due to short read lengths and shallow sequencing. It can also effectively identify colistin resistance caused by large fragment sequence insertions in the mgrB gene region of Klebsiella pneumoniae.

[0023] 2. The method of the present invention is also applicable to the detection of IS sequence insertions in other microbial species or phenotype-related genes. Attached Figure Description

[0024] Figure 1 This is a schematic diagram of the alignment pattern of the target genome region and the construction of representative sequences that can indicate IS sequence insertion;

[0025] Figure 2 This is a frequency distribution map of the mgrB gene of Klebsiella pneumoniae and the IS sequence insertion site in its 150bp upstream region;

[0026] Figure 3 The curves show the impact of the extracted IS sequence length on the use of assembly-based and read-based alignment analysis procedures for sequence typing. Detailed Implementation

[0027] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0028] like Figure 1 As shown, the method for identifying large sequence insertions in a target genomic region according to the present invention includes the following steps:

[0029] Step 1: Collect genomic and phenotypic information data of the target pathogen strain. For the pathogen and phenotype of interest (e.g., drug resistance phenotype), collect genomic and phenotypic data of the strain from public databases (NCBI NDARO, BV-BRC, etc.) or local partner hospitals. Then, according to the strain's Biosample number, remove redundancy from the collected genomes and perform quality control and screening. Strain quality control filtering criteria: remove genomes with poor assembly quality (contig number > 500, genome length exceeding the average length by 1.5 times or falling below the average length by 0.5 times, predicted ORF number exceeding the average number by 1.5 times or falling below the average number by 0.5 times); remove genomes that, after alignment, do not carry the target species' 16S rRNA gene; and remove genomes with multiple inconsistent phenotypic data points.

[0030] Step 2: Select the wild-type target pathogen genome region as the reference sequence. Using Blasn software, align and analyze the genomes of all strains collected in Step 1 with the reference sequence to find the insertion site and the IS sequence at the insertion site. Then, truncate a short sequence of 10-100 bp from the end of the IS sequence and manually insert it into the corresponding insertion site in the reference sequence to construct a representative sequence that can indicate the insertion of large fragment sequences in the target genome region. The alignment results with the reference sequence are divided into the following five modes:

[0031] M0: Only one alignment exists, and it covers the entire region of the reference sequence. In this case, no IS sequence is inserted.

[0032] M1-1: There is only one alignment, which covers the right region of the reference sequence. The alignment breakpoint is the insertion site of the IS sequence, and the unaligned query contig sequence to the left of the alignment breakpoint is the IS sequence.

[0033] M1-2: There is only one alignment, which covers the left region of the reference sequence. The alignment breakpoint is the insertion site of the IS sequence, and the unaligned query contig sequence to the right of the alignment breakpoint is the IS sequence.

[0034] M2-1: There are two alignments for the same contig, covering the left L1 and right R2 regions of the reference sequence, respectively. If there is an overlap between the two alignments, there is one insertion site. Further, for the overlapping region, if the number of SNPs and InDels in the L1 alignment is higher than that in the R2 alignment, the insertion site is located at the alignment breakpoint on the left side of the overlapping region; otherwise, the insertion site is located at the alignment breakpoint on the right side of the overlapping region. If there is no overlap between the two alignments, there are two insertion sites. Further, the insertion sites are located at the alignment breakpoint on the right side of the L1 alignment and the alignment breakpoint on the left side of the R2 alignment. The unaligned sequence between L1 and R2 in the query is the IS sequence.

[0035] M2-2: Two different contigs each have one alignment, resulting in two alignments, covering the left and right regions of the reference sequence respectively. If there is an overlapping region between the two alignments, there is one insertion site. Further, regarding the overlapping region, if the number of SNPs and InDels in the L1 alignment is higher than that in the R2 alignment, the insertion site is located at the alignment breakpoint on the left side of the overlapping region; otherwise, the insertion site is located at the alignment breakpoint on the right side of the overlapping region. If there is no overlapping region between the two alignments, there are two insertion sites. The insertion sites are located at the alignment breakpoint on the right side of the L1 alignment and the alignment breakpoint on the left side of the R2 alignment. The unaligned region on the right side of the query contig corresponding to the L1 alignment and the unaligned region on the left side of the query contig corresponding to the R2 alignment constitute the IS sequence.

[0036] Step 3: Establish a genotyping and identification process for representative sequences of the target genome region based on genome contigs and WGS short sequencing read alignment, and evaluate the optimal length of the short sequences inserted to construct the representative sequences of the target genome region, as well as the accuracy of the identification process. Details are as follows:

[0037] Assembly-based alignment and identification process: Using Blastn software, all strain genome contig sequences collected in step 1 are directly aligned with the wild-type reference sequence of the target genome region and artificially constructed representative sequences. The HIT with the highest score is selected as the final alignment result, which is identified as the final target gene or region sequence type. At the same time, the m0 alignment format file is parsed to obtain SNPs and InDels information on the target gene or region.

[0038] Read-based alignment and identification process: First, based on all genomes collected in step 1, ART software is used to simulate 10X WGS read sequences for each strain. Then, KMA software is used to align the simulated 10X WGS read sequences with wild-type reference sequences of the target genome region and artificially constructed representative sequences, and the final results are identified and annotated. The detailed identification and annotation process is as follows: First, for each query read, only the best alignment result is retained. Then, the number of read sequences aligned with each reference sequence is counted and sorted from highest to lowest read count. If the number of reads in the second-highest ranking is less than the number of reads in the highest ranking, the read is finally identified as the sequence type corresponding to the highest ranking reference sequence. If the number of reads in the second-highest ranking is equal to the number of reads in the highest ranking, the identification resolution is recorded as family level, meaning that a specific target genome region sequence type cannot be identified at this time.

[0039] Using assembly-based and read-based alignment identification workflows, the target genomic region sequence types of all strains collected in step 1 were identified and analyzed. Then, using the alignment breakpoint analysis results obtained in step 2 as a reference standard, the identification accuracy of the assembly-based and read-based workflows and the optimal short sequence insertion length for constructing representative sequences were compared and confirmed.

[0040] Taking the colistin resistance phenotype of Klebsiella pneumoniae as an example, the genomic data of Klebsiella pneumoniae strains and their colistin susceptibility test results were collected from clinical settings to verify and evaluate the performance of this method. The representative sequences of the target genomic regions constructed in step 2 were integrated with the CARD public resistance library. Then, WGS analysis was performed on the collected clinical culture isolates to evaluate the performance of the read-based identification process.

[0041] Example 1: Identification of large fragment sequence insertion characteristics in the mgrB gene region of Klebsiella pneumoniae and prediction of colistin resistance.

[0042] Step 1: Search and download the genome and antibiotic susceptibility phenotype data of Klebsiella pneumoniae from public databases.

[0043] Download from the NCBI NDARO database: Open the website, enter "Klebsiella pneumoniae" in the search bar to retrieve information about Klebsiella pneumoniae, then in the Matched Isolates sub-window, click "Choose columns" and select "ASTpheotypes" to display the information in this column. Next, download the table data for the entire window, organize the Klebsiella pneumoniae strains with drug susceptibility phenotype data, and download the genome sequences in batches from the NCBI genome database based on the Assembly ID information.

[0044] Download from the BV-BRC database platform: On a networked Linux server, use the command `wget -c` to download all strain drug susceptibility data files from the BV-BRC database. Then, find all rows containing "Klebsiella pneumoniae" in the `genome_name` column, and use the command `wget -qNc` to download the corresponding Klebsiella pneumoniae strain genome data based on the corresponding `genome_id` information.

[0045] Next, the collected strain genomes and drug susceptibility data were subjected to quality control filtering. The specific quality control filtering criteria included: a) removing genomes with multiple but inconsistent drug susceptibility test (AST) results; b) genomes with an N50 value less than 5000 bp or more than 2000 contigs; c) genomes with a length significantly deviating from the average length of the target species (more than 1.5 times or less than 0.5 times); d) predicted gene counts significantly lower (less than 0.5 times the average predicted count of the target species) or failure to meet gene density requirements (the ratio of predicted CDS counts to genome length per kilobase pair is not between 0.5 and 1.5); e) average nucleotide similarity (ANI) with the target species reference genome less than 0.95, or classification annotations obtained by aligning contigs to the NCBI NT database that are inconsistent with the target species; f) failure to detect the 16S rDNA sequence of the target species.

[0046] Based on Biosample, the genomes of the strains downloaded from the NCBI NDARO and BV-BRC databases were de-merged for redundancy, resulting in a total of 2765 non-redundant Klebsiella pneumoniae genomes and their corresponding drug susceptibility test (AST) results. Specific strain drug resistance data are shown in Table 1.

[0047] Table 1

[0048] strain name Drug Classification Drug Name Total number of strains Antimicrobial susceptibility is the number of drug-resistant strains. Antimicrobial susceptibility refers to the number of susceptible bacterial strains. Klebsiella pneumoniae Colistin Colistin 2765 531 2234

[0049] Step 2: Identify the polymorphism of the mgrB gene and upstream promoter region sequence of Klebsiella pneumoniae, as well as the distribution of IS sequence insertion sites.

[0050] 2.1 The mgrB gene and its upstream 150bp region from the colistin-sensitive wild-type Klebsiella pneumoniae genome were selected as the reference sequence. All other genomes were aligned with this reference sequence to obtain the alignment results and variation information of the hit regions. All 2765 collected Klebsiella pneumoniae genome sequences were aligned with the mgrB gene and its upstream 150bp reference sequence using blastn (version 2.9.0+, parameters: -evalue1e-5 -num_threads 6 -outfmt 0 -num_alignments 10000 -task blastn). Hits with an identity less than 90% or a reference gene coverage less than 60% were filtered out. The best hit (first hit) was selected as the final alignment result for each contig region. The m0 format alignment information was parsed to obtain the qstart, qend, sstart, send, and variation information from the query contig and the reference sequence.

[0051] 2.2 Analysis of alignment patterns and IS sequence insertion site distribution in the mgrB gene region. Based on the number of alignment hits, alignment breakpoints, and alignment coverage on the reference sequence between the strain genome and the reference sequence, the alignment results of each strain genome were summarized into the following five patterns (see...). Figure 1 ), and determine the insertion site and IS sequence.

[0052] 1) Mode M0: Only one alignment exists, and it covers the entire region of the reference sequence. In this case, no IS sequence is inserted.

[0053] 2) Pattern M1-1: Only one alignment exists, and it mainly covers the right side of the reference sequence. The alignment breakpoint is the insertion site of the IS sequence, and the unaligned query contig sequence to the left of the alignment breakpoint is the IS sequence.

[0054] 3) Mode M1-2: Only one alignment exists, and it mainly covers the left part of the reference sequence. The alignment breakpoint is the insertion site of the IS sequence, and the unaligned query contig sequence to the right of the alignment breakpoint is the IS sequence.

[0055] 4) Pattern M2-1: Two alignments exist within the same contig, covering the left (L1) and right (R2) regions of the reference sequence, respectively. If there is an overlap between the two alignments, there is one insertion site. Further, regarding the overlap region, if the number of SNPs and InDels in the L1 alignment is higher than that in the R2 alignment, the insertion site is located at the alignment breakpoint to the left of the overlap region; otherwise, the insertion site is located at the alignment breakpoint to the right of the overlap region. If there is no overlap between the two alignments, there are two insertion sites. The insertion sites are located at the alignment breakpoint to the right of the L1 alignment and the alignment breakpoint to the left of the R2 alignment. The unaligned sequence between L1 and R2 on the query is the IS sequence.

[0056] 5) Pattern M2-2: Two different contigs each have one alignment, resulting in two alignments that cover portions of the left and right sides of the reference sequence, respectively. If there is an overlap between the two alignments, there is one insertion site. Further, regarding the overlap region, if the number of SNPs and InDels in the L1 alignment is higher than that in the R2 alignment, the insertion site is located at the alignment breakpoint to the left of the overlap region; otherwise, the insertion site is located at the alignment breakpoint to the right of the overlap region. If there is no overlap between the two alignments, there are two insertion sites. Specifically, the insertion sites are located at the alignment breakpoint to the right of the L1 alignment and the alignment breakpoint to the left of the R2 alignment. The unaligned region on the right side of the query contig corresponding to the L1 alignment and the unaligned region on the left side of the query contig corresponding to the R2 alignment constitute the IS sequence.

[0057] Regarding the presence of IS sequence insertions in the mgrB gene region, the distribution of insertion sites in this region across all strain genomes was statistically analyzed. The results showed that high-frequency IS sequence insertion sites included 220bp (mgrB 70bp), 106bp (mgrB -44bp), 282bp (mgrB 132bp), and 226bp (mgrB 76bp), etc. (e.g.) Figure 2 ).

[0058] Step 3: Artificially construct representative sequences of the mgrB gene region indicating the presence of IS insertions. Using the reference sequence of the mgrB gene region as a template, for the insertion site, extract short sequences of a certain length (e.g., 10bp, 12bp, 15bp, 18bp, 20bp, 21bp, 24bp, 27bp, 30bp, 50bp, 100bp, etc.) from the unaligned IS sequences on the query contig near the insertion site, and manually insert them into the corresponding positions in the template. This yields representative sequences that indicate the presence of large-fragment IS insertions in the mgrB gene region. Specifically, for M2-1 or M2-2, since there are two alignments, two representative sequences can be obtained accordingly.

[0059] Step 4: For the constructed representative sequence of the mgrB gene region, two alignment and identification processes were established: Assembly-based and Read-based. The identification accuracy of the Read-based process and the impact of the length of the manually inserted fragment were evaluated.

[0060] 4.1 Establishing an Assembly-Based Alignment and Identification Process. Using the genomic sequence (contig sequence) of a Klebsiella pneumoniae strain as input, Blasten (version 2.9.0+, parameters: -evalue 1e-5 -num_threads 6 -outfmt 0 -num_alignments 10000 -task blastn) was used to directly align the contig sequence with a reference sequence library. The best alignment result (highest bit-score) for each contig was used as the final mgrB gene region sequence typing result for that strain. Simultaneously, by parsing the m0 format file generated by Blasten (-outfmt 0), SNPs (single nucleotide polymorphisms) and InDels (insertion / deletion variants) information within the mgrB gene region were extracted.

[0061] 4.2 Establishing a Read-based Alignment and Identification Process. 1) Simulated Sequencing Data: Based on the collected genome sequences of all Klebsiella pneumoniae strains, ART software (v2.5.8, parameters: -ss NS50 -l 50 -f 10 -na) was used to simulate 10X whole-genome sequencing (WGS) reads (Single-End, 50bp). 2) Read Alignment and Genotyping: KMA software (v1.4.9) was used to align the simulated reads with the reference sequence library. 3) Result Annotation: First, only the best alignment result was retained for each query read. Then, the number of reads matching each reference sequence was counted and sorted in descending order. If the number of reads for the Top1 reference sequence is higher than that for the Top2 (i.e., Top2 < Top1), it is identified as the sequence genotype corresponding to Top1. If the number of reads for Top1 and Top2 are equal, the specific sequence genotype information cannot be accurately identified at this time and is marked as "family-level genotyping," i.e., it cannot be accurate to the subtype.

[0062] 4.3 Based on representative sequences of the mgrB gene region constructed using insertion fragments of different lengths, the genomes of all collected strains were analyzed using assembly-based and read-based alignment and identification procedures. The consistency rates between the assembly-based procedure and the initial insertion site analysis results, as well as between the assembly-based and read-based identification results, were statistically analyzed. The results show (e.g.) Figure 3 The assembly-based approach consistently maintained 100% agreement with the initial insertion site analysis results. In contrast, the read-based approach showed the highest agreement rate (99.02%) when manually inserting an 18bp sequence to construct a representative sequence, but the agreement rate decreased when the manually inserted sequence length exceeded 30bp. Therefore, an 18bp IS sequence was chosen to construct the representative sequence for the mgrB region.

[0063] Step 5: Collect clinically cultured isolates and perform WGS sequencing verification. Fifty-four clinically cultured Klebsiella pneumoniae isolates were collected from domestic hospitals and sent to WGS for sequencing. Of these, 21 were colistin-resistant and 31 were susceptible. Simultaneously, the constructed mgrB region representative sequence, the mgrB region reference sequence, and the integrated resistance public library genes (including the mcr1 gene, and Klebsiella pneumoniae Phop / PhoQ / pmrB / pmrA / crrB genes) were merged as the final reference gene sequence. Then, a read-based analysis was used to identify the mgrB gene region sequence genotyping of these 54 strains, predicting colistin resistance caused by IS insertion, and comparing the results with the culture drug susceptibility results. The results showed that two strains had large fragment insertions, consistent with the culture drug susceptibility results. (See Table 2).

[0064] Table 2

[0065]

[0066] This invention analyzes large-fragment insertion sequences in the genome, extracts short sequences of a certain length from the ends of the inserted sequences, and artificially inserts them into the original wild-type genome reference sequence template. This constructs representative sequences that can indicate the presence of large-fragment sequence insertions in the genome. A corresponding read-based alignment analysis and identification process is established, enabling rapid identification of large-fragment sequence insertion variations in the genome of target pathogens based on strain WGS or clinical specimen mNGS sequencing reads. Taking the Klebsiella pneumoniae colistin resistance gene mgrB as an example, a large amount of colistin-resistant strain genome data is first collected from public libraries. Through direct alignment analysis of genome long contigs, the alignment pattern with the mgrB gene and its upstream promoter region (i.e., the -150bp region) is determined, thereby inferring the large-fragment IS sequence insertion site and its IS terminal sequence in the mgrB gene region. Then, a short sequence of a certain length was extracted from the unaligned IS sequence on the query contig at the insertion site and manually inserted into the mgrB gene region to construct a representative resistance sequence indicating the presence of an IS sequence insertion on mgrB. Simultaneously, a corresponding short read-based rapid alignment and identification process was established. Next, based on all the collected strains, 10X WGS data was simulated using NGS data simulation software to verify the consistency between the read-based alignment and identification process and the initial insertion site analysis results. Finally, WGS sequencing of some clinically cultured isolates also effectively identified the corresponding IS fragment insertion characteristics of colistin-resistant strains.

[0067] In the above embodiments, implementation can be achieved, in whole or in part, through software, hardware, firmware, or any combination thereof. When implemented, in whole or in part, as a computer program product, the computer program product includes one or more computer instructions. When the computer program instructions are loaded or executed on a computer, all or part of the processes or functions described in the embodiments of the present invention are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., coaxial cable, fiber optic, digital subscriber line (DSL) or wireless (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape) or an optical medium.

[0068] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for identifying large sequence insertions in a target genomic region, characterized in that, Alignment using short sequencing reads includes the following steps: S1. Collect the genome and phenotypic information of the target pathogen strain, and perform quality control to remove genomes with poor assembly quality; S2. Select the wild-type target pathogen genome region as the reference sequence, compare the genomes of all strains collected in S1 with the reference sequence, find the insertion site and the IS sequence at the insertion site; cut a short sequence of 10-100 bp from the end of the IS sequence and manually insert it into the corresponding insertion site in the reference sequence to construct a representative sequence that can indicate the insertion of large fragment sequences in the target genome region. S3. Establish a genotyping and identification process for representative sequences of target genomic regions based on the alignment of genome contigs and WGS short sequencing reads, and evaluate and determine the optimal length of the short sequences inserted to construct representative sequences of target genomic regions and the accuracy of the identification process; In S2, the alignment results with the reference sequence are divided into the following five modes: M0: Only one alignment exists, and it covers the entire region of the reference sequence. In this case, no IS sequence is inserted. M1-1: There is only one alignment, which covers the right region of the reference sequence. The alignment breakpoint is the insertion site of the IS sequence, and the unaligned query contig sequence to the left of the alignment breakpoint is the IS sequence. M1-2: There is only one alignment, which covers the left region of the reference sequence. The alignment breakpoint is the insertion site of the IS sequence, and the unaligned query contig sequence to the right of the alignment breakpoint is the IS sequence. M2-1: There are two alignments for the same contig, covering the left L1 and right R2 regions of the reference sequence, respectively. If there is an overlap between the two alignments, there is one insertion site. Further, for the overlapping region, if the number of SNPs and InDels in the L1 alignment is higher than that in the R2 alignment, the insertion site is located at the alignment breakpoint on the left side of the overlapping region; otherwise, the insertion site is located at the alignment breakpoint on the right side of the overlapping region. If there is no overlap between the two alignments, there are two insertion sites. Further, the insertion sites are located at the alignment breakpoint on the right side of the L1 alignment and the alignment breakpoint on the left side of the R2 alignment. The unaligned sequence between L1 and R2 in the query is the IS sequence. M2-2: Two different contigs each have one alignment, resulting in two alignments, covering the left L1 and right R2 regions of the reference sequence, respectively. If there is an overlap between the two alignments, there is one insertion site. Further, regarding the overlap region, if the number of SNPs and InDels in the L1 alignment is higher than that in the R2 alignment, the insertion site is located at the alignment breakpoint on the left side of the overlap region; otherwise, the insertion site is located at the alignment breakpoint on the right side of the overlap region. If there is no overlap between the two alignments, there are two insertion sites. Further, the insertion sites are located at the alignment breakpoint on the right side of the L1 alignment and the alignment breakpoint on the left side of the R2 alignment. The unaligned region on the right side of the query contig corresponding to the L1 alignment and the unaligned region on the left side of the query contig corresponding to the R2 alignment constitute the IS sequence.

2. The method for identifying large sequence insertions in a target genomic region according to claim 1, characterized in that, In S1, genome and phenotypic information of the target pathogen strain are collected from the NCBI NDARO, BV-BRC public database or local hospitals. Genomes with poor assembly quality, such as those with more than 500 genome contigs, genome length exceeding 1.5 times the average length or less than 0.5 times the average length, are removed.

3. The method for identifying large sequence insertions in a target genomic region according to claim 1, characterized in that, In S3, Assembly-based and Read-based alignment and identification procedures were used to identify and analyze the target genomic region sequence types of all strains obtained in S1. The alignment breakpoint analysis results obtained in S2 were used as a reference standard to compare and confirm the identification accuracy of Assembly-based and Read-based procedures and the optimal short sequence insertion length when constructing representative sequences.

4. The method for identifying large sequence insertions in a target genomic region according to claim 3, characterized in that, Assembly-based alignment and identification process: Using Blastn software, the contig sequences of all S1 strain genomes are directly aligned with the wild-type reference sequence of the target genome region and artificially constructed representative sequences. The Hit with the highest score is selected as the final alignment result, which is identified as the final target gene or region sequence type. At the same time, the m0 alignment format file is parsed to obtain SNPs and InDels information on the target gene or region.

5. The method for identifying large sequence insertions in a target genomic region according to claim 3, characterized in that, Read-based alignment and identification process: First, based on all genomes collected by S1, ART software is used to simulate the 10X WGS reads sequence of each strain; then, KMA software is used to align the simulated 10X WGS reads sequence of the strain with the wild-type reference sequence of the target genome region and the artificially constructed representative sequence, and the final results are identified and annotated.

6. The method for identifying large sequence insertions in a target genomic region according to claim 5, characterized in that, The identification and annotation process is as follows: First, for each query read, only the best alignment result is retained. Then, the number of read sequences aligned with each reference sequence is counted and sorted from high to low according to the number of aligned reads. If the number of reads in the second-highest ranking is less than the number of reads in the first-highest ranking, then the sequence type corresponding to the first-highest reference sequence is finally identified. If the number of reads in the second-highest ranking is equal to the number of reads in the first-highest ranking, then the identification resolution is recorded as the family level, that is, at this time, the specific target genomic region sequence type cannot be identified.

7. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored computer program, wherein, when the computer program is executed, it controls the device on which the computer-readable storage medium is located to perform the method as described in any one of claims 1-6.

8. A computer device, characterized in that, The computer device includes a memory, a processor, and a program stored in and executable on the memory, the program being executed by the processor to implement the steps of the method as described in any one of claims 1-6.

Images