A method and system for detecting a CYP21A2 gene mutation and a storage medium

Abstract

The application provides a CYP21A2 gene mutation detection method, system and storage medium, comprising the following steps: identifying difference sites and homologous intervals by reference genome self alignment, and constructing a paralog characteristic site fingerprint library; constructing a Masker reference genome according to the homologous intervals on the reference genome, and obtaining a candidate variation list; using the physical link information of a double-end sequencing read pair and the characteristic site fingerprint library to determine the true or false gene source of each candidate variation; quantitatively calculating the true gene copy number based on re-alignment sequencing depth, background reference library and characteristic site allele frequency; and determining whether a recombination event occurs between the true gene and the false gene based on the physical link information of the double-end sequencing read pair across multiple characteristic sites. The application effectively distinguishes different copy states such as gene deletion, single copy, normal double copy and duplication, and assists in clinically determining gene deletion or duplication events.

Description

Technical Field

[0001] This invention relates to the field of molecular diagnostics and bioinformatics analysis of genetic diseases, specifically to a method, system, and storage medium for detecting CYP21A2 gene mutations. Background Technology

[0002] The CYP21A2 gene, located in the 6p21.3 region of human chromosome 21, encodes 21-hydroxylase, a key gene for maintaining endocrine homeostasis. When the CYP21A2 gene malfunctions, leading to a deficiency of 21-hydroxylase, it is inherited in an autosomal recessive manner. Classical cases can result in adrenal crisis, which can be life-threatening. Hyperandrogenemia can cause masculinization in women, leading to accelerated bone age progression, short stature, abnormal puberty development, and impaired fertility.

[0003] The challenge of molecular diagnosis of this gene locus stems from its unique genomic structure. In the 6p21.3 region, the functional gene CYP21A2 is closely adjacent to its highly homologous pseudogene CYP21A1P, exhibiting a tandem repeat distribution, with exon sequence homology exceeding 98%. Due to the high sequence similarity between the true and pseudogenes, unequal crossing over easily occurs during meiosis, leading to the transfer of harmful mutant sequences from the pseudogene to the functional gene CYP21A2, triggering gene conversion. Furthermore, this region is also prone to gene recombination and deletion; these multiple variation types collectively constitute the core technical challenge of CYP21A2 gene detection.

[0004] Whole-exome sequencing (WES) has been widely adopted and applied in the clinical diagnosis of genetic diseases due to its ability to perform unbiased capture sequencing of almost all protein-coding regions in the genome in a single operation. However, under current next-generation sequencing (NGS) conditions, limited by the PE150 read length, conventional alignment algorithms commonly exhibit multiple alignments when handling highly homologous regions of CYP21A2 and CYP21A1P. A large number of valid reads are incorrectly assigned to pseudogene regions or discarded after being marked as low-quality alignments, resulting in a sharp decrease in the effective coverage depth of the target region, difficulty in determining the mutation source, and genotyping errors. Currently, three main solutions are used clinically to address these problems, but all have significant limitations. One approach is WES analysis based on the standard GATK workflow, using a standard reference genome combined with BWA-MEM for read alignment, followed by SNV and Indel detection using general variant detection tools such as GATKHaplotypeCaller or VarScan. This approach suffers from systematic missed detections in highly homologous regions and treats each read as an independent statistical event, lacking the ability to mine information about the physical links between adjacent variant sites. The second approach is long-fragment PCR combined with Sanger sequencing or NGS, which uses specific primers to amplify across homologous regions. This method not only increases experimental steps and requires extremely high-quality DNA templates, but also significantly increases the complexity of constructing NGS libraries after fragmentation, making it difficult to meet the efficiency requirements of large-scale clinical screening. Thirdly, ultra-long reads can be obtained using third-generation sequencing technologies such as PacBioSMRT or Oxford Nanopore, which can theoretically directly span the full-length CYP21A2 region. However, Oxford Nanopore still has a high random error rate in the homopolymer region. Although PacBio HiFi has high accuracy, it is expensive, has limited throughput, and is not widely used. Currently, most clinical laboratories still mainly use second-generation sequencing platforms.

[0005] In summary, existing technologies present inherent contradictions in terms of sensitivity for detecting mutations in highly homologous regions, ability to trace true and false gene variants, and accessibility for large-scale clinical application. Therefore, developing a novel bioinformatics analysis method that can fully utilize existing WES data and overcome the bottleneck of homologous sequence alignment is an urgent clinical need and a key direction for technological breakthroughs in achieving non-invasive, accurate, and one-stop molecular diagnosis of the CYP21A2 gene. Summary of the Invention

[0006] To overcome the shortcomings of existing technologies, this invention proposes a method for detecting CYP21A2 gene mutations, comprising the following steps: S1. By comparing with the reference genome, identify the differential sites and homologous regions between the CYP21A2 true gene and the CYP21A1P pseudogene, and construct a fingerprint database of paralogous characteristic sites. S2. Based on the homology interval, the reference genome is masked to construct the Masker reference genome. The whole exome sequencing reads are re-aligned to the Masker reference genome to obtain a candidate variant list. S3. Using the physical link information of paired-end sequencing read pairs, combined with the feature site fingerprint database, determine the true and false gene origins of each variant site in the candidate variant list, and obtain the variant support information of the CYP21A2 true gene and the CYP21A1P false gene respectively. S4. Based on the re-alignment sequencing depth and background reference library, and combined with the allele frequency of the characteristic sites in the characteristic site fingerprint database, the CYP21A2 true gene copy number is quantitatively estimated. S5. Based on the physical link information of paired-end sequencing reads spanning multiple characteristic sites, assemble linked haplotypes and determine whether recombination events have occurred between true and false genes.

[0007] Further, the reference genome self-alignment in step S1 is implemented as follows: the target locus region sequence is divided into overlapping sub-fragments using a sliding window algorithm, the sub-fragments are back-aligned to the reference genome, fragments with high-scoring alignments at multiple positions are screened, multiple sequence alignments are performed on the identified homologous regions, and stable differential sites between CYP21A2 and CYP21A1P are extracted to form the feature site fingerprint database. The sliding window algorithm has a window length of 500 to 5000 bases and an overlap step of 10% to 90% of the window length.

[0008] Furthermore, the masking process described in step S2 is a hard masking process that replaces the sequence of the CYP21A1P pseudogene homologous region in the reference genome with the character 'N', so that the read segment that was originally ambiguous due to multiple alignments is uniquely aligned to the CYP21A2 true gene region. The masking process can also be replaced by a soft masking process, which involves converting the reference sequence of the CYP21A1P pseudogene homologous region into lowercase letters and using the difference in processing uppercase and lowercase sequences by the alignment software to reduce the probability of misaligning the read segment to the pseudogene region.

