A method and system for identifying a tuberculosis gene signature

Abstract

The present application relates to the technical field of gene detection, and discloses a tuberculosis gene feature recognition method and system, the method comprising: obtaining a set of read lengths of mycobacterium tuberculosis to be detected, calibrating base quality based on context entropy value to obtain a quality score matrix; constructing a graph reference index containing a parallel path of a core genome and a repeat sequence polymorphism based on the quality score matrix, performing multi-path alignment to generate an optimal mapping path and a quality value; performing haplotype phase typing in combination with the mapping quality and the score matrix, separating subclone read lengths into different phase blocks, and calculating the frequency of feature sites in each block; combining a sequencing background error model, recognizing low-frequency signals, extracting heterogeneity drug resistance associated features and lineage specificity features; then synthesizing specific capture oligonucleotide probes, and target enrichment containing feature fragments; sequencing the enrichment product, integrating mapping position information, and outputting tuberculosis gene feature recognition results; the present application can improve the accuracy of tuberculosis gene feature recognition.

Description

Technical Field

[0001] This invention relates to the field of gene detection technology, and in particular to a method and system for identifying tuberculosis gene features. Background Technology

[0002] Currently, traditional tuberculosis gene feature identification schemes lack a base quality calibration mechanism based on contextual entropy. They mostly use a linear single reference genome for sequence alignment, which is difficult to adapt to repetitive sequences and polymorphic variation regions. They cannot complete subclonal haplotype phase typing and separation, lack a targeted sequencing background error model to distinguish sequencing noise from real low-frequency variation signals, and there is no probe-targeted enrichment and secondary sequencing co-validation process. As a result, the detection effect of heterogeneous drug resistance and lineage-related low-frequency gene features is poor.

[0003] For example, traditional gene identification methods directly use the raw base quality values ​​output by the sequencer without correcting for systematic sequencing biases caused by upstream and downstream base context; linear reference genome alignment results in high mismatch rates in repetitive polymorphic regions; mixed reads from different subclones cannot be split into phase blocks; random background noise from sequencing is directly identified as variant sites; low-frequency drug resistance mutations are easily masked by noise; and relying solely on the output of a single sequencing run lacks a targeted enrichment verification step, resulting in a high proportion of false positives and false negatives in the final identification results.

[0004] Therefore, existing Mycobacterium tuberculosis gene feature identification technologies have shortcomings such as lack of contextual entropy correction for base quality, poor adaptability of linear reference genomes to polymorphic repetitive regions, inability to phase separate subclonal reads, difficulty in distinguishing sequencing noise from true low-frequency features, and lack of targeted enrichment and collaborative verification processes. This invention adopts an integrated identification architecture that combines contextual entropy base quality calibration, multi-path graph reference genome indexing, weighted consistency association graph haplotype phase typing, sequencing background error model to filter noise, probe set targeted enrichment, and secondary sequencing collaborative verification. Summary of the Invention

[0005] This invention provides a method and system for identifying tuberculosis gene features to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides a method for identifying tuberculosis gene features, comprising:

[0007] S1. Obtain the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and perform weight calibration on the base quality in the sequencing read set based on the entropy value of the context environment to obtain the quality score matrix.

[0008] S2. Based on the quality score matrix, construct a graph reference genome index containing parallel paths of polymorphism in the core genome and repetitive sequence regions. Perform multi-path alignment between the sequencing read set and the graph reference genome index to generate the optimal mapping path and its mapping quality value.

[0009] S3. Based on the mapping quality value and quality score matrix, perform haplotype phase typing on the sequencing read set, separate different subclone reads into different phase blocks, and calculate the frequency distribution of characteristic sites in each phase block;

[0010] S4. Based on the frequency distribution and combined with the pre-constructed sequencing background error model, identify low-frequency feature signals, and extract heterogeneous drug resistance association features and lineage-specific features from the low-frequency feature signals;

[0011] S5. Based on the heterogeneity of drug resistance association characteristics and lineage specificity characteristics, synthesize a set of specific capture oligonucleotide probes for hybridization with target fragments carrying the aforementioned characteristics to achieve targeted enrichment;

[0012] S6. Sequencing of the enriched target nucleic acid fragments, integration of mapping location information using the validation sequencing results, and outputting tuberculosis gene feature identification results.

[0013] Preferably, the step of obtaining the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and weighting and calibrating the base quality in the sequencing read set based on the entropy value of the context environment to obtain a quality score matrix includes:

[0014] Based on each base site of each read in the sequencing read set, the base type combination of the upstream and downstream neighboring bases of the base site is extracted as contextual information, and the initial base quality value of the base site is obtained.

[0015] Based on the context information, the base identification error frequency of all base sites with the same context information in the sequencing read set is statistically analyzed, and the empirical error rate of the context is calculated based on the base identification error frequency.

[0016] The entropy weighting factor of the context is determined based on the deviation between the empirical error rate and the expected error rate characterized by the initial base quality value.

[0017] The initial base quality value of the base site is weighted and corrected according to the entropy weight factor corresponding to the context of the base site to obtain the calibrated base quality value. The calibrated base quality values ​​are then summarized to construct a quality score matrix.

[0018] Preferably, the step of constructing a graph reference genome index containing parallel paths for polymorphisms in the core genome and repetitive sequence regions based on the quality score matrix, and performing multi-path alignment between the sequencing read set and the graph reference genome index to generate the optimal mapping path and its mapping quality value includes:

[0019] Based on the distribution of base quality values ​​in the quality score matrix across the genome, the sequencing read set is divided into high-confidence alignment regions and low-confidence alignment regions.

[0020] For the high-confidence alignment regions, a single linear reference path is constructed; for the low-confidence alignment regions, known polymorphic sequence variants are obtained and multiple parallel reference paths are constructed. All paths are then integrated into a graph reference genome index.

[0021] Each read in the sequencing read set is aligned with each reference path in the graph reference genome index to obtain an initial alignment score for each read on the path.

[0022] Based on the quality values ​​of the read length at each alignment position in the quality score matrix, the initial alignment score is weighted and corrected to obtain the weighted alignment total score for each path;

[0023] The path with the highest weighted total score is selected as the optimal mapping path, and the mapping quality value of the optimal mapping path is generated based on the difference in weighted total score between the optimal path and the second-best path.

[0024] Preferably, for the high-confidence alignment region, a single linear reference path is constructed; for the low-confidence alignment region, known polymorphic sequence variants are obtained and multiple parallel reference paths are constructed; all paths are integrated into a graph reference genome index, including:

[0025] Based on the quality value of each read locus in the quality score matrix, the average quality value of the reads covered by each genomic region is calculated. Regions with an average quality value higher than a preset reliability threshold are classified as high-confidence alignment regions, and regions with an average quality value lower than the reliability threshold are classified as low-confidence alignment regions.

[0026] For the high-confidence alignment regions, a single linear reference path is constructed using the reference genome sequence of the corresponding region as a template;

[0027] For the low-confidence alignment region, all known polymorphic sequence variants in the corresponding region are obtained, and parallel reference paths are constructed for each polymorphic sequence variant. Each parallel path establishes a connection node with the adjacent region of the single linear reference path at the end of the path, generating a preliminary graph index containing a branch structure.

[0028] The single linear reference path and the parallel reference path are arranged in genomic coordinate order. The paths are connected end to end through the connection nodes between each path to generate a connected graph structure that covers the core genome and repetitive sequence regions and includes polymorphic branch paths. The connected graph structure is used as the graph reference genome index.

[0029] Preferably, the step of performing haplotype phase typing on the sequencing read set based on the mapped quality value and quality score matrix, separating different subclonal reads into different phase blocks, and calculating the frequency distribution of characteristic sites in each phase block includes:

[0030] Based on the mapped quality value, the effective read lengths for genotyping are selected, and the calibrated quality value is used as the base confidence weight to generate a weighted base sequence for each effective read length.

[0031] Based on the ratio of the sum of the weighted product of two reads with identical bases at the common site to the total sum of the weighted product, a weighted consistency association graph is constructed, and the connected regions in the graph are defined as the initial phase blocks.

[0032] The base divergence of candidate sites within the initial phase block is detected. Based on the cumulative weight of the base divergence, low-weight erroneous bases are discarded, and the consensus sequence of the corresponding block is reconstructed based on the remaining read length.

[0033] Based on the consensus sequence, the read length assignment is corrected, and the read lengths that still cannot be assigned to the existing block after correction are formed into new phase blocks separately, thus obtaining the final phase block;

[0034] The cumulative confidence weight of each base type at each feature site within the final phase block is calculated, and the proportion of the cumulative weight of each base type to the total cumulative weight of the corresponding site is used as the frequency distribution of the corresponding feature site.

[0035] Preferably, the step of identifying low-frequency feature signals based on frequency distribution and in conjunction with a pre-constructed sequencing background error model, and extracting heterogeneous drug resistance association features and lineage-specific features from the low-frequency feature signals, includes:

[0036] The frequency distribution of feature sites in the phase block is compared one by one with the background error frequency threshold of the corresponding site context in the pre-constructed sequencing background error model. Base types whose frequency values ​​exceed the preset threshold are selected as candidate low-frequency signal bases, and the excess amount of the candidate low-frequency signal bases is recorded as the signal intensity value.

[0037] Based on the number of read lengths supported by candidate low-frequency signal bases in the corresponding phase blocks and the signal intensity value, the micro-enrichment of candidate low-frequency signal bases in the subclonal population is calculated.

[0038] Candidate low-frequency signal bases with micro-enrichment levels higher than a preset standard are identified as true low-frequency feature signals, thus obtaining a set of true low-frequency feature signals.

[0039] Based on the genomic coordinates of each point in the real low-frequency feature signal set, signals located in known drug resistance gene regions are classified as heterogeneous drug resistance associated features, and signals located in phylogenetic marker gene regions are classified as lineage-specific features.

[0040] Preferably, the pre-constructed sequencing background error model includes:

[0041] The conserved region sequence of the Mycobacterium tuberculosis genome was obtained as the background training base, and the reads that were consistent with the background training base were extracted from historical sequencing data as the training read set.

[0042] The occurrence frequency of each type of base that is inconsistent with the reference base of the training read set is counted, and the occurrence frequency is divided by the total number of times the corresponding site is covered by the training read to obtain the empirical error rate of each base substitution type at the corresponding site.

[0043] Using the upstream and downstream adjacent base sequences of each site as context markers, sites with the same context markers are grouped into the same context category;

[0044] The empirical error rates of each base substitution type at all sites within the same environmental category are summarized, and the summarized results are used as the background error frequency threshold for the corresponding environmental category. The correspondence between environmental categories and background error frequency thresholds is established, and a sequencing background error model is constructed.

[0045] Preferably, the synthesis of a specific capture oligonucleotide probe set based on heterogeneous drug resistance association characteristics and lineage-specific characteristics, for hybridization with target fragments carrying the aforementioned characteristics to achieve targeted enrichment, includes:

[0046] Based on the sequence composition of the heterogeneous drug resistance association features and the lineage-specific features, high-affinity capture probes are designed for the drug resistance association features and discriminative capture probes are designed for the lineage-specific features, and the two types of probes are mixed to form an initial probe set;

[0047] The binding free energy of each probe in the initial probe set with the corresponding flanking sequence of the characteristic site is evaluated. The probe with the lowest binding free energy is selected as the optimal capture probe for the characteristic site and replaces the corresponding probe in the initial probe set to generate a specific capture oligonucleotide probe set.

[0048] The specific capture oligonucleotide probe set is liquid-phase hybridized with the nucleic acid fragment of the sample to be identified, so that the target fragment carrying heterogeneous drug resistance association characteristics or lineage specific characteristics forms a hybridization complex with the corresponding probe. The hybridization complex is separated from the hybridization system to obtain the targeted enriched nucleic acid fragment.

[0049] Preferably, the step of sequencing the enriched target nucleic acid fragments, integrating the mapping location information using the verification sequencing results, and outputting tuberculosis gene feature identification results includes:

[0050] The targeted enriched nucleic acid fragments were sequenced to obtain a set of targeted sequencing reads. The set of targeted sequencing reads was then re-aligned to the reference genome index to obtain the mapping position of each targeted read on the genome.

[0051] The mapping position is compared with the genomic coordinates of the heterogeneous drug resistance association feature and the lineage-specific feature. The reads that overlap with any feature coordinate are selected as effective feature reads. The number of times each feature site is covered by the effective feature reads is counted as the verification confidence of the corresponding site.

[0052] The verification confidence of each feature site is compared with the previously obtained frequency distribution and signal intensity value. Feature sites that meet the consistency criteria are marked as confirmed features. All confirmed features and corresponding classification information are summarized to generate tuberculosis gene feature identification results.

[0053] To address the above problems, the present invention also provides a tuberculosis gene feature recognition system, the system comprising:

[0054] The entropy weight quality calibration module is used to obtain the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and to perform weight calibration on the base quality in the sequencing read set based on the entropy value of the context environment to obtain the quality score matrix.

[0055] The graph alignment and mapping module is used to construct a graph reference genome index containing parallel paths of polymorphism in the core genome and repetitive sequence regions based on the quality score matrix, and to perform multi-path alignment between the sequencing read set and the graph reference genome index to generate the optimal mapping path and its mapping quality value.

[0056] The haplotype typing module is used to perform haplotype phase typing on the sequencing read set according to the mapped quality value and quality score matrix, separate different subclone reads into different phase blocks, and calculate the frequency distribution of characteristic sites in each phase block;

[0057] The low-frequency feature extraction module is used to identify low-frequency feature signals based on frequency distribution and in conjunction with a pre-constructed sequencing background error model, and to extract heterogeneous drug resistance association features and lineage-specific features from the low-frequency feature signals.

[0058] The probe-targeted enrichment module is used to synthesize a set of specific capture oligonucleotide probes based on heterogeneous drug resistance association characteristics and lineage-specific characteristics, which are used to hybridize with target fragments carrying the aforementioned characteristics to achieve targeted enrichment.

[0059] The sequencing feature recognition module is used to sequence the enriched target nucleic acid fragments, integrate the mapping position information using the verification sequencing results, and output the tuberculosis gene feature recognition results.

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

[0061] 1. This invention extracts upstream and downstream contextual information of bases, calculates entropy weighting factors based on empirical error rates to complete weighted correction of base quality, and constructs a quality score matrix; it distinguishes genomic confidence regions based on base quality distribution, builds a graph reference genome index containing polymorphic parallel branch paths, performs weighted multipath sequence alignment to output the optimal mapping path and mapping quality, realizes fine correction of sequencing base quality and high-fit sequence alignment of polymorphic and repetitive genomic regions, and improves the data reliability of the sequencing read length mapping process.

[0062] 2. This invention relies on mapping quality and calibrated base quality to complete weighted haplotype phase typing, realize subclonal read phase block division and statistically analyze site frequency distribution; combined with sequencing background error model to identify real low-frequency signals, separate heterogeneous drug resistance and lineage-specific gene features, and obtain capture probe set through probe binding free energy optimization to complete nucleic acid fragment targeted enrichment. Then, the enriched sequencing data is used for multi-dimensional consistency verification, and the confirmed tuberculosis gene feature identification results are output, which fully realizes the analysis and secondary verification of low-frequency variant features and improves the confidence level of gene feature identification results. Attached Figure Description

[0063] Figure 1 This is a flowchart illustrating a tuberculosis gene feature identification method according to an embodiment of the present invention;

[0064] Figure 2 This is a functional block diagram of a tuberculosis gene feature recognition system provided in an embodiment of the present invention;