[0009] Furthermore, the determination of true and false gene origins in step S3 includes the following steps: Step S31: Establishment of molecular source tracing anchors: Using the feature site fingerprint database constructed in step S1 as a reference, extract the specific alleles corresponding to the CYP21A2 true gene and the CYP21A1P pseudo gene at each feature site to form a set of molecular source tracing anchors. Step S32: Physical link read screening. For each candidate variant site, among the re-aligned reads obtained in step S2, reads that can be used for tracing are screened in the following two ways: Step S321: Direct coverage, the comparison interval of the read segment itself simultaneously covers at least one anchor point in the candidate mutation site and the molecular tracing anchor point set; Step S322: Physical link extension, the alignment interval of the read itself covers the candidate variant site, and the alignment interval of its paired reads covers at least one anchor point in the molecular origin anchor point set; based on the physical co-origin relationship of the paired reads established by paired sequencing, the candidate variant and the anchor point are included in the same alignment coordinate range and associated. The method of step S321 is preferred; when the number of available read segments in step S321 is lower than the preset coverage depth threshold, the read segments screened by the method of step S322 are added. Regarding the switching threshold between the direct coverage method and the physical link extension method, the present invention sets the default value of the preset coverage depth threshold to 10×. That is, when the number of available traceability reads obtained by screening through the direct coverage method step S321 at a certain candidate mutation site is less than 10, the physical link extension method step S322 is initiated to supplement and include paired reads to increase the number of available traceability reads.

[0010] The threshold is set based on the following: In conventional WES samples with a sequencing depth of at least 50×, the effective coverage depth of the CYP21A2 target region after Masker realignment is typically between 30× and 80×. The number of source reads obtained by direct coverage can meet the minimum requirement of 10 reads at most sites. Only in exon margin regions or local regions with low capture efficiency may the number of directly covered reads be lower than this threshold. In this case, enabling physical linking extension can effectively supplement source information without introducing additional attribution errors. In practical applications, the preset coverage depth threshold can be adjusted within the range of 5× to 20× according to the overall sequencing depth of the sample and the uniformity of target region capture. For low-depth samples with an average sequencing depth of less than 30×, it is recommended to lower the threshold to 5× and prioritize physical linking extension to ensure a sufficient number of source reads.

[0011] Step S33: Source attribution determination. Read the allele bases of the segment screened in step S32 at the molecular source anchor point and compare them one by one with the specific alleles of the CYP21A2 true gene and CYP21A1P pseudo gene in the characteristic site fingerprint database. When a segment only covers a single anchor point, the attribution is determined directly based on the matching result of that anchor point. When a segment covers two or more anchor points and all anchor points point to the same attribution conclusion, the attribution is determined based on the consistent conclusion. When a segment covers two or more anchor points and the anchor points point to inconsistent attribution conclusions, the segment with more than half of the anchor points supporting the same attribution conclusion shall be used; if neither side exceeds half, the segment shall be marked as having uncertain attribution and excluded from subsequent statistics. Based on the above matching results, the reads and their candidate variants are assigned to either the CYP21A2 true gene or the CYP21A1P pseudo gene.

[0012] Furthermore, in step S3, after the molecular source tracing anchor set is established in step S31 and before the physical link read segment screening is initiated in step S32, the candidate variant list is pre-classified based on the characteristic site fingerprint database, and differentiated molecular source tracing anchor strategies are selected for different categories: Characteristic site mutation type: The genomic coordinates of the candidate variant match a differential site in the characteristic site fingerprint database, and the allele type is consistent with the known allele status of the differential site; when tracing the origin of such mutations, other characteristic sites besides the differential site are used as molecular tracing anchors, excluding the site itself from participating in the determination of the origin of the mutation. Regarding the criteria for determining the mutation type at a characteristic site, the operation definition of "the genomic coordinates of the candidate variant match a differential site in the characteristic site fingerprint database" in this invention is as follows: Under the same reference genome version, if the starting coordinates of the candidate variant are completely consistent with the recorded coordinates of a differential site in the characteristic site fingerprint database, and the allele type of the candidate variant (including reference bases and substitution bases) is consistent with the known allele status of the differential site recorded in the characteristic site fingerprint database, then it is determined that the coordinates match.

[0013] When the reference genome version used for the sample being tested differs from the version used to construct the fingerprint database (e.g., the analysis workflow uses GRCh37 while the fingerprint database is constructed based on GRCh38), the fingerprint database coordinates should first be converted to the corresponding version using a coordinate transformation tool (such as LiftOver) before consistency determination. After coordinate transformation, the initial coordinates should still be completely consistent as the judgment criterion, without setting a coordinate offset tolerance, to avoid misjudgments caused by different indel representation standards (left alignment and right alignment). For indel-type feature sites, allele consistency determination should uniformly adopt the left alignment standard to standardize candidate variants and fingerprint database records before comparison.

[0014] Novel site mutations: The genomic coordinates and allele types of candidate variants are not found in the characteristic site fingerprint database; when tracing the origin of such mutations, any one or more characteristic sites in the characteristic site fingerprint database are used as molecular tracing anchors. The source attribution rules corresponding to their respective categories are applied to the two types of mutations.

[0015] Furthermore, the method for estimating the CYP21A2 true gene copy number in step S4 is as follows: The total copy number of CYP21A2 and CYP21A1P was obtained by performing in-sample data volume normalization, GC bias correction, and secondary correction based on the background reference library on the weight ratio depth. ; Extract the pseudogene allele frequencies at all available feature sites in the feature site fingerprint database, and use the median of the frequencies. As a robust estimator; The CYP21A2 true gene copy number can be calculated using the following formula. : .

[0016] Further, the assembly method of the linked haplotype in step S5 is as follows: screen paired reads spanning at least two feature sites, record the allele status at each site in sequence according to the feature site, and gradually splice and extend to form a haplotype sequence covering the full length of CYP21A2; compare the obtained haplotype sequence with the CYP21A2 reference haplotype position by position, and if the combination of feature sites in a continuous interval is consistent with the CYP21A1P pseudogene haplotype, it is determined that a recombination event has occurred in the interval, and the genomic intervals corresponding to adjacent feature sites are marked as the estimated range of recombination breakpoints.

[0017] Furthermore, the feature site fingerprint database can also be constructed using one of the following alternative methods: (a) K-mer analysis: Specific K-mers appearing only in the CYP21A2 eugene and specific K-mers appearing only in the CYP21A1P pseudogene were counted, and the specific K-mers were used as differential sites, with K values ​​ranging from 31 to 63. (b) Static database method: The differential sites of CYP21A2 and CYP21A1P already marked in the public variant database are used as the initial feature site fingerprint database, and local calibration is performed in combination with sample data.

[0018] Furthermore, the default value of the allele status attribution threshold is 50%, that is, when the number of supporting reads of a certain allele at a characteristic site accounts for more than 50% of the total number of reads covered by that site, the read is determined to originate from the corresponding gene; the attribution threshold can be adjusted within the range of 40% to 70% according to the sample sequencing depth and quality control indicators.