[0065] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0066] It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0067] This application provides a method for identifying tuberculosis gene features. The executing entity of this tuberculosis gene feature identification method includes, but is not limited to, at least one of the following electronic devices that can be configured to execute the method provided in this application embodiment: a server, a terminal, etc. In other words, the tuberculosis gene feature identification method can be executed by software or hardware installed on a terminal device or a server device. The server includes, but is not limited to, a single server, a server cluster, a cloud server, or a cloud server cluster. The server can be an independent server or a cloud server that provides basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communication, middleware services, domain name services, security services, content delivery networks (CDNs), and big data and artificial intelligence platforms.

[0068] Reference Figure 1 The diagram shown is a flowchart illustrating a tuberculosis gene feature identification method according to an embodiment of the present invention. In this embodiment, the tuberculosis gene feature identification method includes:

[0069] S1. Obtain the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and perform weight calibration on the base quality in the sequencing read set based on the entropy value of the context environment to obtain the quality score matrix.

[0070] In this embodiment of the invention, the step of obtaining the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and weighting and calibrating the base quality in the sequencing read set based on the entropy value of the context environment to obtain a quality score matrix includes:

[0071] Based on each base site of each read in the sequencing read set, the base type combination of the upstream and downstream neighboring bases of the base site is extracted as contextual information, and the initial base quality value of the base site is obtained.

[0072] Based on the context information, the base identification error frequency of all base sites with the same context information in the sequencing read set is statistically analyzed, and the empirical error rate of the context is calculated based on the base identification error frequency.

[0073] The entropy weighting factor of the context is determined based on the deviation between the empirical error rate and the expected error rate characterized by the initial base quality value.

[0074] The initial base quality value of the base site is weighted and corrected according to the entropy weight factor corresponding to the context of the base site to obtain the calibrated base quality value. The calibrated base quality values ​​are then summarized to construct a quality score matrix.

[0075] Specifically, the process iterates through each base site of every sequencing read in the Mycobacterium tuberculosis sequencing read set of the sample to be identified. Centered on a single base site, the base types of the one base immediately upstream and one base immediately downstream of that site are fixed and extracted. The triadic combination formed by the central base and the upstream and downstream bases is used as the unique contextual information of that base site. This range is determined based on statistical analysis of the systematic bias of the Illumina sequencing platform. One bp upstream and downstream is sufficient to capture the influence of major neighboring bases on the mismatch rate. At the same time, the numerical initial base quality value of that base site is directly retrieved from the original output of the sequencing platform to complete the basic data extraction of a single base site. For example, if a base site is adenine A, the upstream base is guanine G, and the downstream base is thymine T, then the contextual information of that site is GAT, and the original initial base quality value corresponding to that site is obtained simultaneously. .

[0076] The contextual information of all base sites in the sequencing read set is classified and grouped. All base sites with completely identical base type combinations (i.e., triplet sequences with completely identical contextual information) are grouped into the same data group, and the total number of all base sites in each group is counted one by one. The number of incorrectly identified base sites, M, is determined by the misidentification of base types in the sequencing genome, where the base type does not match the corresponding base type in the Mycobacterium tuberculosis standard reference genome. This is determined using the formula... Calculate the empirical error rate for each context environment, and complete the statistical calculation of the error rate for all context environments. For example, if 1000 base sites with the context environment of GAT are collected, and 12 of these sites have base recognition errors, then the empirical error rate for that context environment is calculated. .

[0077] For each context, the frequency of occurrence of each base type (A, T, C, G) in that context is statistically analyzed. According to Shannon's entropy formula Calculate the entropy value of this context, which measures the degree of disorder in the base composition within that environment; simultaneously, convert it according to the sequencing quality standard conversion formula. Calculate the expected error rate corresponding to the initial quality value of a single base site. Then, using the weighting factor formula ,in The entropy weighting factor for that context is calculated based on the maximum entropy value across all contexts (used for normalization). This factor comprehensively reflects the degree of deviation of the empirical error rate from the expected error rate, as well as the uncertainty of the base distribution in the environment itself. Compared with a simple difference, it can more comprehensively correct for systematic biases.

[0078] For each base site in the sequencing read set, a unique entropy weighting factor corresponding to its specific context is matched, and then a calibration formula is used. The initial base quality values ​​are weighted and corrected, where... This is the calibrated base quality value. This represents the original initial base mass value. To determine the entropy weighting factor corresponding to the context, after calibrating and calculating the quality value of all bases at each site, the calibrated base quality values ​​of all sites are summarized in the original site arrangement order of the sequencing read length, and then organized into a matrix format with row and column correspondence and one-to-one site matching, finally generating a complete quality score matrix.

[0079] This step extracts the upstream and downstream context of bases and combines it with empirical error rates to obtain an entropy weighting factor, thus completing the base quality weighted calibration. This corrects the systematic bias of the sequencer's original quality values ​​caused by the base environment, improves the reliability of input data for subsequent alignment and genotyping, reduces noise interference, lays a high-quality data foundation for subsequent graph genome alignment and low-frequency variant identification, and reduces false positives and false negatives.

[0080] S2. Based on the quality score matrix, construct a graph reference genome index containing parallel paths of polymorphism in the core genome and repetitive sequence regions. Perform multi-path alignment between the sequencing read set and the graph reference genome index to generate the optimal mapping path and its mapping quality value.

[0081] In this embodiment of the invention, the step of constructing a graph reference genome index containing parallel paths of polymorphisms in the core genome and repetitive sequence regions based on the quality score matrix, and performing multi-path alignment between the sequencing read set and the graph reference genome index to generate the optimal mapping path and its mapping quality value includes:

[0082] Based on the distribution of base quality values ​​in the quality score matrix across the genome, the sequencing read set is divided into high-confidence alignment regions and low-confidence alignment regions.

[0083] For the high-confidence alignment regions, a single linear reference path is constructed; for the low-confidence alignment regions, known polymorphic sequence variants are obtained and multiple parallel reference paths are constructed. All paths are then integrated into a graph reference genome index.

[0084] Each read in the sequencing read set is aligned with each reference path in the graph reference genome index to obtain an initial alignment score for each read on the path.

[0085] Based on the quality values ​​of the read length at each alignment position in the quality score matrix, the initial alignment score is weighted and corrected to obtain the weighted alignment total score for each path;

[0086] The path with the highest weighted total score is selected as the optimal mapping path, and the mapping quality value of the optimal mapping path is generated based on the difference in weighted total score between the optimal path and the second-best path.

[0087] For the high-confidence alignment regions, a single linear reference path is constructed; for the low-confidence alignment regions, known polymorphic sequence variants are obtained and multiple parallel reference paths are constructed. All paths are integrated into a graph reference genome index, including:

[0088] Based on the quality value of each read locus in the quality score matrix, the average quality value of the reads covered by each genomic region is calculated. Regions with an average quality value higher than a preset reliability threshold are classified as high-confidence alignment regions, and regions with an average quality value lower than the reliability threshold are classified as low-confidence alignment regions.

[0089] For the high-confidence alignment regions, a single linear reference path is constructed using the reference genome sequence of the corresponding region as a template;

[0090] For the low-confidence alignment region, all known polymorphic sequence variants in the corresponding region are obtained, and parallel reference paths are constructed for each polymorphic sequence variant. Each parallel path establishes a connection node with the adjacent region of the single linear reference path at the end of the path, generating a preliminary graph index containing a branch structure.

[0091] The single linear reference path and the parallel reference path are arranged in genomic coordinate order. The paths are connected end to end through the connection nodes between each path to generate a connected graph structure that covers the core genome and repetitive sequence regions and includes polymorphic branch paths. The connected graph structure is used as the graph reference genome index.

[0092] Specifically, the quality values ​​of all bases in each coordinate region of the genome within the quality score matrix are aggregated and calculated one by one, using the formula for calculating the regional average quality value. In the formula For the genome The average quality value of each region This represents the total number of base sites covered by this genomic region. For the first in this region The calibrated base quality values ​​corresponding to each base site in the quality score matrix were used to predetermine the reliability threshold as follows: Sequencing data from 10 publicly available Mycobacterium tuberculosis standard strains (SRR numbers XXXXX~XXXXX, sourced from the NCBI SRA database) were selected. The average quality value of each genomic region was calculated. The 5th percentile of the average quality value of the regions verified as correctly aligned by known variant sites across all samples was used as the threshold. This threshold was statistically found to be 20, therefore the predefined reliability threshold was fixed at 20. Genomic regions with an average quality value greater than 20 were designated as high-confidence alignment regions, and genomic regions with an average quality value less than 20 were designated as low-confidence alignment regions.

[0093] The complete base sequence of the high-confidence alignment region corresponding to the standard reference genome of Mycobacterium tuberculosis is retrieved. Using this standard base sequence as the only template, a single linear reference path with no branches, no variations, and continuous base sites is generated according to the fixed arrangement order from the 5' end to the 3' end of the genome coordinates. The base arrangement, length, and coordinate interval of this path are completely consistent with the corresponding region of the standard reference genome.

[0094] The public Mycobacterium tuberculosis polymorphism database was retrieved. This database integrates whole-genome sequencing variant data from clinical isolates in public databases such as NCBI SRA, EBI-ENA, and GISAID, as well as verified polymorphic sites included in professional tuberculosis databases such as TBDB and PATRIC. Specifically, it was obtained by periodically downloading updates from a public FTP site or searching online via API. For low-confidence alignment regions, all reported polymorphic sequence variants verified by Sanger sequencing within these regions were retrieved from the database, including SNPs, short InDels, and copy number variations of repetitive sequences. For each independent polymorphic sequence variant retrieved, an independent parallel reference path was generated. The genomic start and end coordinates of all parallel reference paths were completely consistent with the corresponding low-confidence alignment region. Fixed connection nodes were set at both ends of each parallel reference path to the endpoints of adjacent single linear reference paths, so that each parallel branch path was attached to the main body of the linear reference path, generating a preliminary graph index with a multi-branch topology.

[0095] Following the ascending global coordinate order of the Mycobacterium tuberculosis standard genome, all single linear reference paths and multiple parallel reference paths are ordered by position. Through the fixed connection nodes set in the early stage, adjacent linear paths and branch paths are precisely connected end to end, so that the scattered independent paths form a complete and connected topological network structure. This structure completely covers the core genome region and repetitive sequence polymorphism region of Mycobacterium tuberculosis. This complete connected graph structure is directly defined as the graph reference genome index.

[0096] Each sequencing read within the sequencing read set is sequentially compared globally with all single linear reference paths and polymorphic parallel reference paths contained in the graph reference genome index. A fixed alignment scoring rule is used, which includes base matching score (match +1), mismatch penalty score (mismatch -5), and insertion / deletion penalty score (insertion / deletion -2), to obtain the initial alignment score for each sequencing read on each reference path.

[0097] Retrieve the calibrated base quality values ​​corresponding to each alignment site of the current sequencing read length from the quality score matrix, and calculate using a weighted correction formula. In the formula The weighted comparison score for a single reference path. This represents the total number of base sites involved in the alignment of the current read length. For the first The calibrated quality value of each base site. For the first The initial unit alignment score of each base site is calculated by weighting and accumulating the scores site by site to obtain the independent weighted alignment total score of each sequencing read across all reference paths.

[0098] The weighted ratio scores of all paths corresponding to a single sequencing read are sorted by value, and the reference path corresponding to the largest weighted ratio score after sorting is determined as the optimal mapping path for that sequencing read.

[0099] Using the formula for calculating the quality value of the mapping In the formula The mapping quality value is the value for the optimal mapping path. The weighted alignment score is the total score for the optimal mapping path. The weighted comparison total score of the suboptimal mapping path is calculated by performing a logarithmic transformation based on the difference in total scores between the two paths, generating a unique mapping quality value for the corresponding optimal mapping path. For example, if the weighted total score of the optimal path is 520 and the weighted total score of the suboptimal path is 480, the difference is 40. Substituting these values ​​into the formula will yield the corresponding mapping quality value.

[0100] This step utilizes a quality score matrix to divide confidence regions, constructs a graph reference genome index with polymorphic parallel paths, and performs weighted multi-path alignment. This improves upon the shortcomings of traditional linear reference genomes in poor adaptation to repetitive and polymorphic regions, outputting reliable mapping quality values. This provides an accurate alignment basis for subsequent haplotype phase typing and locus frequency calculations, reducing alignment errors and false positives in subsequent analyses.

[0101] S3. Based on the mapping quality value and quality score matrix, perform haplotype phase typing on the sequencing read set, separate different subclone reads into different phase blocks, and calculate the frequency distribution of characteristic sites in each phase block;

[0102] In this embodiment of the invention, the step of performing haplotype phase typing on the sequencing read set based on the mapped quality value and quality score matrix, separating different subclonal reads into different phase blocks, and calculating the frequency distribution of characteristic sites in each phase block includes:

[0103] Based on the mapped quality value, the effective read lengths for genotyping are selected, and the calibrated quality value is used as the base confidence weight to generate a weighted base sequence for each effective read length.

[0104] Based on the ratio of the sum of the weighted product of two reads with identical bases at the common site to the total sum of the weighted product, a weighted consistency association graph is constructed, and the connected regions in the graph are defined as the initial phase blocks.

[0105] The base divergence of candidate sites within the initial phase block is detected. Based on the cumulative weight of the base divergence, low-weight erroneous bases are discarded, and the consensus sequence of the corresponding block is reconstructed based on the remaining read length.

[0106] Based on the consensus sequence, the read length assignment is corrected, and the read lengths that still cannot be assigned to the existing block after correction are formed into new phase blocks separately, thus obtaining the final phase block;

[0107] The cumulative confidence weight of each base type at each feature site within the final phase block is calculated, and the proportion of the cumulative weight of each base type to the total cumulative weight of the corresponding site is used as the frequency distribution of the corresponding feature site.

[0108] Specifically, a preset fixed mapping quality value screening threshold of 30 was used to screen the optimal mapping path mapping quality values ​​for all sequencing reads. This threshold was determined as follows: using simulated data of two Mycobacterium tuberculosis subclones with known mixing ratios (1:1, 1:5, 1:10), the genotyping accuracy under different mapping quality thresholds was tested. When the threshold was set to 30, the accuracy stabilized above 95%, so it was fixed at 30. Only sequencing reads with mapping quality values ​​greater than or equal to 30 were retained as valid reads for haplotype phase genotyping. The calibrated base quality value of each base site corresponding to the valid read was retrieved from the quality score matrix, and this quality value was used as the base confidence weight for the corresponding base. By combining the base type with the corresponding confidence weight, weight parameters were bound to each valid read site by site, generating a weighted base sequence carrying site confidence parameters.

[0109] Common alignment sites are extracted from any two valid reads. The set of matching sites where the base types of the two reads are completely identical at these common sites is calculated, along with the set of all common sites. For each matching site, the confidence weights of the corresponding bases for the two reads are multiplied and summed to obtain the cumulative product of the consistency weights. Similarly, for each of all common sites, the confidence weights of the corresponding bases for the two reads are multiplied and summed to obtain the cumulative product of the total weights. The ratio of the cumulative product of the consistency weights to the cumulative product of the total weights is used as the weighted consistency correlation between the two reads. All valid read lengths are treated as network nodes, and the calculated weighted consistency correlation degree is used as the weight of the node association edge to construct a complete weighted consistency correlation graph. The specific construction process is as follows: Initialize an undirected graph G(V,E), where V is all valid read length nodes. For each pair of read lengths (u,v), if its weighted consistency correlation degree is greater than or equal to 0.6 (this threshold is obtained by simulating mixed samples for optimization), then add an edge between u and v with a weight equal to the correlation degree. Perform a depth-first search to traverse the connected components of the weighted consistency correlation graph, and delineate all interconnected node sets in the graph as independent initial phase blocks. Each connected component corresponds to a unique initial phase block.