[0019] Furthermore, the rule for determining the attribution of the characteristic site mutation is as follows: (a) If the read does not contain the candidate mutation and the proportion of the true gene allele at another characteristic site exceeds the attribution threshold, the read is determined to originate from the CYP21A2 true gene, and the characteristic site is wild-type in the true gene. (b) If a read contains the candidate mutation and the proportion of the genotype allele at another characteristic site exceeds the attribution threshold, the read is determined to originate from the CYP21A2 genotype, and the characteristic site is a mutant in the genotype. (c) If a read contains the candidate mutation and the proportion of pseudogene alleles at another characteristic site exceeds the attribution threshold, the read is determined to originate from the CYP21A1P pseudogene, and the characteristic site is wild-type in the pseudogene. (d) If the read does not contain the candidate mutation and the proportion of pseudogene alleles at another characteristic site exceeds the attribution threshold, the read is determined to originate from the CYP21A1P pseudogene, and the characteristic site is a mutant in the pseudogene.

[0020] Furthermore, the rule for determining the attribution of the newly occurring mutation site is as follows: (a) If the read does not contain the candidate mutation and the proportion of the true gene allele at the characteristic site exceeds the attribution threshold, the read is determined to originate from the CYP21A2 true gene, and the new site is wild-type in the true gene. (b) If a read contains the candidate mutation and the proportion of the true gene allele at the characteristic site exceeds the attribution threshold, the read is determined to originate from the CYP21A2 true gene, and the new site is a mutant in the true gene. (c) If the read does not contain the candidate mutation and the proportion of pseudogene alleles at the characteristic site exceeds the attribution threshold, the read is determined to originate from the CYP21A1P pseudogene, and the new site is wild-type in the pseudogene. (d) If a read contains the candidate mutation and the proportion of pseudogene alleles at the characteristic site exceeds the attribution threshold, the read is determined to originate from the CYP21A1P pseudogene, and the new site is a mutant in the pseudogene.

[0021] Furthermore, the construction of the background reference library satisfies the following conditions: (a) The number of samples included shall not be less than 30, and the CYP21A2 copy number of the samples shall be confirmed to be normal by an independent verification method; (b) The same DNA extraction method, library construction kit and sequencing platform were used for the included samples and the samples to be tested; (c) The average sequencing depth of the target region in the included samples is not less than 50×, and the exome capture uniformity is not less than 85%; (d) After the background library is constructed, the depth variation coefficient of each target interval in all background samples does not exceed 0.15; (e) When changing sequencing batches, the background library is updated on a rolling basis to maintain consistency with the construction conditions of the sample library to be tested.

[0022] A CYP21A2 gene mutation detection system, applicable to the aforementioned CYP21A2 gene mutation detection method, comprising: The feature site fingerprint database construction module (first module) is used to identify the differential sites and homologous regions between the CYP21A2 egenic gene and the CYP21A1P pseudogene by self-alignment with the reference genome, and to construct a paralogous feature site fingerprint database. The alignment reconstruction and variant detection module (second module) is used to mask the reference genome based on the homologous regions, construct the Masker reference genome, re-align the whole exome sequencing reads to the Masker reference genome, and obtain a candidate variant list. The physical link attribution determination module (third module) is used to determine the true and false gene origins of each candidate variant by utilizing the physical link information of paired-end sequencing read pairs and the feature site fingerprint database. The copy number assessment module (module 4) is used to quantitatively estimate the CYP21A2 eukaryotic copy number based on the re-aligned sequencing depth, background reference library, and allele frequency of characteristic sites. The haplotype assembly and recombination detection module (the fifth module) is used to assemble linked haplotypes based on the physical link information of paired-end sequencing reads spanning multiple characteristic sites, and to determine whether recombination events have occurred between true and false genes.

[0023] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the method described above.

[0024] Compared with the prior art, the beneficial effects of the present invention are: 1. This invention provides a method, system, and storage medium for detecting CYP21A2 gene mutations. By masking the homologous region of the CYP21A1P pseudogene in the reference genome, it fundamentally eliminates the problem of read discarding caused by multiple alignments in conventional alignment algorithms. It forcibly recovers sequencing signals that were originally marked as low alignment quality and discarded to the CYP21A2 eukaryotic region, effectively restoring the coverage depth of key variant regions and significantly reducing the risk of missed detection in clinical samples.

[0025] 2. This invention provides a method, system, and storage medium for detecting CYP21A2 gene mutations. By fully exploring the physical link information carried by the PE150 paired-end sequencing reads, the candidate variant sites are cis-bound to anchor points in the characteristic site fingerprint database, achieving accurate determination of the true and false gene origins. This fundamentally solves the problem of false positives caused by pseudogene contamination in CYP21A2 diagnosis and provides haplotype-supported tracing results for clinical variant interpretation.

[0026] 3. This invention provides a method, system, and storage medium for detecting CYP21A2 gene mutations. By integrating sequencing depth normalization, GC bias correction, background reference library correction, and allele frequency analysis of characteristic sites, a complete algorithm system for quantitatively estimating the true copy number of the CYP21A2 gene is established. This system can effectively distinguish different copy states such as gene deletion, single copy, normal double copy, and duplication, and assist in the clinical determination of gene deletion or duplication events.

[0027] 4. This invention provides a method, system, and storage medium for detecting CYP21A2 gene mutations. By assembling linked haplotypes spanning multiple characteristic sites, it can directly detect recombination events between true and false genes and locate recombination breakpoint intervals, filling the gap in existing WES analysis procedures for CYP21A2 recombination detection.

[0028] 5. This invention provides a method, system, and storage medium for detecting CYP21A2 gene mutations. It is based entirely on conventional WES sequencing data, requiring no additional long-fragment PCR experiments or third-generation sequencing platforms. It can be directly connected to existing automated sequencing production lines, supports large-scale parallel sample processing, and has significant cost advantages and clinical applicability. Attached Figure Description

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

[0030] Figure 1 This is a schematic diagram of the process of this invention; Figure 2 This is a schematic diagram illustrating the four possible attributions of novel site mutations; Figure 3 A schematic diagram illustrating four possible attributions of the characteristic site mutation; Figure 4 This is a schematic diagram of the detection of chain-type assembly and recombination events; Figure 5This is a VAF correction effect verification chart (three-layer comparison frame); Figure 6 This is a four-layer comparison diagram of mutation frequencies in the feature site fingerprint database of Example 1; Figure 7 This is a four-layer comparison diagram of the mutation frequency of the feature site fingerprint database in Example 2; Figure 8 This is a copy number distribution diagram of 78 samples from the background library; Figures 9-10 This is a copy number distribution diagram of 9 samples after correction. Detailed Implementation

[0031] The technical solution of the present invention will be more clearly and completely explained below with reference to the accompanying drawings and through the description of preferred embodiments of the present invention.