[0110] The algorithm iterates through all read-covered variant sites within each initial phase block. Sites with two or more base types are defined as candidate divergence sites. For each candidate divergence site, the confidence weights of all reads within the block are aggregated, and the weights of different base types are summed to obtain the cumulative divergence weight for each base. A fixed error weight threshold of 5 is set. This threshold is determined as follows: for sequencing data of known homozygous strains without mixing, the cumulative weight distribution of non-reference bases within the same block is statistically analyzed. The 99th percentile is lower than 5, so 5 is used as the minimum weight threshold to distinguish between sequencing errors and true variants. Base types with cumulative weights below this threshold within a site are identified as erroneous bases caused by sequencing errors and are discarded directly. Valid base types with cumulative weights above the threshold are retained. After removing read site data corresponding to all erroneous bases, based on the remaining valid base data and corresponding complete read sequences of each candidate site, the dominant base type (the one with the largest weight sum) is statistically analyzed for each site, and these are spliced ​​to generate a unique block consensus sequence for the corresponding initial phase block.

[0111] The weighted base sequences of all valid reads are globally matched and compared with the consensus sequences reconstructed from each initial phase block. Similarity is calculated using edit distance (with a penalty of 1 for replacements, insertions, and deletions), and reads are assigned to the corresponding phase block with the smallest edit distance. If the minimum edit distance is greater than read length × 0.1 (i.e., similarity is less than 90%), the read is considered an anomalous read that cannot be assigned to an existing block. All anomalous reads are collectively collected and used to form independent new phase blocks (if the number of anomalous reads is less than 3, they are discarded to exclude low-quality data). After all blocks have undergone read assignment correction and new blocks have been added, a final phase block with a complete structure and accurate read classification is formed.

[0112] The variant sites used to characterize subclonal differences within all final phase blocks are identified as feature sites. Each feature site in each final phase block is traversed, and the confidence weights of all read lengths corresponding to each base type (A, T, C, G) within that site are calculated. The cumulative confidence weights of each base type are accumulated to obtain the individual cumulative confidence weights. The total cumulative weight of the site is obtained by summing the cumulative confidence weights of all base types at that site. The dimensionless ratio obtained by dividing the cumulative confidence weight of a single base type by the total cumulative weight of the site is the frequency distribution result of the base type corresponding to that feature site.

[0113] This step utilizes mapping quality and calibration base quality to perform weighted haplotype phase typing. Phase blocks are divided using a weighted consistency correlation graph, enabling effective separation of subclonal reads and statistical analysis of site frequency distribution. This process effectively eliminates sequencing error interference, accurately analyzes subclonal composition, and provides reliable frequency input for subsequent sequencing background error models to identify low-frequency drug resistance and lineage characteristics, thus reducing errors in variant detection.

[0114] S4. Based on the frequency distribution and combined with the pre-constructed sequencing background error model, identify low-frequency feature signals, and extract heterogeneous drug resistance association features and lineage-specific features from the low-frequency feature signals;

[0115] In this embodiment of the invention, the step of identifying low-frequency feature signals based on frequency distribution and in conjunction with a pre-constructed sequencing background error model, and extracting heterogeneous drug resistance association features and lineage-specific features from the low-frequency feature signals, includes:

[0116] The frequency distribution of feature sites in the phase block is compared one by one with the background error frequency threshold of the corresponding site context in the pre-constructed sequencing background error model. Base types whose frequency values ​​exceed the preset threshold are selected as candidate low-frequency signal bases, and the excess amount of the candidate low-frequency signal bases is recorded as the signal intensity value.

[0117] Based on the number of read lengths supported by candidate low-frequency signal bases in the corresponding phase blocks and the signal intensity value, the micro-enrichment of candidate low-frequency signal bases in the subclonal population is calculated.

[0118] Candidate low-frequency signal bases with micro-enrichment levels higher than a preset standard are identified as true low-frequency feature signals, thus obtaining a set of true low-frequency feature signals.

[0119] Based on the genomic coordinates of each point in the real low-frequency feature signal set, signals located in known drug resistance gene regions are classified as heterogeneous drug resistance associated features, and signals located in phylogenetic marker gene regions are classified as lineage-specific features.

[0120] The pre-constructed sequencing background error model includes:

[0121] The conserved region sequence of the Mycobacterium tuberculosis genome was obtained as the background training base, and the reads that were consistent with the background training base were extracted from historical sequencing data as the training read set.

[0122] The occurrence frequency of each type of base that is inconsistent with the reference base of the training read set is counted, and the occurrence frequency is divided by the total number of times the corresponding site is covered by the training read to obtain the empirical error rate of each base substitution type at the corresponding site.

[0123] Using the upstream and downstream adjacent base sequences of each site as context markers, sites with the same context markers are grouped into the same context category;

[0124] The empirical error rates of each base substitution type at all sites within the same environmental category are summarized, and the summarized results are used as the background error frequency threshold for the corresponding environmental category. The correspondence between environmental categories and background error frequency thresholds is established, and a sequencing background error model is constructed.

[0125] Specifically, the pre-constructed sequencing background error model is built as follows: Coding regions from the Mycobacterium tuberculosis standard reference genome that have undergone multi-species conservation verification, have no known polymorphisms, and have no repetitive sequences (e.g., conserved domains of the rpoB and katG genes) and non-coding regions (e.g., spacers outside the RRDR region) are selected, totaling 50 regions with a length of approximately 10kb, as the background training base. At least 50 sets of Mycobacterium tuberculosis sequencing data (without drug resistance-related variants) from healthy individuals or standard strains are downloaded from the NCBI SRA database. All reads that can be completely and accurately aligned to the background training base sequence are extracted and aggregated into a training read set. All reads covering each genomic locus of the background training base are traversed within the training read set. The reference base for each locus is fixed as the base standard base. The number of occurrences of all substitution bases inconsistent with the reference base type within each locus is counted individually. The empirical error rate for a single base substitution type = the total number of occurrences of that substitution base at the current locus / the total number of times the current locus is covered by the training reads. Centered on each site, extract the upstream 1bp and downstream 1bp to form a triplet context environment identifier. Sites with the same context environment identifier are grouped into the same environment category. The empirical error rates of each base substitution type of all sites within the same environment category are summarized, and the 95th percentile is taken as the background error frequency threshold of the environment category. The correspondence between the environment category and the background error frequency threshold is established, thus forming the sequencing background error model.

[0126] During application, the base frequency distribution results of all feature sites within the final phase block are retrieved. The triadic context environment identifier corresponding to each feature site is matched, and the background error frequency threshold for the corresponding environment category is retrieved from the pre-constructed sequencing background error model. The actual frequency distribution value of each base type at the feature site is compared one by one with the corresponding background error frequency threshold. Base types whose actual frequency distribution value is strictly greater than the background error frequency threshold are uniformly identified as candidate low-frequency signal bases. Simultaneously, the difference between the actual base frequency distribution value and the corresponding background error frequency threshold is fixedly recorded as the signal intensity value of that candidate low-frequency signal base.

[0127] The number of reads supported by the candidate low-frequency signal base in the corresponding phase block (i.e., the number of effective reads that cover the site and whose bases are consistent with the candidate base). and signal strength value Based on this, the micro-enrichment degree calculation formula is used. ,in The total coverage reads for this site are calculated, with a constant of 0.01 added to the denominator to prevent division by zero. The micro-enrichment of candidate low-frequency signal bases in the subclonal population is then obtained. Candidate low-frequency signal bases with a micro-enrichment higher than a preset standard of 0.05 (determined by the 99th percentile of the enrichment distribution of negative control samples) are identified as true low-frequency feature signals, thus obtaining the true low-frequency feature signal set.

[0128] The precise genomic coordinates corresponding to each signal in the real low-frequency feature signal set were retrieved and compared with the coordinate set of pre-annotated Mycobacterium tuberculosis drug resistance gene regions (derived from the TBDB database and the WHO drug resistance mutation catalog, including the exons and regulatory regions of all known drug resistance genes such as rpoB, katG, inhA, embB, and pncA, with coordinates based on NC_000962.3). Signals whose genomic coordinates completely fall within the drug resistance gene annotation region were uniformly classified as heterogeneous drug resistance association features; signals whose coordinates fall within phylogenetic marker gene regions, such as housekeeping genes like gyrB, tuf, and rrs, as well as SNP barcode sites, with coordinates also based on NC_000962.3, were classified as phylogenetic specific features.