[0032] To ensure that the technical solutions and effects of this invention can be fully understood and implemented by those skilled in the art, the following detailed description of each module of this invention is provided in conjunction with specific embodiments. This invention is described based on the human reference genome GRCh38; however, those skilled in the art should understand that the method described in this invention can also be implemented under other versions of the reference genome (including but not limited to GRCh37 and T2T-CHM13) after coordinate transformation. The choice of reference genome version does not constitute a limitation on the scope of protection of this invention.

[0033] The core technical challenge addressed in this invention lies in the over 98% exon sequence homology between the CYP21A2 true gene and the CYP21A1P pseudogene. This extremely high homology causes numerous multiple alignments when conventional NGS alignment software processes the PE150 reads, resulting in the loss of effective sequencing signals, inability to distinguish mutation origins, and distortion of copy number estimation. This invention systematically solves these problems through the following interconnected technical modules.

[0034] like Figure 1 As shown, the first module is the feature site fingerprint database construction module. It receives a reference genome and WES data as input, and identifies differentially expressed sites and homologous regions between the CYP21A2 true gene and the CYP21A1P pseudogene through self-alignment with the reference genome, constructing a paralogous homologous feature site fingerprint database to provide a molecular anchor basis for all subsequent modules. The second module is the alignment reconstruction and variant detection module. Based on the homologous regions marked in the first module, it constructs a Masker reference genome, realigns the WES reads to this reference genome, and then calls a variant detection tool to obtain a list of candidate variants, fundamentally solving the problem of sequencing signal loss caused by multiple alignments. The third module is the physical link attribution determination module. It comprehensively utilizes the feature site fingerprint database constructed in the first module and the physical link information of the paired-end sequencing reads to determine the true and false gene origins of each candidate variant obtained in the second module. like Figure 1 As shown by the dashed loop on the right, the feature site fingerprint database of the second module is directly reused as a molecular tracing anchor, demonstrating its core and fundamental role throughout the entire process. After the third module completes the attribution determination, the workflow is divided into two parallel branches: the fourth module is the copy number assessment module, which quantitatively estimates the CYP21A2 true gene copy number based on re-alignment sequencing depth, background reference library, and allele frequencies at feature sites, to assist in determining gene deletions or duplications; the fifth module is the haplotype recombination detection module, which assembles linked haplotypes based on the physical linking information of paired reads spanning multiple feature sites to determine whether recombination events have occurred between true and false genes. The analysis results of the fourth and fifth modules are combined to output a complete CYP21A2 gene detection report, including variant sites and genotypes, quantitative copy number estimation, haplotype sequences, and recombination breakpoint intervals.

[0035] The first module involves the construction of a fingerprint database of paralogous characteristic sites. The characteristic site fingerprint database is the data foundation for all subsequent analyses in this invention. Its core objective is to systematically extract stable cis-differentiated sites between CYP21A2 and CYP21A1P from the reference genome as molecular anchors to distinguish between true and false gene origins.

[0036] In practice, a sliding window algorithm is used to segment the reference sequence of chromosome 6p21.3 covering the CYP21A2 and CYP21A1P loci. The window length is set to 500 to 5000 bases, and the overlap step size is set to 10% to 90% of the window length, generating a set of overlapping sub-fragments covering the target region. The resulting set of sub-fragments is then compared back to the original reference genome using the BWA-MEM algorithm. Fragments with alignment quality values ​​below a preset threshold (default threshold is 20) or with multiple high-scoring secondary alignment positions are selected. Special attention is paid to recording fragments whose initial alignment position is in the CYP21A2 region but also have high-scoring alignments in the CYP21A1P region, as well as fragments in the reverse case, thereby completing a comprehensive survey of homologous regions.

[0037] Multiple sequence alignment was performed on the homologous fragments obtained from the above screening to extract stable single-base differential sites (SNP-type characteristic sites) and short fragment insertion / deletion sites (Indel-type characteristic sites) between CYP21A2 and CYP21A1P. After filtering, each characteristic site was required to have a clear true and pseudo-genotype at the reference genome level and to exhibit stable differential expression in the population context. Taking GRCh38 as an example, the number of characteristic sites that can be identified by this module is usually between 30 and 60, covering the key regions of exons 1 to 10 of CYP21A2. Table 2 lists some of the characteristic loci identified by this module under the GRCh38 reference genome in this invention. Example characteristic loci include: at chr6:32038297, the true genotype is C and the pseudogenotype is T, with the corresponding pseudogenotype location at chr6:32005563; at chr6:32038310, the true genotype is GGGT and the pseudogenotype is AGGC, with the corresponding pseudogenotype location at chr6:32005576, etc. These characteristic loci and their corresponding homologous regions together constitute a characteristic locus fingerprint database for subsequent module use. As shown in Table 1: Table 1:

[0038] Table 1 shows the direct output of the first module, which is a specific display of the feature site fingerprint database. It serves as an example of the technical implementation results of this module and also as a reference for anchor point data. In an alternative embodiment of the present invention, the characteristic site fingerprint database can also be constructed using K-mer analysis, that is, specific K-mer sequences (K value preferably 31 to 63) that appear only in the CYP21A2 eukaryotic region and not in the CYP21A1P pseudogene region are statistically analyzed, and these specific K-mers are used as differential site identifiers. This method can provide higher operating efficiency when computational resources are limited. In another alternative embodiment, the CYP21A2 and CYP21A1P differential sites that have been clearly marked in public variant databases (such as dbSNP, ClinVar) can be used as the initial characteristic site fingerprint database, and then locally calibrated in combination with actual sample data. This method is suitable for clinical laboratory scenarios with stable historical data accumulation.

[0039] The second module involves alignment reconstruction and candidate variant detection based on the Masker reference genome. The core objective of this module is to recover sequencing reads that were originally discarded due to multiple alignments into the CYP21A2 eukaryotic region by shielding pseudogene interference regions, thereby establishing a complete data foundation for variant detection.

[0040] Quantitative performance comparison with conventional GATK analysis workflow To quantitatively evaluate the improvement in technical performance of the present invention compared to conventional analysis procedures, the present invention simultaneously ran a conventional GATK analysis procedure (using GRCh38 as the reference genome, using BWA-MEM for read alignment, retaining reads with an alignment quality value of not less than 20, using GATK HaplotypeCaller for variant detection, without any masking or physical link tracing) on ​​the above nine test samples. The detection results of the conventional GATK procedure were compared with the detection results of the present invention using the detection conclusion of Long-PCR combined with NGS sequencing as the gold standard. The comparison results are shown in Table 2, comparing the detection performance with the conventional GATK procedure. Table 2:

[0041] TP: Number of true positives; FN: Number of false negatives; FP: Number of false positives.

[0042] Comparative data shows that the conventional GATK procedure has an average recall rate of only 0.57 in the highly homologous CYP21A2 / CYP21A1P region, meaning that over 40% of the true pathogenic mutations are missed in the conventional procedure. This is mainly because a large number of valid reads are marked as low-quality and discarded due to multiple alignments, resulting in severely insufficient effective coverage depth of the target region. This invention, through the combined effect of Masker reference genome reconstruction and physical link tracing, achieved a 100% recall rate in all nine samples, representing an average improvement of 43 percentage points compared to the conventional GATK procedure.

[0043] The aforementioned performance improvement is not achieved by reducing specificity. The average precision of this invention is 0.90, significantly better than the average precision of 0.61 of the conventional GATK process on the same samples (the latter generates a large number of false positive variants due to the incorporation of pseudogene contamination reads). In terms of the overall performance index F1-score, this invention achieves an average of 0.95, while the conventional GATK process only reaches 0.59, an improvement of 0.36. This exceeds the performance gain that a person skilled in the art could reasonably expect based on existing technology improvements, constituting an unexpected technical effect.

[0044] Table 3 lists the detection performance data of the present invention on nine samples. The recall rate of each sample reached 1.00, and the F1-score was between 0.91 and 1.00, which comprehensively verified the robustness of the present invention in the detection of mutations in the high homology region of CYP21A2. Table 3:

[0045] Sample: Sample; TP: Count of true positives; FN: Count of false negatives; FP: Count of false positives; Recall: Recall rate; Precision: Precision rate; F1-score: F1 score.

[0046] It is worth noting that the aforementioned performance defects of the conventional GATK workflow stem from the inherent limitations of its algorithmic logic—this workflow treats each read as an independent statistical event, lacks targeted utilization of characteristic sites in paralogous genes, and also lacks the ability to trace mutations using the physical linking relationships between paired reads. Therefore, the performance defects of the conventional GATK workflow cannot be overcome through parameter tuning or simple workflow optimization.

[0047] In practice, based on the homologous region coordinates marked in the first module, the homologous region corresponding to the CYP21A1P pseudogene in the GRCh38 reference genome is hard-masked, that is, the reference sequence bases in the corresponding region are replaced with the character 'N', generating the Masker reference genome. The raw paired-end reads (PE150 mode) from the whole-exome sequencing of the sample are re-aligned to the Masker reference genome using BWA-MEM, generating a re-aligned BAM / CRAM file. GATK or an equivalent variant detection tool is then used to detect SNVs and Indels in the re-aligned file to obtain a candidate variant list.

[0048] In an alternative embodiment of the present invention, the above-mentioned masking process can also adopt a soft masking method, that is, converting the reference sequence of the homologous region of the CYP21A1P pseudogene into a lowercase letter representation, and using the difference in the default processing of uppercase and lowercase sequences by the alignment software to reduce the probability of read segment misalignment, which is suitable for computing environments where there are constraints on the modification of the reference genome.

[0049] After candidate variant detection is completed, the candidate variant list is classified and labeled according to the characteristic site fingerprint database established in the first module. Variants that are completely identical in location and allele type to those recorded in the characteristic site fingerprint database are defined as "characteristic site mutations." Variants located within the homologous regions of CYP21A2 and CYP21A1P but inconsistent with the characteristic site fingerprint database are defined as "de novo mutations." The true and false gene origins of these two types of variants are determined in the third module using corresponding attribution rules.

[0050] The third module involves determining the true and false gene origins of variants based on physical link information. Its basic logic is that read pairs generated by PE150 paired-end sequencing physically originate from the same DNA fragment; therefore, the allele status at any site within the same read pair is naturally in a cis relationship. By analyzing read pairs spanning candidate variant sites and characteristic sites, the known true and false gene identifiers at the characteristic sites can be used to infer the gene origin of the candidate variant sites, thus achieving the determination of variant origin.

[0051] For attribution determination of "novel mutations," reads covering both candidate variant sites and at least one characteristic site are screened. The number of supporting reads from eugenic alleles (wild-type) and pseudogenic alleles (mutant) at the covered characteristic sites is counted, and the proportion of each to the total covered reads at that characteristic site is calculated. When the proportion of eugenic alleles at the characteristic site exceeds the attribution threshold, the read is attributed to the CYP21A2 eugenic gene; when the proportion of pseudogenic alleles at the characteristic site exceeds the attribution threshold, the read is attributed to the CYP21A1P pseudogenic gene.

[0052] "Newly Occurring Site Mutation" Site Correction Figure 2 As shown: a. Does not contain this mutation, passes through at least one characteristic site, and has more than 50% wild-type at the characteristic site, and is a true wild-type; b. Contains the mutation, passes through at least one characteristic site, and the wild type at the characteristic site exceeds 50%, indicating a true mutant; c. Does not contain the mutation, passes through at least one characteristic site, and the characteristic site mutant is detected in more than 50% of cases; pseudogene wild type. d. Contains the mutation, passes through at least one characteristic site, and the characteristic site mutation is detected in more than 50% of cases, which is a pseudogene mutation; the above four classification scenarios are illustrated, where (a) to (d) correspond to the four determination results of true wild type, true mutant, pseudo wild type and pseudo mutant, respectively.

[0053] By combining the number of reads containing candidate mutations in the reads attributed to the true gene source, the true allele frequency and genotype of the mutation at the new site on the CYP21A2 true gene are determined; by combining the reads attributed to the pseudogene source, the corresponding situation of the mutation on the CYP21A1P pseudogene is determined.

[0054] For attribution determination of "characteristic site mutations," since the candidate mutation site itself is a characteristic site, it cannot be used as an attribution anchor. Therefore, it is necessary to screen reads that simultaneously cover the candidate characteristic site and at least another independent characteristic site. The allele status at the other characteristic site is used as the attribution basis. The determination rules are consistent with those for de novo mutations. The four attribution scenarios for characteristic site mutations are illustrated as follows: Figure 3 As shown, a. Does not contain this mutation, passes through at least another characteristic site, and the wild type at the characteristic site exceeds 50%, and the wild type is present in the true gene; b. Contains the mutation, passes through at least another characteristic site, and the wild type at the characteristic site exceeds 50%, indicating a mutant on the true gene; c. Contains the mutation, passes through at least another characteristic site, and the characteristic site mutant is detected in more than 50% of cases; pseudogene wild type. d. Does not contain this mutation, passes through at least another characteristic site, and the characteristic site mutation is detected in more than 50% of cases; pseudogene mutation. and Figure 2 The essential difference is that the source tracing anchor must be selected from another characteristic site besides the candidate mutation site, in order to avoid the circular dependence of that site on its own source tracing.

[0055] Regarding the setting of the attribution threshold, this invention sets the default value of the attribution threshold to 50%. That is, at a characteristic site, if the proportion of the number of supporting reads of a certain allele to the total number of reads covered by that site exceeds 50%, the read is determined to originate from the corresponding gene. The following provides a detailed explanation of the basis for selecting this threshold.