[0129] This step involves constructing a sequencing background error model adapted to Mycobacterium tuberculosis, distinguishing sequencing noise from true variants by combining site context, obtaining true low-frequency signals through micro-enrichment screening, and differentiating heterogeneous drug resistance and lineage-specific characteristics. This reduces false positives caused by noise, provides reliable targets for subsequent probe design and targeted enrichment, reduces false positives in the overall identification process, and improves the accuracy of tuberculosis gene feature detection.

[0130] S5. Based on the heterogeneity of drug resistance association characteristics and lineage specificity characteristics, synthesize a set of specific capture oligonucleotide probes for hybridization with target fragments carrying the aforementioned characteristics to achieve targeted enrichment;

[0131] In this embodiment of the invention, the synthesis of a specific capture oligonucleotide probe set based on heterogeneous drug resistance association characteristics and lineage-specific characteristics, for hybridization with target fragments carrying the aforementioned characteristics to achieve targeted enrichment, includes:

[0132] Based on the sequence composition of the heterogeneous drug resistance association features and the lineage-specific features, high-affinity capture probes are designed for the drug resistance association features and discriminative capture probes are designed for the lineage-specific features, and the two types of probes are mixed to form an initial probe set;

[0133] The binding free energy of each probe in the initial probe set with the corresponding flanking sequence of the characteristic site is evaluated. The probe with the lowest binding free energy is selected as the optimal capture probe for the characteristic site and replaces the corresponding probe in the initial probe set to generate a specific capture oligonucleotide probe set.

[0134] The specific capture oligonucleotide probe set is liquid-phase hybridized with the nucleic acid fragment of the sample to be identified, so that the target fragment carrying heterogeneous drug resistance association characteristics or lineage specific characteristics forms a hybridization complex with the corresponding probe. The hybridization complex is separated from the hybridization system to obtain the targeted enriched nucleic acid fragment.

[0135] Specifically, using the complete mutation site corresponding to the heterogeneous drug resistance association characteristics and the 20bp conserved flanking sequences upstream and downstream as the targeting template, a 40bp high-affinity capture probe was designed using the base complementary pairing principle. The probe sequence completely covers the drug resistance characteristic mutation site without mismatches or degenerate bases, ensuring stable binding ability to the drug resistance mutation fragment. Using the strain-specific differential sequence corresponding to the lineage-specific characteristics as the targeting template, a 40bp distinguishing capture probe was designed. The probe accurately covers the lineage characteristic site and can specifically distinguish different strain lineage sequences. All the designed high-affinity capture probes and distinguishing capture probes were mixed at equimolar concentrations to generate an initial probe set containing both types of probes. For example, a high-affinity capture probe was designed for the rifampicin resistance mutation site of the rpoB gene, and a distinguishing capture probe was designed for the SNP site of the Mycobacterium tuberculosis lineage. The mixture constitutes the initial probe set.

[0136] The formula for calculating the thermodynamic free energy of nucleic acid molecules The binding free energy of each oligonucleotide probe in the initial probe set was evaluated, where... The binding free energy between the probe and the target sequence. This is the enthalpy change value of base pairing. This represents the entropy change value of base pairing. The standard reaction temperature for liquid-phase hybridization was fixed at 65℃. The complete flanking sequence of the corresponding characteristic site of the probe was retrieved as the binding template. The binding free energy of a single probe was calculated by substituting the corresponding thermodynamic parameters. The lower the free energy value, the stronger the probe binding stability.

[0137] For each heterogeneous drug resistance-associated characteristic site and lineage-specific characteristic site, the probe with the lowest binding free energy value among multiple candidate probes corresponding to that site is selected as the optimal capture probe for that characteristic site. All the original probes in the initial probe set corresponding to the site are replaced with the matching optimal capture probe. After completing the optimal replacement of all probes, all probe sequences and concentration parameters are uniformly standardized to generate a specific capture oligonucleotide probe set with the highest sequence matching degree and the best binding stability.

[0138] The prepared specific capture oligonucleotide probe set was added to a standard liquid hybridization buffer system and thoroughly mixed with the nucleic acid fragments obtained from the nucleic acid extraction and fragmentation of the sample to be identified. The hybridization reaction was carried out at a constant temperature of 65°C for 12 hours. This allowed the target fragments carrying heterogeneous drug resistance-associated features or lineage-specific features in the nucleic acid fragments to bind stably to the corresponding optimal capture probes through precise base complementarity pairing, forming a structurally stable probe-nucleic acid hybridization complex. Ordinary nucleic acid fragments without corresponding features did not generate specific binding and remained free in the hybridization system solution.

[0139] The reaction system after hybridization was separated and purified by streptavidin magnetic bead solid-phase adsorption separation. The magnetic beads were pre-modified with biotin and could specifically bind to the biotin-labeled groups at the probe ends. All hybridization complexes carrying the probe were specifically captured by magnetic adsorption. Unbound free nucleic acid fragments and impurity molecules in the system were discarded. After washing with buffer to remove non-specific adsorbed impurities, the nucleic acid fragments in the hybridization complex were eluted by high-temperature denaturation and dissociation. Finally, the purified nucleic acid fragments containing only the target feature sequence were obtained.

[0140] This step involves designing and optimizing a probe set based on identified gene characteristics, and using liquid-phase hybridization to achieve targeted enrichment of target nucleic acid fragments. This increases the effective proportion of heterogeneous drug resistance and lineage-specific fragments, dilutes irrelevant nucleic acid interference, provides high-quality templates for subsequent secondary validation sequencing, reduces sequencing detection difficulty, further reduces false negatives in the overall protocol, and ensures the reliability of the final identification results.

[0141] S6. Sequencing of the enriched target nucleic acid fragments, integration of mapping location information using the validation sequencing results, and outputting tuberculosis gene feature identification results.

[0142] In this embodiment of the invention, the step of sequencing the enriched target nucleic acid fragments, integrating the mapping position information using the verification sequencing results, and outputting tuberculosis gene feature identification results includes:

[0143] The targeted enriched nucleic acid fragments were sequenced to obtain a set of targeted sequencing reads. The set of targeted sequencing reads was then re-aligned to the reference genome index to obtain the mapping position of each targeted read on the genome.

[0144] The mapping position is compared with the genomic coordinates of the heterogeneous drug resistance association feature and the lineage-specific feature. The reads that overlap with any feature coordinate are selected as effective feature reads. The number of times each feature site is covered by the effective feature reads is counted as the verification confidence of the corresponding site.

[0145] The verification confidence of each feature site is compared with the previously obtained frequency distribution and signal intensity value. Feature sites that meet the consistency criteria are marked as confirmed features. All confirmed features and corresponding classification information are summarized to generate tuberculosis gene feature identification results.

[0146] Specifically, the targeted enriched nucleic acid fragments obtained in the aforementioned steps were sequenced at 150bp paired ends using the Illumina NovaSeq 6000 platform. Following the standard library preparation and sequencing workflow, base signal acquisition, signal correction, base identification, and quality control filtering were performed. The bcl2fastq software was used to... A value ≥30 is used as the filtering standard to remove invalid reads that do not meet the sequencing quality standards, ultimately generating a complete and clean set of targeted sequencing reads. The previously constructed graph reference genome index is used as the sole alignment template. The BWA-MEM algorithm is employed with default parameters, performing global multi-path alignments one by one between all reads in the targeted sequencing read set and the linear and parallel polymorphic reference paths in the graph reference genome index. The optimal mapping path with the highest weighted alignment score for each read is consistently selected as the final alignment result, accurately locating and recording the genomic mapping position and aligned base sequence corresponding to each targeted read.

[0147] We retrieved all precise genomic coordinate intervals corresponding to the previously identified heterogeneous drug resistance association features and lineage-specific features. The genomic mapping positions of each target read, including the start and end coordinates of the read coverage, were compared with all feature coordinate intervals. Target reads whose mapping intervals overlapped with any feature genomic coordinate at least by 1 bp were uniformly determined and screened as valid feature reads. The total number of times each feature genomic coordinate site was covered by valid feature reads was counted for each site, and this total number of coverages was directly defined as the validation confidence level for the corresponding feature site.