[0056] In a diploid genome, the theoretical read ratio of each of the two alleles at a normal heterozygous site is approximately 50%. However, the differentially expressed sites in a signature fingerprint database do not reflect the mutational heterozygosity of the sample, but rather the inherent sequence differences between true and false genes. For a normal diploid individual (two copies each of CYP21A2 and CYP21A1P), the read from the true gene theoretically carries the true gene allele, and the read from the false gene theoretically carries the false gene allele. Under ideal unbiased alignment conditions, the theoretical frequency of the true gene allele at the signature site is approximately 50%. In this invention, after Masker reference genome re-alignment, the reads are preferentially directed to the CYP21A2 true gene region, resulting in a systematically higher frequency of the true gene allele at the signature site than 50%, thus giving the 50% threshold good discriminative power in practice.

[0057] To verify the rationality of the 50% threshold, this invention analyzed 78 background samples with normal copy numbers (CYP21A2 CN=2), statistically analyzing the distribution of true gene allele frequencies at each characteristic locus. The results showed that the median allele frequency was concentrated in the 55% to 70% range, while the median pseudogene allele frequency was concentrated in the 30% to 45% range, with a clear distribution boundary between the two at 50%. Furthermore, this invention retrospectively analyzed positive control samples known to carry pathogenic CYP21A2 mutations (verified using the gold standard method of Long-PCR+NGS). Under the condition of setting the attribution threshold at 50%, the accuracy rate of attributing true gene variant sites reached 100%, and the accuracy rate of excluding pseudogene contamination sites reached over 97%, with no systematic misjudgments due to the threshold setting. Sensitivity analysis of the thresholds showed that when the attribution thresholds were set to 40%, 50%, and 60%, the attribution accuracy at all three thresholds was not significantly different in samples with a sequencing depth of at least 50×, remaining above 96%. In the low-depth region (sequencing depth below 30×), the 50% threshold exhibited a higher recall rate than the 60% threshold and a lower false positive rate than the 40% threshold, demonstrating optimal overall performance. Therefore, this invention sets 50% as the default attribution threshold. In practical applications, this threshold can be adjusted within the range of 40% to 70% based on the sample sequencing depth and quality control indicators. Figure 5 As shown, the correction effect of this module is intuitively verified with actual sample data: The three frames from top to bottom in the figure respectively show the characteristic VAF detected by the gold standard method (Long-PCR+NGS), the original VAF after Masker re-alignment of this invention, and the VAF after physical link tracing correction. The corrected results are highly consistent with the gold standard, verifying the accuracy of the attribution determination algorithm of this module.

[0058] The fourth module involves the quantitative assessment of the CYP21A2 eukaryotic copy number. This module establishes a quantitative estimation algorithm for the CYP21A2 eukaryotic copy number by integrating deep normalized data with allele frequency information at characteristic loci. The derivation process of the algorithm is explained in detail below.

[0059] Let the copy number of the CYP21A2 eukaryotic gene in the sample be . The copy number of the CYP21A1P pseudogene is Total number of copies At a specific feature locus in the feature locus fingerprint database, the frequency of pseudogene alleles... Defined as the proportion of reads carrying pseudogene alleles to the total number of reads effectively covered at that locus. Under ideal conditions of uniform sequencing depth, It should be equal to the proportion of pseudogene copy number to the total copy number, that is: ; Therefore, the CYP21A2 true gene copy number can be estimated using the following formula: ; Since the frequency estimation of a single feature site is affected by the fluctuation of sequencing depth, this invention uses the median of the pseudogene allele frequencies at all available feature sites in the feature site fingerprint database. As a robust estimator, replacing single-site frequencies, the final formula is: ; in, The sequencing depth data obtained through the following steps are as follows: First, the target region is divided into the realigned sequencing depth data. Based on the average depth of each exon capture interval within the target region, the data volume within the sample is normalized according to the total number of effective reads to eliminate differences in sequencing volume between samples. Then, the correlation between GC content and sequencing depth is corrected using a local weighted regression method to eliminate systematic errors introduced by GC bias. Finally, a background reference library is introduced for secondary correction. The sample depth after the above two corrections is compared with the depth distribution of the same region in the background library. After standardization using the Z-score method, the total copy number is obtained by linear conversion using the median copy number in the background library (theoretically 4 in normal diploid samples, i.e., 2 copies each of true and false genes) as the anchor point. The absolute estimate.

[0060] Regarding the construction of the background reference library, this invention has the following specifications to ensure the accuracy of copy number estimation.

[0061] Regarding sample selection, the background reference library should include no fewer than 30 samples, with an optimal sample size of 50 to 200. Included samples should meet the following criteria: First, the CYP21A2 copy number should be confirmed to be normal (i.e., CYP21A2 CN=2, CYP21A1P CN=2) by independent validation methods (including but not limited to MLPA, long-fragment PCR, or third-generation sequencing); second, the same DNA extraction method, library construction kit, and sequencing platform should be used as those used for the test samples; third, the average sequencing depth of the target region of the sample should be no less than 50×; fourth, the uniformity of exome capture region coverage (defined as the proportion of regions with a coverage depth of no less than 20% of the average depth to all target regions) should be no less than 85%.

[0062] In terms of depth processing, each sample included in the background library is normalized with the same amount of data as the sample to be tested and GC bias correction is performed in sequence. The corrected depth value of each target interval is calculated, and the median depth value of each interval in all background samples is taken to form a background depth reference matrix.

[0063] In terms of quality control assessment, after the background library is constructed, the coefficient of variation (CV) of each background sample in the target regions of CYP21A2 and CYP21A1P should be calculated. The CV value for each region should not exceed 0.15 to ensure the depth stability of the background library in the target region. When new samples are added or sequencing batches are changed, the background library should be updated regularly to maintain consistency between the background library and the test samples in terms of library construction conditions. Figure 8 As shown, the distribution of the estimated total copy number of each partition in the target region is presented after normalization, GC correction, and secondary correction of the background library for 78 background samples. The estimated values ​​of each partition are highly concentrated around CN=4, and the intra-batch standard deviation does not exceed 0.3, which verifies the construction quality of the background library and the stability of the depth correction algorithm, and provides a reliable reference benchmark for the copy number estimation of subsequent test samples.

[0064] Taking Embodiment 1 of the present invention as an example, in 78 background samples, after the above-mentioned normalization, GC correction, and secondary correction of the background library, the estimated total copy number of each partition of the target region is highly concentrated around CN=4, and the standard deviation of the estimated values ​​of each sample within the batch does not exceed 0.3, verifying the construction quality of the background library and the stability of the correction algorithm. Based on this background library, the estimated total copy number of the test sample in Embodiment 1 is: ; The feature site fingerprint database contains 36 available feature sites, and the median frequency of pseudogene alleles at each feature site is... Substitute into the above formula: ; The results indicate that the CYP21A2 eukaryotic copy number in Sample 1 is approximately 2, consistent with the actual situation, and there are no gene deletions or duplications. The estimated total copy number of the sample in Example 2 is... There are 35 available feature sites. Substitute into the formula: ; The results also indicated that the CYP21A2 true gene copy number was approximately 2, with no deletions or duplications, consistent with the gold standard results.