[0148] A pre-defined fixed consistency criterion is used for multi-dimensional data collaborative verification. This criterion is defined as follows: A mixed sample containing known low-frequency variants with frequencies of 1%, 2%, 5%, and 10% is selected. After processing through the first five steps of this workflow, targeted sequencing is performed. The positive predictive value and recall rate of the detection results are used as indicators. The optimal combination of three-dimensional parameters is determined through grid search, and is ultimately set as follows: a validation confidence level of not less than 20, a site frequency distribution value higher than the corresponding background error frequency threshold, and a signal intensity value greater than 0.05. The validation confidence level, the previously calculated frequency distribution, and the previously calculated signal intensity value corresponding to each feature site are simultaneously substituted into this criterion for three-dimensional collaborative comparison and verification. Feature sites that meet all three parameters are fixedly marked as confirmed features.

[0149] All verified features are collected, and the classification label corresponding to each verified feature is fully preserved: heterogeneous drug resistance association feature or lineage-specific feature, corresponding genomic coordinates, frequency distribution value, signal intensity value, and verification confidence parameter. All feature data are uniformly organized according to the ascending order of genomic coordinates, and all valid data are structured and summarized to finally generate tuberculosis gene feature identification results containing complete feature information and classification information.

[0150] This step involves targeted sequencing of enriched fragments, alignment with the genome index, and multi-dimensional collaborative verification using confidence levels combined with previous frequency distribution and signal intensity data to achieve secondary confirmation of characteristics. This overcomes the limitations of single-sequencing analysis, filters out false-positive variant sites, outputs reliable tuberculosis gene identification results, and improves the overall accuracy of drug resistance and pedigree characteristic identification.

[0151] like Figure 2 The diagram shown is a functional block diagram of a tuberculosis gene feature recognition system provided in an embodiment of the present invention.

[0152] The tuberculosis gene feature recognition system described in this invention can be installed in an electronic device. Depending on the functions implemented, the tuberculosis gene feature recognition system may include an entropy weight quality calibration module, a graph alignment mapping module, a haplotype typing module, a low-frequency feature extraction module, a probe targeting enrichment module, and a sequencing feature recognition module. The modules described in this invention can also be referred to as units, which are a series of computer program segments that can be executed by the processor of an electronic device and perform a fixed function, stored in the memory of the electronic device.

[0153] In this embodiment, the functions of each module / unit are as follows:

[0154] The entropy weight quality calibration module is used to obtain the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and to perform weight calibration on the base quality in the sequencing read set based on the entropy value of the context environment to obtain a quality score matrix.

[0155] The graph alignment and mapping module is used to construct a graph reference genome index containing parallel paths of polymorphism in the core genome and repetitive sequence regions based on the quality score matrix, perform multi-path alignment between the sequencing read set and the graph reference genome index, and generate the optimal mapping path and its mapping quality value.

[0156] The haplotype typing module is used to perform haplotype phase typing on the sequencing read set according to the mapped quality value and quality score matrix, separate different subclone reads into different phase blocks, and calculate the frequency distribution of characteristic sites in each phase block;

[0157] The low-frequency feature extraction module is used to identify low-frequency feature signals based on frequency distribution and in conjunction with a pre-constructed sequencing background error model, and to extract heterogeneous drug resistance association features and lineage-specific features from the low-frequency feature signals.

[0158] The probe-targeted enrichment module is used to synthesize a set of specific capture oligonucleotide probes based on heterogeneous drug resistance association characteristics and lineage-specific characteristics, which are used to hybridize with target fragments carrying the aforementioned characteristics to achieve targeted enrichment.

[0159] The sequencing feature recognition module is used to sequence the enriched target nucleic acid fragments, integrate the mapping position information using the verification sequencing results, and output the tuberculosis gene feature recognition results.

[0160] In the several embodiments provided by this invention, it should be understood that the disclosed methods and systems can be implemented in other ways. For example, the system embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and other division methods may be used in actual implementation.

[0161] The modules described as separate components may or may not be physically separate. The components shown as modules may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs.

[0162] Furthermore, the functional modules in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or in the form of hardware plus software functional modules.

[0163] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the present invention.

[0164] This application embodiment can acquire and process relevant data based on artificial intelligence technology. Artificial intelligence is the theory, method, technology, and application system that uses digital computers or machines controlled by digital computers to simulate, extend, and expand human intelligence, perceive the environment, acquire knowledge, and use that knowledge to obtain optimal results.

[0165] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for identifying tuberculosis gene features, characterized in that, The method includes: S1. Obtain the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and perform weight calibration on the base quality in the sequencing read set based on the entropy value of the context environment to obtain the quality score matrix. S2. Based on the quality score matrix, construct a graph reference genome index containing parallel paths of polymorphism in the core genome and repetitive sequence regions. Perform multi-path alignment between the sequencing read set and the graph reference genome index to generate the optimal mapping path and its mapping quality value. S3. Based on the mapping quality value and quality score matrix, perform haplotype phase typing on the sequencing read set, separate different subclone reads into different phase blocks, and calculate the frequency distribution of characteristic sites in each phase block; S4. Based on the frequency distribution and combined with the pre-constructed sequencing background error model, identify low-frequency feature signals, and extract heterogeneous drug resistance association features and lineage-specific features from the low-frequency feature signals; S5. Based on the heterogeneity of drug resistance association characteristics and lineage specificity characteristics, synthesize a set of specific capture oligonucleotide probes for hybridization with target fragments carrying the aforementioned characteristics to achieve targeted enrichment; S6. Sequencing of the enriched target nucleic acid fragments, integrating mapping location information using the validation sequencing results, and outputting tuberculosis gene feature identification results.

2. The tuberculosis gene feature identification method as described in claim 1, characterized in that, The process involves obtaining the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and then weighting and calibrating the base quality in the sequencing read set based on the entropy value of the context environment to obtain a quality score matrix, including: Based on each base site of each read in the sequencing read set, the base type combination of the upstream and downstream neighboring bases of the base site is extracted as contextual information, and the initial base quality value of the base site is obtained. Based on the context information, the base identification error frequency of all base sites with the same context information in the sequencing read set is statistically analyzed, and the empirical error rate of the context is calculated based on the base identification error frequency. The entropy weighting factor of the context is determined based on the deviation between the empirical error rate and the expected error rate characterized by the initial base quality value. The initial base quality value of the base site is weighted and corrected according to the entropy weight factor corresponding to the context of the base site to obtain the calibrated base quality value. The calibrated base quality values ​​are then summarized to construct a quality score matrix.

3. The tuberculosis gene feature identification method as described in claim 1, characterized in that, Based on the quality score matrix, a graph reference genome index is constructed, containing parallel paths for polymorphism in the core genome and repetitive sequence regions. The sequencing read sets are then compared with the graph reference genome index via multiple paths to generate the optimal mapping path and its mapping quality value, including: Based on the distribution of base quality values ​​in the quality score matrix across the genome, the sequencing read set is divided into high-confidence alignment regions and low-confidence alignment regions. For the high-confidence alignment regions, a single linear reference path is constructed; for the low-confidence alignment regions, known polymorphic sequence variants are obtained and multiple parallel reference paths are constructed. All paths are then integrated into a graph reference genome index. Each read in the sequencing read set is aligned with each reference path in the graph reference genome index to obtain an initial alignment score for each read on the path. Based on the quality values ​​of the read length at each alignment position in the quality score matrix, the initial alignment score is weighted and corrected to obtain the weighted alignment total score for each path; The path with the highest weighted total score is selected as the optimal mapping path, and the mapping quality value of the optimal mapping path is generated based on the difference in weighted total score between the optimal path and the second-best path.

4. The tuberculosis gene feature identification method as described in claim 3, characterized in that, For the high-confidence alignment regions, a single linear reference path is constructed; for the low-confidence alignment regions, known polymorphic sequence variants are obtained and multiple parallel reference paths are constructed. All paths are integrated into a graph reference genome index, including: Based on the quality value of each read locus in the quality score matrix, the average quality value of the reads covered by each genomic region is calculated. Regions with an average quality value higher than a preset reliability threshold are classified as high-confidence alignment regions, and regions with an average quality value lower than the reliability threshold are classified as low-confidence alignment regions. For the high-confidence alignment regions, a single linear reference path is constructed using the reference genome sequence of the corresponding region as a template; For the low-confidence alignment region, all known polymorphic sequence variants in the corresponding region are obtained, and parallel reference paths are constructed for each polymorphic sequence variant. Each parallel path establishes a connection node with the adjacent region of the single linear reference path at the end of the path, generating a preliminary graph index containing a branch structure. The single linear reference path and the parallel reference path are arranged in genomic coordinate order. The paths are connected end to end through the connection nodes between each path to generate a connected graph structure that covers the core genome and repetitive sequence regions and includes polymorphic branch paths. The connected graph structure is used as the graph reference genome index.

5. The tuberculosis gene feature identification method as described in claim 1, characterized in that, The step involves performing haplotype phase typing on the sequencing read set based on the mapped quality value and quality score matrix, separating different subclonal reads into different phase blocks, and calculating the frequency distribution of characteristic sites in each phase block, including: Based on the mapped quality value, the effective read lengths for genotyping are selected, and the calibrated quality value is used as the base confidence weight to generate a weighted base sequence for each effective read length. Based on the ratio of the sum of the weighted product of two reads with identical bases at the common site to the total sum of the weighted product, a weighted consistency association graph is constructed, and the connected regions in the graph are defined as the initial phase blocks. The base divergence of candidate sites within the initial phase block is detected. Based on the cumulative weight of the base divergence, low-weight erroneous bases are discarded, and the consensus sequence of the corresponding block is reconstructed based on the remaining read length. Based on the consensus sequence, the read length assignment is corrected, and the read lengths that still cannot be assigned to the existing block after correction are formed into new phase blocks separately, thus obtaining the final phase block; The cumulative confidence weight of each base type at each feature site within the final phase block is calculated, and the proportion of the cumulative weight of each base type to the total cumulative weight of the corresponding site is used as the frequency distribution of the corresponding feature site.

6. The tuberculosis gene feature identification method as described in claim 1, characterized in that, Based on the frequency distribution and combined with a pre-constructed sequencing background error model, low-frequency feature signals are identified, and heterogeneous drug resistance association features and lineage-specific features are extracted from the low-frequency feature signals, including: The frequency distribution of feature sites in the phase block is compared one by one with the background error frequency threshold of the corresponding site context in the pre-constructed sequencing background error model. Base types whose frequency values ​​exceed the preset threshold are selected as candidate low-frequency signal bases, and the excess amount of the candidate low-frequency signal bases is recorded as the signal intensity value. Based on the number of read lengths supported by candidate low-frequency signal bases in the corresponding phase blocks and the signal intensity value, the micro-enrichment of candidate low-frequency signal bases in the subclonal population is calculated. Candidate low-frequency signal bases with micro-enrichment levels higher than a preset standard are identified as true low-frequency feature signals, thus obtaining a set of true low-frequency feature signals. Based on the genomic coordinates of each point in the real low-frequency feature signal set, signals located in known drug resistance gene regions are classified as heterogeneous drug resistance associated features, and signals located in phylogenetic marker gene regions are classified as lineage-specific features.

7. The tuberculosis gene feature identification method as described in claim 6, characterized in that, The pre-constructed sequencing background error model includes: The conserved region sequence of the Mycobacterium tuberculosis genome was obtained as the background training base, and reads that were consistent with the background training base were extracted from historical sequencing data as the training read set. The occurrence frequency of each type of base that is inconsistent with the reference base of the training read set is counted, and the occurrence frequency is divided by the total number of times the corresponding site is covered by the training read to obtain the empirical error rate of each base substitution type at the corresponding site. Using the upstream and downstream adjacent base sequences of each site as context markers, sites with the same context markers are grouped into the same context category; The empirical error rates of each base substitution type at all sites within the same environmental category are summarized, and the summarized results are used as the background error frequency threshold for the corresponding environmental category. The correspondence between environmental categories and background error frequency thresholds is established, and a sequencing background error model is constructed.

8. The tuberculosis gene feature identification method as described in claim 1, characterized in that, The method involves synthesizing a specific capture oligonucleotide probe set based on heterogeneous drug resistance association characteristics and lineage-specific characteristics, for hybridization with target fragments carrying the aforementioned characteristics to achieve targeted enrichment, including: Based on the sequence composition of the heterogeneous drug resistance association features and the lineage-specific features, high-affinity capture probes are designed for the drug resistance association features and discriminative capture probes are designed for the lineage-specific features, and the two types of probes are mixed to form an initial probe set; The binding free energy of each probe in the initial probe set with the corresponding flanking sequence of the characteristic site is evaluated. The probe with the lowest binding free energy is selected as the optimal capture probe for the characteristic site and replaces the corresponding probe in the initial probe set to generate a specific capture oligonucleotide probe set. The specific capture oligonucleotide probe set is liquid-phase hybridized with the nucleic acid fragment of the sample to be identified, so that the target fragment carrying heterogeneous drug resistance association characteristics or lineage specific characteristics forms a hybridization complex with the corresponding probe. The hybridization complex is separated from the hybridization system to obtain the targeted enriched nucleic acid fragment.

9. The tuberculosis gene feature identification method as described in claim 1, characterized in that, The process involves sequencing the enriched target nucleic acid fragments, integrating the mapping location information using the validation sequencing results, and outputting tuberculosis gene feature identification results, including: The targeted enriched nucleic acid fragments were sequenced to obtain a set of targeted sequencing reads. The set of targeted sequencing reads was then re-aligned to the reference genome index to obtain the mapping position of each targeted read on the genome. The mapping position is compared with the genomic coordinates of the heterogeneous drug resistance association feature and the lineage-specific feature. The reads that overlap with any feature coordinate are selected as effective feature reads. The number of times each feature site is covered by the effective feature reads is counted as the verification confidence of the corresponding site. The verification confidence of each feature site is compared with the previously obtained frequency distribution and signal intensity value. Feature sites that meet the consistency criteria are marked as confirmed features. All confirmed features and corresponding classification information are summarized to generate tuberculosis gene feature identification results.

10. A tuberculosis gene feature identification system, characterized in that, The system for implementing the tuberculosis gene feature identification method according to claim 1, the system comprising: The entropy weight quality calibration module is used to obtain the sequencing read set of Mycobacterium tuberculosis in the sample to be identified, and to perform weight calibration on the base quality in the sequencing read set based on the entropy value of the context environment to obtain the quality score matrix. The graph alignment and mapping module is used to construct a graph reference genome index containing parallel paths of polymorphism in the core genome and repetitive sequence regions based on the quality score matrix, and to perform multi-path alignment between the sequencing read set and the graph reference genome index to generate the optimal mapping path and its mapping quality value. The haplotype typing module is used to perform haplotype phase typing on the sequencing read set according to the mapped quality value and quality score matrix, separate different subclone reads into different phase blocks, and calculate the frequency distribution of characteristic sites in each phase block; The low-frequency feature extraction module is used to identify low-frequency feature signals based on frequency distribution and in conjunction with a pre-constructed sequencing background error model, and to extract heterogeneous drug resistance association features and lineage-specific features from the low-frequency feature signals. The probe-targeted enrichment module is used to synthesize a set of specific capture oligonucleotide probes based on heterogeneous drug resistance association characteristics and lineage-specific characteristics, which are used to hybridize with target fragments carrying the aforementioned characteristics to achieve targeted enrichment. The sequencing feature recognition module is used to sequence the enriched target nucleic acid fragments, integrate the mapping position information using the verification sequencing results, and output the tuberculosis gene feature recognition results.

Images