[0065] Figure 6 The results of a four-layer comparison of mutation frequencies of more than 30 characteristic loci fingerprint database sites in Example 1 are presented: from top to bottom, they are Long-PCR+NGS gold standard VAF, VAF after true gene correction of the present invention, original WES re-aligned full loci VAF, and VAF after full loci correction of the present invention. The correspondence of the four layers of data intuitively presents the difference in effect before and after correction of the present invention.

[0066] Table 4:

[0067] Table 4 lists detailed data for all detected mutation sites in Example 1, including the number of reference reads, the number of mutant reads, the inferred genotype, and the corresponding detection results of the gold standard method for each site in both the true and false genes. All true gene mutant genotypes detected by this invention are completely consistent with the gold standard conclusions, verifying the accuracy of the physical link tracing algorithm of this invention in actual samples.

[0068] Figure 7 This is a comparison chart of the four-layer mutation frequencies corresponding to Example 2, and its data structure is similar to... Figure 6 The consistency further verifies the stability and reproducibility of the results of this invention across different samples.

[0069] Table 5

[0070] Table 5 lists detailed data for all detected mutation sites in Example 2, with the data structure consistent with Table 3. Example 2 detected a total of 14 mutation sites. All true gene mutation genotypes were consistent with the gold standard conclusion, and the pseudogene allele read numbers were accurately separated after attribution correction, without producing any systematic false positives. This further verifies the reproducibility of this invention across samples with different mutation profiles.

[0071] like Figure 9 , Figure 10 As shown, the total copy number distribution of the nine test samples after background library correction is summarized and displayed. The estimated values ​​for each sample are all concentrated around CN=4, which is consistent with... Figure 7 The background sample distribution shown has good consistency, further verifying the stability and applicability of the copy number estimation algorithm of the present invention in batch samples.

[0072] Module 5 involves the detection of linked haplotype assembly and recombination events. The CYP21A2 gene is located in the RCCX module of the MHC class III region, which is a high-risk area for non-allelic homologous recombination (NAHR) due to its extremely high homology with the pseudogene CYP21A1P. Recombination events can lead to partial or complete pseudogenotyping of the CYP21A2 gene, which is one of the important molecular mechanisms of congenital adrenal hyperplasia. Accurate detection of such events has significant clinical value.

[0073] In practice, paired reads that span at least two feature sites are selected from the re-alignment files processed by the third module, and the allele status sequence of each read pair at each feature site covered is recorded. For each read pair, the feature site sequence it covers is compared with the CYP21A2 reference haplotype (all feature sites show true alleles) and the CYP21A1P reference haplotype (all feature sites show pseudoalleles). Using the local haplotype information carried by each read pair as the basic unit, the haplotypes are gradually extended and spliced ​​according to the position of the feature sites to construct a linked haplotype sequence covering the entire length of CYP21A2.

[0074] In recombination determination, if the assembled linked haplotypes exhibit different allele combinations before and after a certain characteristic locus interval—that is, the characteristic locus combination before the interval is consistent with the CYP21A2 eukaryotic reference haplotype, while the characteristic locus combination after the interval is consistent with the CYP21A1P pseudogene reference haplotype, or vice versa—then a true or pseudogene recombination event is determined to have occurred within that interval, and the adjacent characteristic locus intervals where the combinational shift occurred are marked as the estimated range of the recombination breakpoint. Figure 3 Taking the example shown, if haplotype assembly is performed using paired reads covering feature sites 1 to 4, the expected result is the haplotype "ATTT" consistent with the true gene. If the final assembly result shows an allelic type shift between feature sites 2 and 3, resulting in the haplotype "ATCA", then it is determined that a recombination event has occurred within the genomic region corresponding to feature sites 2 to 3. This region is then marked as the estimated range of the recombination breakpoint in the chr6 coordinate system. Figure 4 For example, by using paired reads r1, r2, and r3 covering feature sites 1 to 4 to assemble step by step, it is expected to assemble a haplotype 'ATTT' that is consistent with the true gene; if the assembly result changes to 'ATCA', it is determined that a recombination event has occurred in the genomic region between feature site 2 and feature site 3.

[0075] In Example 1, 15 uniformly covered feature loci from the feature locus fingerprint database were selected for haplotype assembly. The feature locus state consistent with the genotype reference haplotype was recorded as 1, and the state consistent with the pseudogene reference haplotype was recorded as 2. The final assembly result was "111112111111111", meaning that the 5th locus among the 15 feature loci was in a pseudogene allele state, while the rest were consistent with the genotype. Combined with the copy number result from the fourth module, it was determined that Sample 1 did not exhibit recombination. The state change of the 5th locus was related to the heterozygous mutation carried by that locus and did not constitute a basis for recombination determination. In Example 2, the assembly result of the same 15 feature loci was "111111111111111", and all feature loci were consistent with the genotype. Therefore, it was determined that Sample 2 did not exhibit recombination.

[0076] Based on the analysis results of modules three, four, and five, this invention can output the following complete CYP21A2 gene detection report for each sample: all SNV / Indel variant sites and their genotypes detected on the CYP21A2 genotype (expressed as the number of genotype covered reads, the number of genotype mutant reads, and allele frequencies); variant information at corresponding sites on the CYP21A1P pseudogene (used to help distinguish between gene conversion and true pathogenic mutations); quantitative estimation of the CYP21A2 genotype copy number and deletion / duplication determination conclusions; linkage haplotype sequence and recombination event detection conclusions, including the estimated recombination breakpoint interval (if any). The above comprehensive output provides comprehensive and accurate genomic information support for clinical genetic counseling and the molecular diagnosis of congenital adrenal hyperplasia.

[0077] The above-described specific embodiments are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Various modifications, substitutions, and improvements made by those skilled in the art to the technical solutions of the present invention based on the provided textual description and drawings, without departing from the design concept and spirit of the present invention, should all fall within the scope of protection of the present invention. The scope of protection of the present invention is determined by the claims.

Claims

1. A method for detecting CYP21A2 gene mutations, characterized in that, Includes the following steps: S1. By comparing with the reference genome, identify the differential sites and homologous regions between the CYP21A2 true gene and the CYP21A1P pseudogene, and construct a fingerprint database of paralogous characteristic sites. S2. Based on the homology interval, the region corresponding to the CYP21A1P pseudogene in the reference genome is masked to construct the Masker reference genome. The whole exome sequencing reads are re-aligned to the Masker reference genome to obtain a candidate variant list. S3. Using the physical link information of paired-end sequencing read pairs and the characteristic site fingerprint database, determine the true and false gene origins of each candidate variant; the determination of the true and false gene origins in step S3 includes the following steps: Step S31: Using the feature site fingerprint database constructed in step S1 as a reference, extract the specific alleles corresponding to the CYP21A2 true gene and the CYP21A1P pseudo gene at each feature site to form a set of molecular tracing anchor points. Step S32: For each candidate variant site, select the reads that can be used for tracing the source using the following two methods: Step S321: The alignment interval of the read segment itself simultaneously covers at least one anchor point in the candidate mutation site and the molecular origination anchor point set; Step S322: The alignment interval of the read itself covers the candidate variant site, and the alignment interval of its paired reads covers at least one anchor point in the molecular tracing anchor point set. Based on the physical co-origin relationship of the paired reads established by paired sequencing, the candidate variant and the anchor point are included in the same alignment coordinate range and associated. The method of step S321 is preferred; when the number of available read segments in step S321 is lower than the preset coverage depth threshold, the read segments screened by the method of step S322 are added. Step S33: Read the allele bases presented at the molecular tracing anchor point of the segment screened in step S32, and compare them one by one with the specific alleles of the CYP21A2 true gene and CYP21A1P pseudo gene in the characteristic site fingerprint database. When a segment only covers a single anchor point, the attribution is determined directly based on the matching result of that anchor point. When a segment covers two or more anchor points and the assignment conclusions of each anchor point are consistent, the assignment is determined based on the consistent conclusions. When the conclusions of the anchor points are inconsistent, the anchor point that supports the same conclusion shall be used if the number of anchor points is more than half of the total number of valid anchor points. If no one side has more than half the votes, the segment is marked as having uncertain ownership and is excluded from subsequent statistics; Based on the above results, the reads and the candidate variants they carry are assigned to either the CYP21A2 true gene or the CYP21A1P pseudogene. In step S3, after the molecular source tracing anchor set is established and before the physical link read screening is started, the candidate variant list is pre-classified based on the characteristic site fingerprint database, and differentiated molecular source tracing anchor strategies are selected for different categories: The genomic coordinates of the candidate variant match a differentially expressed site in the characteristic site fingerprint database, and the allele type is consistent with the known allele status of that differentially expressed site. When tracing the origin of such mutations, other characteristic sites besides the differential site are used as molecular tracing anchors, excluding the site itself from participating in the determination of the origin of the mutation. The genomic coordinates and allele types of the candidate variants were not found in the feature site fingerprint database; When tracing the origin of such mutations, one or more characteristic sites in the characteristic site fingerprint database are used as molecular origin tracing anchors; Apply the source attribution determination rules corresponding to their respective categories to the two types of mutations; S4. Based on the re-alignment sequencing depth, background reference library, and allele frequencies at characteristic sites, quantitatively estimate the CYP21A2 eukaryotic copy number; the estimation method for the CYP21A2 eukaryotic copy number in step S4 is as follows: The total copy number of CYP21A2 and CYP21A1P was obtained by performing in-sample data volume normalization, GC bias correction, and secondary correction based on the background reference library on the weight ratio depth. ; Extract the pseudogene allele frequencies at all available feature sites in the feature site fingerprint database, and use the median of the frequencies. As a robust estimator; The CYP21A2 true gene copy number can be calculated using the following formula. : ； S5. Based on the physical link information of paired-end sequencing reads spanning multiple characteristic sites, assemble linked haplotypes and determine whether recombination events have occurred between true and false genes.

2. The method for detecting CYP21A2 gene mutations according to claim 1, characterized in that, The reference genome self-alignment in step S1 is implemented as follows: the target locus region sequence is divided into overlapping sub-fragments using a sliding window algorithm, with a window length of 500 to 5000 bases and an overlap step size of 10% to 90% of the window length; the sub-fragments are then compared back to the reference genome to screen for fragments with high-resolution alignments at multiple positions. Multiple sequence alignments are performed on the identified homologous regions to extract stable differential sites between CYP21A2 and CYP21A1P, forming the feature site fingerprint database.

3. The method for detecting CYP21A2 gene mutations according to claim 1, characterized in that, The masking process described in step S2 is a hard masking process that replaces the sequence of the CYP21A1P pseudogene homologous region in the reference genome with the character 'N', so that the reads that were originally ambiguous due to multiple alignments are now uniquely aligned to the CYP21A2 true gene region.

4. The method for detecting CYP21A2 gene mutations according to claim 1, characterized in that, The construction of the background reference library satisfies the following conditions: The number of samples included was no less than 30, and the CYP21A2 copy number of the samples was confirmed to be normal by an independent validation method. The same DNA extraction method, library construction kit, and sequencing platform were used for both the included samples and the samples to be tested. The average sequencing depth of the target region in the included samples was no less than 50×, and the exome capture uniformity was no less than 85%. After the background library is constructed, the depth variation coefficient of each target interval in all background samples does not exceed 0.15; When changing sequencing batches, the background library is updated on a rolling basis to maintain consistency with the construction conditions of the sample library to be tested.

5. The method for detecting CYP21A2 gene mutations according to claim 1, characterized in that, The method for determining the assembly and recombination of linked monomers in step S5 is as follows: Select paired reads that span at least two characteristic sites, record the allele status at each site in order of the characteristic site position, and gradually splice and extend them to form a linked haplotype sequence covering the full length of CYP21A2. The obtained linked haplotype sequences were compared position by position with the CYP21A2 true gene reference haplotype and the CYP21A1P pseudo gene reference haplotype; If a linked haplotype exhibits allele combinations from different sources before and after a certain characteristic locus interval, i.e., the characteristic locus combination before the interval is consistent with the CYP21A2 true gene reference haplotype while the one after the interval is consistent with the CYP21A1P pseudo gene reference haplotype, or vice versa, then it is determined that a true or pseudo gene recombination event has occurred within that interval, and the genomic interval corresponding to the adjacent characteristic locus where the recombination change occurred is marked as the estimated range of the recombination breakpoint.

6. A CYP21A2 gene mutation detection system, applicable to the CYP21A2 gene mutation detection method according to any one of claims 1-5, characterized in that, include: The feature site fingerprint database construction module is used to identify differential sites and homologous regions between the CYP21A2 eukaryote and the CYP21A1P pseudogene by self-alignment with the reference genome, and to construct a paralogous feature site fingerprint database. The alignment reconstruction and variant detection module is used to mask the region corresponding to the CYP21A1P pseudogene in the reference genome based on the homology interval, construct the Masker reference genome, re-align the whole exome sequencing reads to the Masker reference genome, and obtain a list of candidate variants. The physical link attribution determination module is used to determine the true and false gene origins of each candidate variant by utilizing the physical link information of paired sequencing read pairs and the feature site fingerprint database. The copy number assessment module is used to quantitatively estimate the CYP21A2 eukaryotic copy number based on the re-alignment sequencing depth, background reference library, and allele frequency of characteristic sites. The haplotype assembly and recombination detection module is used to assemble linked haplotypes based on the physical link information of paired-end sequencing reads spanning multiple characteristic sites, and to determine whether recombination events have occurred between true and false genes.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1-5.

Images