A screening method of a variety SNP marker, a computer device and a readable storage medium

Abstract

This invention discloses a method, computer device, and readable storage medium for screening SNP markers of plant varieties, relating to the field of bioinformatics. The method for screening SNP markers of plant varieties provided by this invention can automatically and efficiently screen a set of core single nucleotide polymorphism (SNP) sites with the smallest number and strongest discriminative ability from large-scale high-throughput sequencing variant site data (VCF), for the identification of the plant variety to which a sample belongs.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, and more specifically, to a method for screening SNP markers for a variety, a computer device, and a readable storage medium. Background Technology

[0002] In fields such as plant and animal variety protection, seed purity identification, variety rights protection, germplasm resource management, and judicial identification, a reliable and efficient molecular marker system is needed to uniquely identify varieties. With the popularization of high-throughput sequencing technology, population resequencing can obtain millions of SNP loci, but how to select a set of "core SNPs" from them—that can cover the main regions of the genome and distinguish all varieties with as few loci as possible—is the core technical challenge in this field.

[0003] One method for variety differentiation is based on fixed molecular marker identification (SSR marker method): using developed simple sequence repeat (SSR) microsatellite markers for variety identification, and detecting fragment length polymorphism by capillary electrophoresis. Disadvantages include: limited number of markers, uneven genome coverage; long and costly development cycle for new markers; and poor reproducibility between different laboratories.

[0004] Another method for species differentiation is manual marker screening based on whole-genome resequencing: researchers manually select highly polymorphic loci from whole-genome variation data for identification based on experience. Disadvantages include: strong subjectivity and lack of standardized procedures; difficulty in ensuring uniform distribution of loci across the genome; inability to quantitatively assess the overall discriminative power of the selected locus set; and poor adaptability, with results varying significantly across different species.

[0005] In summary, existing technologies rely on manual experience to select marker sites, lack systematic automatic screening strategies, and struggle to balance the polymorphism information of sites, uniform genome coverage, and sequence uniqueness. Furthermore, they cannot automatically generate core markers for rapid traceability of varieties.

[0006] In view of this, the present invention is proposed. Summary of the Invention

[0007] The purpose of this invention is to provide a method, computer device, and readable storage medium for screening SNP markers of varieties to uniquely distinguish all varieties in a population using the fewest possible SNP loci.

[0008] This invention is implemented as follows: In a first aspect, the present invention provides a method for screening SNP markers of a variety, comprising the following steps: S1: Adaptive filtering of variant sites: The following quality indicators for variant sites are calculated manually or automatically: site quality value, deletion rate and sequencing depth. The filtering threshold for each indicator is determined. Variant sites are filtered according to at least one specified quality indicator, and sites that pass the threshold for all specified indicators are retained as a candidate site set. S2: Sort the candidate site set: The candidate loci were ranked using genetic diversity indices; the genetic diversity indices were selected from at least one of the following: homozygosity product, genotype diversity index, nucleotide polymorphism Pi, expected heterozygosity He, polymorphism information content PIC, and Shannon information index. S3: Chromosome Round-Robin Sort: Within a set of candidate loci with the same genetic diversity index, a chromosome polling mechanism is used to arrange them, so that the candidate loci preferably come from different chromosomes alternately, thereby improving the uniformity of the distribution of loci on the genome. This allows the candidate loci to enter the alternately arranged candidate queue for sorting. When the available loci on some chromosomes are exhausted in advance, the chromosomes that are exhausted in advance are "skipped" in the polling. The adjacent loci from the remaining available loci on the same chromosome are used as candidate loci and enter the alternately arranged candidate queue for sorting. S4: Incremental evaluation of discriminative ability and site filling with adaptive weight reduction: Establish a genotype fingerprint dictionary for each variety. Initially, the fingerprint of each variety is an empty string. Evaluate each candidate site in the order sorted in step S3. For each candidate site, append its genotype information to the fingerprint string of each variety. Calculate the number of varieties that can be uniquely distinguished after updating the genotype fingerprint dictionary. According to the site inclusion rules, include sites that can increase the distinguishability of varieties into the genotype fingerprint dictionary. If all varieties are still not distinguished after traversing all candidate sites under the current weight coefficient, reduce the weight coefficient and use the new weight coefficient to traverse the candidate sites again. S5: Removal of multiple copies or highly homogeneous regions.

[0009] Secondly, the present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the above-mentioned method for screening variety SNP markers.

[0010] Thirdly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-described method for screening variety SNP markers.

[0011] The present invention has the following beneficial effects: The method for screening SNP markers of varieties in this invention significantly reduces the number of loci obtained: by using incremental evaluation of discriminative ability and an adaptive weight reduction strategy, the set of loci screened by this invention has minimal redundancy, achieving unique identification of all varieties with the fewest number of loci, and reducing the cost of subsequent molecular detection.

[0012] Fully automated site selection: Compared to manually selecting marker sites or simple filtering based on a single threshold, this invention achieves a fully automated process from quality filtering to site selection, reducing human intervention and subjective bias.

[0013] Uniform genome coverage: The chromosome polling mechanism ensures that core SNPs are evenly distributed across the genome, avoiding the concentration of markers in a few chromosomal regions and improving the genome-wide representativeness of the identification results.

[0014] Exclude nonspecific sites: By removing sites with multiple copies or highly homologous regions, the specificity and accuracy of the selected sites in actual molecular detection are ensured. Attached Figure Description

[0015] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 This is a flowchart illustrating the technical path of the method for screening SNP markers for varieties according to the present invention. Figure 2 Selecting density maps for screening methods of core markers of varieties; Figure 3 Genotype QR code images selected for the core marker screening method for varieties; Figure 4 This is a schematic diagram of the physical structure of an electronic device; Figure 5 This is a flowchart illustrating the three-level screening strategy of the SNP marker screening method for varieties according to the present invention. Detailed Implementation

[0017] Reference will now be made to detailed embodiments of the present invention, one or more of which are described below. Each example is provided for explanation and not for limitation of the invention. In fact, it will be apparent to those skilled in the art that various modifications and variations can be made to the invention without departing from its scope or spirit. For example, features described or illustrated as part of one embodiment may be used in another embodiment to produce further embodiments.

[0018] In a first aspect, the present invention provides a method for screening SNP markers of a variety, comprising the following steps: S1: Adaptive filtering of variant sites: The following quality indicators for variant sites are calculated manually or automatically: site quality value, deletion rate and sequencing depth. The filtering threshold for each indicator is determined. Variant sites are filtered according to at least one specified quality indicator, and sites that pass the threshold for all specified indicators are retained as a candidate site set. S2: Sort the candidate site set: The candidate loci were ranked using genetic diversity indices; the genetic diversity indices were selected from at least one of the following: homozygosity product, genotype diversity index, nucleotide polymorphism Pi, expected heterozygosity He, polymorphism information content PIC, and Shannon information index. S3: Chromosome Round-Robin Sort: Within a set of candidate loci with the same genetic diversity indicators, a chromosome round-robin mechanism is used to arrange them, so that candidate loci preferably come from different chromosomes to improve the uniformity of locus distribution on the genome. This allows candidate loci to enter an alternately arranged candidate queue for sorting. When the available loci on some chromosomes are exhausted in advance, the chromosomes that are exhausted in advance are "skipped" in the round-robin process. The adjacent loci among the remaining available loci from the same chromosome are used as candidate loci and enter an alternately arranged candidate queue for sorting. At this time, it is possible for adjacent candidate loci to come from the same chromosome. S4: Incremental evaluation of discriminative ability and site filling with adaptive weight reduction: Establish a genotype fingerprint dictionary for each variety. Initially, the fingerprint of each variety is an empty string. Evaluate each candidate site in the order sorted in step S3. For each candidate site, append its genotype information to the fingerprint string of each variety. Calculate the number of varieties that can be uniquely distinguished after updating the genotype fingerprint dictionary. According to the site inclusion rules, include sites that can increase the distinguishability of varieties into the genotype fingerprint dictionary. If all varieties are still not distinguished after traversing all candidate sites under the current weight coefficient, reduce the weight coefficient and use the new weight coefficient to traverse the candidate sites again. S5: Removal of multiple copies or highly homogeneous regions.

[0019] Step S1 uses adaptive filtering to retain sites that meet preset thresholds for specified indicators, filtering out outliers for each indicator, such as sites exceeding the outlier boundary. Traditional methods typically perform threshold filtering based on a single indicator (such as polymorphism information content, PIC), failing to comprehensively consider multi-dimensional quality indicators such as site quality, deletion rate, sequencing depth, and polymorphism information content. Furthermore, the filtering threshold relies on manual setting and lacks adaptability. This invention utilizes statistical methods such as IQR to automatically calculate the filtering threshold based on the data's inherent distribution characteristics, eliminating the need for manual empirical parameters. It also supports manual overriding of specific indicators, balancing automation and flexibility. Compared to manually selecting marker sites or simple filtering based on a single threshold, this invention achieves a fully automated process from quality filtering to site selection, reducing human intervention and subjective bias. It also allows users to manually specify thresholds for specific indicators to override automatically calculated values. Sites are jointly filtered according to specified indicator combinations, retaining only sites where all indicators pass the threshold as the candidate site set.

[0020] In existing technologies, locus ranking indicators are disconnected from actual distinguishing needs. Current methods typically use traditional population genetic indicators such as nucleotide polymorphism (Pi) or polymorphism information content (PIC) to prioritize loci. However, these indicators reflect allele-level diversity and do not directly measure the actual contribution of a locus to variety differentiation. For example, a locus with extremely high heterozygosity but rare homozygous genotypes may have a high Pi value, but since most varieties exhibit heterozygosity at this locus, its actual distinguishing ability is extremely low.

[0021] In step S2 of this invention, when the genetic diversity index is selected from homozygous product and genotype diversity index, a composite score is constructed with homozygous product as the primary ranking key and genotype diversity as the secondary ranking key. Sampling is performed by chromosome round-robin at the same score level to ensure that loci with high variety differentiation potential are prioritized for evaluation and that genome coverage is uniform. This ensures the organic unity between genome coverage uniformity and locus differentiation potential ranking.

[0022] The method for screening variety SNP markers provided by this invention is implemented at three levels: First, in step S2, at the genetic diversity index ranking level, candidate loci are prioritized according to the composite score of each locus, so that loci with higher scores enter the screening process first; wherein, in a preferred embodiment, the composite score uses the homozygous product as the primary ranking basis and the genotype diversity index as the secondary ranking basis.

[0023] Secondly, in step S3, at the chromosome polling and sorting level, for candidate sites with the same composite score, queues are formed according to their respective chromosomes, and a chromosome polling mechanism is adopted to allow candidate sites at the same score level to enter the subsequent evaluation process in a cross-chromosome alternation manner, so as to improve the uniformity of the final results distribution on the genome.

[0024] Finally, in step S4, in the adaptive weighted decreasing discrimination ability increment assessment and site filling level, for each site that enters the assessment in the aforementioned order, combined with the genotype fingerprint dictionary or sample grouping status formed by the currently selected sites, the number of newly distinguishable varieties after the site is added is calculated, and the site is determined to be included in the core marker set according to the preset discrimination increment rules.

[0025] The flowchart is shown below. Figure 5 As shown.

[0026] The example is shown below: Assuming that after filtering in step S1, 6 candidate loci are obtained, their composite scores and chromosome distributions are as follows: | Candidate site | Composite score | Chromosome location | |---|---:|---| | P1 | 0.240000666666 | Chr1 | | P2 | 0.240000666666 | Chr1 | | P3 | 0.240000666666 | Chr2 | | P4 | 0.240000666666 | Chr3 | | P5 | 0.210000555555 | Chr2 | | P6 | 0.180000444444 | Chr1 | Among them, the composite score preferably uses the homozygous product as the primary ranking criterion and the genotype diversity index as the secondary ranking criterion.

[0027] Step S2: Ranking hierarchy of genetic diversity indicators: First, candidate loci are processed according to their composite scores from high to low. Therefore, the locus group with a composite score of `0.240000666666` is processed first. - P1 (Chr1); - P2 (Chr1); - P3 (Chr2); - P4 (Chr3); Then process the lower-scoring P5 and P6.

[0028] Step S3: Chromosome polling and sorting hierarchy For locus groups with the same composite score, they are then grouped by chromosome, for example: - Chr1 queue: P1, P2; - Chr2 queue: P3; - Chr3 queue: P4; Assuming the preset chromosome order is `Chr1 → Chr2 → Chr3`, the order of the sites to be evaluated after chromosome round-robin sorting is as follows: - Round 1: P1 → P3 → P4; - Round 2: P2 (Chr2 and Chr3 have no remaining sites, so skip them); Therefore, the order of the sites to be evaluated in this scoring group is: >**P1 → P3 → P4 → P2**.

[0029] The purpose of this step is to prioritize candidate loci in the subsequent evaluation process under the same composite scoring conditions, in a cross-chromosomal alternation manner, rather than having P1 and P2 on Chr1 processed consecutively first.

[0030] Step S4: Incremental evaluation of discriminative ability with adaptive weight reduction and site filling hierarchy Assume there are currently 8 varieties, and initially the sample grouping status corresponding to the variety genotype fingerprint dictionary is group 1.

[0031] (1) First evaluate P1 With the addition of P1, the original 1 group of 8 varieties can be split into 3 groups.

[0032] If, under the current weighting stage, the number of newly added distinctions meets the inclusion rules, then: - P1 is included in the core tag set - Variety genotype fingerprint dictionary updated; - The current sample grouping status has changed from group 1 to group 3.

[0033] (2) Reassess P3: Based on the inclusion of P1, continue to evaluate P3.

[0034] Assuming that after P3 is added, the sample grouping status remains 3 groups, and no new distinguishable varieties are added, then: - P3 is not included in the core tag set; - The genotype fingerprint dictionary for the variety is not updated; - The current sample grouping status remains at 3 groups; (3) Reassess P4: Continue to evaluate P4.

[0035] Assuming that after P4 is added, the current 3 groups can be further divided into 5 groups, and the inclusion rules under the current weight stage are satisfied, then: - P4 Proceed to the homology / uniqueness verification step; If P4 passes the homology / uniqueness verification, then: - P4 is included in the core tag set; - Variety genotype fingerprint dictionary updated; - The current sample grouping status has changed from 3 groups to 5 groups. (4) Reassess P2 Continue to evaluate P2.

[0036] Assuming that the addition of P2 can further divide the current 5 groups into 6 groups, meeting the requirements for incremental discrimination ability assessment; however, its flanking sequences were found to have multiple high similarity hits after whole-genome alignment, failing the homology / uniqueness verification, then: - P2 is not included in the core tag set; - The current sample grouping status remains unchanged at 5 groups.

[0037] Proceed to the next rating level: After completing the assessment of the locus group with a composite score of `0.240000666666`, proceed to process the next scoring level, P5 and P6, and repeat the above process.

[0038] In other embodiments, the genetic diversity index is selected from a composite score of homozygous product-nucleotide polymorphism (Pi), or a composite score of homozygous product-expected heterozygosity (He), or a composite score of homozygous product-polymorphism information content (PIC) or Shannon information index. Sampling by chromosome round-robin at the same scoring level can also, to some extent, achieve priority assessment of loci with high variety differentiation potential and uniform genome coverage. The preferred genetic diversity index is a composite score of homozygous product and genotype diversity index.

[0039] Step S3, chromosome round-robin sorting, aims to sort the sets of candidate loci with the same genetic diversity indexes by placing those with unevenly distributed loci on the genome at the back and those with evenly distributed loci at the front. This ensures that, in the case of the same score, adjacent candidate loci come from different chromosomes, thereby guaranteeing the uniform distribution of the final screening results on the genome and avoiding the clustering of loci in a few highly variable regions.

[0040] The core idea of step S4 is: in the initial stage, rapidly accumulate loci with high discrimination efficiency using strict criteria, and in the later stage, gradually relax the inclusion criteria to supplement loci that are effective for a small number of difficult-to-distinguish varieties, so as to automatically balance the contradiction between "the minimum number of loci" and "the maximum discrimination capability" during the iterative process. The weight coefficient not only controls the stringency of the dynamic threshold, but also realizes "progressive relaxation of the overall screening strategy" by affecting the starting boundary of the later stage, and finally achieves approximately optimal selection of the minimum locus set, thereby balancing the minimum number of loci and the maximum discrimination capability. However, existing methods do not consider the redundancy of discrimination information between loci when screening loci, that is, whether adding a new locus can further increase the number of distinguishable varieties on the basis of the existing loci. This leads to a large number of redundant loci in the selected locus set, and more markers are required under the same discrimination capability. The present invention significantly reduces the number of markers under the same discrimination capability.

[0041] The prior art does not exclude multi-copy / homologous sequence regions: the loci selected by the prior art may be located in repeated regions or multi-copy regions of the genome, which will cause false positives due to non-specific amplification in actual molecular detection (such as primer design and probe hybridization). In the present invention, the sequence uniqueness check at the bioinformatics level is embedded into the locus screening loop through step S5 to exclude loci in multi-copy regions, so as to ensure the specificity and accuracy of the selected loci in actual molecular detection.

[0042] In a preferred embodiment of the present invention, the locus inclusion rules in step S4 include: Let the number of varieties that have been distinguished by the current candidate locus set be D, the total number of varieties be N, the number of distinguishable varieties after adding a new locus be D_new, the current weight coefficient be W, the dynamic threshold M = max(1, W × N), the stage dividing line S1 = 0.1 × N, S2 = 0.9 × N + 0.1 × W × N; S1 refers to the discrimination dividing line between the initial stage and the middle stage, and S2 refers to the discrimination dividing line between the middle stage and the later stage.

[0043] 0.1 in the above stage dividing lines is an empirical value, under which fewer SNPs are obtained. 0.9 is also an empirical value, under which more SNP loci that can improve the discrimination can be included.

[0044] In the initial stage, where D < S1: when D is 0, if D_new - D ≥ 2, this locus can be included in the variety genotype fingerprint dictionary; when D > 0, any locus that can increase at least one distinguishable variety is included in the variety genotype fingerprint dictionary; In the middle stage, where S1 ≤ D ≤ S2: if D_new - D > M, this locus can be included in the variety genotype fingerprint dictionary; In the later stage, when D>S2: any locus that can add at least one distinguishable variety is included in the variety genotype fingerprint dictionary.

[0045] In a preferred embodiment of the present invention, W is initially set to 0.7 (empirical value). If all candidate sites are not distinguished under the current weight coefficient, the weight coefficient is reduced and a new weight coefficient is used to re-traverse the candidate sites. In the later stage, the inclusion criteria are gradually relaxed to supplement sites that are effective for a few difficult-to-distinguish varieties. W is not initially set to 1, considering that no sites will be included when W is 1.

[0046] In a preferred embodiment of the present invention, the weighting coefficient W is gradually reduced from 0.7 to 0.

[0047] The weighting coefficient W not only controls the strictness of the dynamic threshold, but also achieves a "gradual relaxation of the overall screening strategy" by influencing the starting boundary of the later stage, ultimately achieving a near-optimal minimum site set selection.

[0048] In a preferred embodiment of the present invention, if all candidate sites in the candidate site set fail the incremental evaluation of discriminative ability, the added genotype information is removed from the fingerprint strings of all varieties from the variety genotype fingerprint dictionary.

[0049] In a preferred embodiment of the present invention, step S5, removal of multiple copies or highly homologous regions, includes: extracting flanking sequences on the reference genome for each candidate site that has passed the discrimination ability assessment, and constructing local sequence fragments; comparing the local sequence fragments with the whole genome for sequence similarity, and removing candidate sites with multiple highly similar alignment positions; In a preferred embodiment of the present invention, the criterion for determining the existence of multiple highly similar alignment positions is as follows: if the similarity of the second similar (i.e., the second hit) sequence alignment is greater than or equal to a preset threshold and the alignment length is greater than or equal to a preset proportion of the total length of the flanking sequences, then it is determined that there are multiple highly similar alignment positions.

[0050] In a preferred embodiment of the present invention, after step S5 is completed, steps S4 and S5 are executed in multiple rounds in sequence. The screening results of each round are automatically added to the "known sites" and "skipped sites" lists of the next round, so that subsequent rounds can expand from different candidate sites to obtain the core SNP site set.

[0051] To improve the robustness of the screening results, multiple rounds of iteration can be performed on the screening process in steps S4 and S5 (e.g., 2, 3, 4, 5, 6, 7, 8, 9, or 10 rounds). The screening results of each round are automatically added to the "known sites" and "skipped sites" lists of the next round, allowing subsequent rounds to expand from different candidate sites and accumulate a richer set of core SNP sites.

[0052] In a preferred embodiment of the present invention, in step S2, the genetic diversity index is selected from the homozygous product and the genotype diversity index, and the two are combined and scored; the composite score is constructed with the homozygous product as the primary sorting key and the genotype diversity index as the secondary sorting key. The formula for calculating the homozygous product is: HP = f_ref × f_alt; f_ref refers to the frequency of the homozygous reference genotype in the sample, and f_alt refers to the frequency of the homozygous variant genotype in the sample; The formula for calculating the genotype diversity index GD is: GD = 1 - Σ(fi²); where fi is the frequency of the i-th genotype; The scoring method for the composite score includes: retaining 3-8 decimal places for the HP value as the high digit, and concatenating 4-8 decimal places of the GD value as the low digit to form a numerical composite score. First, the scores are sorted in descending order of HP, and if HP is the same, they are sorted in descending order of GD.

[0053] A composite score is constructed using the homozygous product (HP = f_ref × f_alt) as the primary sorting key and the genotype diversity index (GD = 1 - Σ(fi²)) as the secondary sorting key. This method breaks through the traditional approach of using nucleotide polymorphism (Pi) or minimum allele frequency (MAF) for locus sorting, directly focusing on the most valuable homozygous genotype oppositions in variety identification scenarios. This improves the quality of the candidate locus cohort from the source, resulting in a smaller number of core loci selected and higher differentiation efficiency.

[0054] In a preferred embodiment of the present invention, in step S3, if the scoring levels are the same, the candidate sites are first grouped according to their respective chromosomes, and the maximum number of candidate sites on each chromosome is used as the number of iterations. In each round, one site is selected from each chromosome that still has remaining candidate sites, forming a candidate queue that preferentially alternates across chromosomes. At the same time, at the chromosome polling and sorting level, a single chromosome selection upper limit is set based on the number of currently selected sites on each chromosome. Before a candidate site enters the incremental evaluation of discriminative ability, it is first determined whether its chromosome has reached the current upper limit. If it has, the candidate site is skipped in the current upper limit round. After the current upper limit round has scanned the entire candidate site list, the single chromosome selection upper limit is increased, and the candidate site list is scanned again from high score to low score, so that the candidate sites that were previously skipped at this scoring level have another chance to enter the incremental evaluation of discriminative ability.

[0055] Loci with the same scoring level are grouped according to their respective chromosomes, and the number of loci on each chromosome is recorded. This can be done using tools such as Excel or others.

[0056] If the score levels are the same, the maximum number of loci is used as the number of iterations. In each round, one locus is sequentially selected from each chromosome, forming an alternating candidate queue. Thus, within any set of candidate loci, adjacent candidate loci come from different chromosomes. For example, if chr1 has 10 loci with the same score level, chr2 has 7 loci with the same score level, and chr3 has 1 locus with the same score level, then the loop is repeated 10 times. The first iteration repeats the loop for locus 1 of chr1, locus 1 of chr2, and locus 1 of chr3. The second iteration repeats the loop for locus 2 of chr1, locus 2 of chr2, and chr3 is skipped, and so on.

[0057] In a preferred embodiment of the present invention, in step S1, the filtering threshold for the deletion rate is automatically determined by using the interquartile range statistical method; the filtering threshold for the sequencing depth is determined by the average value of the sequencing depth; and the filtering threshold for the site quality value is determined by taking half of the average quality value of all sites. In a preferred embodiment of the present invention, the average depth threshold is set to 0.5 to 2 times the mean; sites with an average depth within this range are included in the candidate site set. The missing rate threshold is set to the upper bound of the IQR; sites with a missing rate not exceeding the upper bound of the IQR are included in the candidate site set, while those exceeding it are removed. Sites with a quality value not lower than the filtering threshold (i.e., the site quality threshold) are included in the candidate site set; sites below this threshold are removed. Sites with an average depth lower than the lower limit or higher than the upper limit of the threshold are considered for removal.

[0058] Calculate the first quartile (Q1) and third quartile (Q3) of each indicator. The outlier boundary is defined as [Q1-1.5×IQR, Q3+1.5×IQR], where IQR=Q3-Q1.

[0059] In the main embodiment of the invention, the number of loci is automatically determined with the goal of distinguishing all varieties, ultimately obtaining the smallest subset of loci, which is smaller than the specifications of common detection chips. In other embodiments, a fixed upper limit for the number of loci can be preset (e.g., 500 or 1000, corresponding to common detection chip specifications), and the maximum subset of loci is obtained through multiple rounds of execution under this constraint.

[0060] The core SNP markers of a variety obtained by the screening method for variety SNP markers provided by this invention can be used for: (1) Seed purity detection: The consistency between seed batches and standard varieties is compared using core SNP fingerprints to quickly identify seed purity and authenticity.

[0061] (2) Kinship analysis: Parentage testing or family analysis is conducted by comparing the similarity of molecular fingerprints between different individuals.

[0062] (3) Judicial identification for plant variety rights protection: Molecular fingerprint QR codes can be used as molecular evidence of variety uniqueness for legal identification in plant variety rights infringement cases.

[0063] (4) Germplasm resource database construction: Generate standardized molecular fingerprint codes for each material in the germplasm resource bank to facilitate the classification management and rapid retrieval of resources.

[0064] Secondly, the present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the above-mentioned method for screening variety SNP markers.

[0065] Specifically, the electronic device may include a memory, a processor, a bus, and a communication interface, which are electrically connected directly or indirectly to each other to enable data transmission or interaction. For example, these components may be electrically connected to each other via one or more buses or signal lines. The processor may process information and / or data related to target identification to perform one or more functions described in this application.

[0066] The memory can be, but is not limited to, Random Access Memory (RAM), Read Only Memory (ROM), Programmable Read-Only Memory (PROM), Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read-Only Memory (EEPROM), etc.

[0067] A processor can be an integrated circuit chip with signal processing capabilities. This processor can be a general-purpose processor, including a Central Processing Unit (CPU), a Network Processor (NP), etc.; it can also be a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), a Field-Programmable Gate Array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0068] Thirdly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-described method for screening variety SNP markers.

[0069] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below. Where specific conditions are not specified in the embodiments, conventional conditions or conditions recommended by the manufacturer shall apply. Reagents or instruments whose manufacturers are not specified are all conventional products that can be purchased commercially.

[0070] The features and performance of the present invention will be further described in detail below with reference to embodiments.

[0071] Example 1 This embodiment provides a method for screening core SNP markers of a variety, the flowchart of which is shown below. Figure 1 As shown, it includes the following steps: Step 1: Adaptive Filtering. Multidimensional Quality Calculation and Adaptive Filtering of Variant Sites. First, obtain the variant detection result file (VCF format) from the population high-throughput whole-genome resequencing data. For each SNP site in the file, calculate the following quality metrics: Site quality score (Qual): Reflects the reliability of variant detection; Miss Rate: The proportion of genotypes missing at this locus across all varieties; Average sequencing depth (AvgDepth): The average coverage depth of all varieties at this site; For the above-mentioned locus quality values ​​and missing rate indicators, the interquartile range (IQR) statistical method was used to automatically determine the filtering threshold. Specifically, the first quartile (Q1) and the third quartile (Q3) of each indicator were calculated, and the outlier boundary was defined as [Q1-1.5×IQR, Q3+1.5×IQR], where IQR=Q3-Q1.

[0072] The thresholds for each indicator are determined as follows: the site quality threshold is half of the mean quality of all sites; the missing rate threshold is the upper bound of the IQR; and the average depth threshold is in the range of 0.5 to 2 times the mean.

[0073] Loci for which all metrics pass the threshold are retained as candidate loci. Loci with a missing rate not exceeding the upper bound of the IQR are included in the candidate loci set; those with a missing rate exceeding the IQR upper bound are removed. Loci with a quality value not lower than the filtering threshold (i.e., the loci quality threshold) are included in the candidate loci set; those with a quality value lower than the loci quality threshold are removed; and loci with an average depth lower than the lower limit of the threshold or higher than the upper limit of the threshold are considered for removal.

[0074] It also supports users manually specifying thresholds for specific indicators to override automatically calculated values. It performs joint filtering on loci based on specified indicator combinations, retaining only loci for which all indicators pass the threshold as the candidate locus set.

[0075] Step 2: Calculation of Homozygous Product-Genogenetic Diversity Composite Score and Chromosome Round-Robin Sort. For the filtered candidate loci, calculate the homozygous product-genotype diversity composite score for each locus, as follows: (1) Genotype statistics and standardization: Read the genotype information of all varieties at each locus, skip the missing genotypes, and standardize the genotypes (sort the alleles by number and unify the phase representation so that the same genotype with different writing orders is regarded as the same genotype, for example, `0 / 1`, `1 / 0`, `0|1`, `1|0`, after standardization they are all regarded as the same genotype), and count the number of varieties for each genotype.

[0076] (2) Homozygous product (HP) calculation: Calculate the frequency f_ref of the homozygous reference genotype (e.g., 0 / 0) and the frequency f_alt of the homozygous variant genotype (e.g., 1 / 1) in the sample. The homozygous product HP = f_ref × f_alt. The design principle of this index is: the core of variety differentiation lies in the fact that different varieties exhibit different homozygous genotypes at this locus. When the frequencies of both homozygous genotypes are high and relatively balanced, the HP value is the largest, meaning that this locus is most likely to produce stable genotypic differences among different varieties.

[0077] Compared to traditional nucleotide polymorphisms (Pi) or minimum allele frequency (MAF), the homozygous product directly focuses on the most valuable homozygous genotype oppositions in variety identification, avoiding the overvaluation of sites with high heterozygosity but insufficient homozygous differentiation.

[0078] (3) Genotype diversity index (GD) calculation: Calculate the frequency distribution of all occurrences of genotypes (including homozygous and heterozygous), where genotype diversity GD = 1 - Σ(fi²), and fi is the frequency of the i-th genotype. This index is based on the Simpson diversity index and reflects the richness of genotype composition of a locus in the population, serving as a supplementary ranking basis for the homozygous product.

[0079] (4) Construction of composite score: A composite ranking score is constructed using homozygous product (HP) as the primary ranking key and genotypic diversity (GD) as the secondary ranking key. Specifically, the HP value is rounded to 6 decimal places as the high-order digit, and the decimal part (6 digits) of the GD value is appended to it as the low-order digit, forming a numerical composite score. This score ensures that the scores are first ranked in descending order of HP (prioritizing sites with the most homozygous differentiation), and when HP values ​​are the same, they are ranked in descending order of GD (prioritizing sites with higher genotypic diversity).

[0080] (5) Chromosome round-robin sorting: Within the same locus group with the same composite score, a chromosome round-robin mechanism is used for sorting. The specific method is as follows: If the scoring levels are the same, the candidate loci are first grouped by the chromosome where they are located, and the maximum value among the number of candidate loci on each chromosome is taken as the number of cycles. In each round, one locus is sequentially taken from each chromosome that still has remaining candidate loci to form a candidate queue that preferentially alternates across chromosomes; meanwhile, at the chromosome polling sorting level, a single-chromosome selection upper limit is set in combination with the current number of selected loci of each chromosome. Before a candidate locus enters the discriminative power increment evaluation, it is first judged whether the chromosome where it is located has reached the current upper limit. If the upper limit is reached, this candidate locus is skipped first in the current upper limit round; after the current upper limit round completes scanning the entire candidate locus list, the single-chromosome selection upper limit is increased, and the candidate locus list is rescanned from high score to low score, so that the candidate loci previously skipped under this scoring level get another chance to enter the discriminative power increment evaluation.

[0081] Step 3: Discriminative power increment evaluation and locus filling with adaptive weight decrement. A variety genotype fingerprint dictionary is established, and in the initial state, the fingerprint of each variety is an empty string. Each candidate locus is sequentially evaluated in the order sorted in Step 2. For each candidate locus, its genotype information is appended to the fingerprint string of each variety, and the number of varieties that can be uniquely distinguished after updating the variety genotype fingerprint dictionary (i.e., the number of varieties with mutually different fingerprint strings) is calculated.

[0082] An adaptive weight decrement strategy is adopted to evaluate the inclusion value of this locus, and the specific evaluation rules are as follows: Let the current number of distinguished varieties be D, the total number of varieties be N, the number of distinguishable varieties after addition be D_new, the current weight coefficient be W (initially 0.7), the dynamic threshold M = max(1, W × N), the stage dividing line S1 = 0.1 × N, and S2 = 0.9 × N + 0.1 × W × N. S1 refers to the distinction dividing line between the initial stage and the middle stage, and S2 refers to the distinction dividing line between the middle stage and the later stage. 0.1 in the above stage dividing line is an empirical value, under which fewer SNPs are contained. 0.9 is also an empirical value, under which more SNP loci that can improve the discrimination degree can be included.

[0083] Initial stage (D<S1): When D is 0, the locus is required to generate at least 2 different fingerprints (D_new - D ≥ 2) to be included; when D>0, any locus that can add at least one distinguishable variety is included; Middle stage (S1 ≤ D ≤ S2): The locus is included only when the number of newly added distinguishable varieties exceeds the dynamic threshold M (D_new - D>M), and this threshold is dynamically adjusted along with the weight coefficient W; Later stage (D>S2): Any locus that can add at least one distinguishable variety is included.

[0084] If all candidate sites are not distinguished under the current weighting coefficient (W) (and each site has been backtracked once), the weighting coefficient is reduced (from 0.7 to 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.01 until 0), and the candidate sites are re-traversed using the new weighting coefficient. As the weighting coefficient decreases, the dynamic threshold M decreases, and the starting boundary S2 in the later stage moves forward, gradually relaxing the screening criteria. The core idea of ​​this strategy is to quickly accumulate highly distinguishing sites with strict criteria in the early stage, and then gradually relax the inclusion criteria in the later stage to supplement sites that are effective for a few difficult-to-distinguish varieties, thereby achieving a balance between minimizing the number of sites and maximizing distinguishing ability.

[0085] If a candidate locus fails the discrimination assessment, the added genotype information will be removed from the fingerprint strings of all varieties (removal of a single locus) to ensure that the evaluation benchmark of subsequent loci is not affected.

[0086] Step 4: Removal of multi-copy or highly homologous regions. Sequence uniqueness verification of candidate sites: For each candidate site assessed for discriminative ability in Step 3, extract its flanking sequences (e.g., 150 bases upstream and downstream) on the reference genome to construct a local sequence fragment. Perform sequence similarity alignment (BLAST alignment) on this sequence against the whole genome to check the number of alignments that match the sequence in the genome.

[0087] If a sequence has multiple highly similar alignment sites in the genome (specific criteria: the similarity of the second alignment is greater than or equal to a preset threshold (e.g., 96%) and the alignment length is greater than or equal to a preset proportion of the total length of the flanking sequences (e.g., 96%)), then the site is determined to be located in a multi-copy or highly homologous region and is unsuitable as a specific molecular marker. It is excluded, and the next candidate site is evaluated. Only sites with a unique alignment location in the genome are retained as the final core SNP.

[0088] Step 5: Multiple rounds of execution to increase stability. To improve the robustness of the screening results, multiple rounds of iteration can be performed on the screening process of Steps 3-4 (e.g., 10 rounds). The screening results of each round are automatically added to the "known sites" and "skipped sites" lists of the next round, allowing subsequent rounds to expand from different candidate sites and accumulate a richer set of core SNP sites.

[0089] Example 2 This embodiment uses the SNP site screening method provided in Example 1 to screen 416 varieties of a certain Brassica species for SNPs.

[0090] Experimental materials: Whole genome resequencing data of 416 varieties of a certain Brassica species were used. The initially obtained SNP variant sites were filtered by biallelic filtering and used as input.

[0091] Adaptive filtering results: The system automatically determines the filtering threshold as follows: site quality value greater than 15375, deletion rate equal to 0 (no deletion), and average sequencing depth between 4 and 15 times.

[0092] Core SNP screening results: After one round of iterative screening, 29 core SNP loci were selected from a massive pool of candidate loci, which can uniquely distinguish all 416 varieties. Each locus passed BLAST genome-wide uniqueness verification (excluding non-chromosomal regions marked with "contig" generated during chromosome splicing). To increase stability, the process was repeated 10 times, with each round selecting between 29 and 42 loci, and each round could distinguish all 416 varieties.

[0093] Figure 2 A density map was selected for the core marker screening method for varieties, with C1-C9 representing chromosome number and the x-axis representing sequence length. The results showed that the selected core SNPs were evenly distributed across 9 chromosomes. The genotype QR codes selected by the core marker screening method for varieties are referenced below. Figure 3 As shown.

[0094] Example 3 The SNP screening method provided by this invention screens out fewer core SNPs, achieving comparable or even better variety differentiation.

[0095] Specifically, this embodiment uses the SNP site screening method provided in Embodiment 1 to screen core SNPs for peach (Prunus persica) varieties, and sets up a control experiment with the identification of unknown sample varieties as the application scenario to verify that the method of the present invention can achieve comparable or even better variety differentiation with fewer sites compared to existing sampling strategies.

[0096] 1. Experimental materials and application scenarios: This embodiment uses whole-genome resequencing data of peaches, including 337 samples. Of these, 335 samples have known varietal origins, forming a known varietal comparison library; the remaining two samples have unknown varietal origins and require identification, designated as Unknown Sample 1 and Unknown Sample 2, respectively. The objective is to determine the varietal affiliation of the unknown samples by calculating the genetic similarity at molecular marker sites between the two unknown samples and samples in the known varietal library, without relying on any known varietal tags. For ease of result presentation, in Appendix Table 1 below, known varietal samples are uniformly labeled as "breed" with a number, and the two unknown samples are labeled as Unknown1 and Unknown2, respectively.

[0097] 2. Construction of the control site set.

[0098] Based on the same peach SNP dataset (3,056,332 SNP loci) filtered by biallelic genes, three strategies were used to construct loci sets for genetic similarity calculation: (1) All loci group (benchmark group): All 3,056,332 SNP loci after filtering were used. Because it contains the most complete locus information, the identification results of this group are used as the benchmark (gold standard) for determining the true variety of unknown samples.

[0099] (2) Core site group of the present invention: Using the method described in Example 1 of this application, namely, sorting based on homozygous product-genotype diversity composite score, chromosome uniform polling, incremental evaluation of discriminative ability with adaptive weight reduction and sequence uniqueness verification, 28 core SNP sites were obtained from all sites after one round of screening.

[0100] (3) Systematic sampling site group (existing protocol control group): The conventional systematic sampling strategy of "representing the whole with the part" in the art is adopted, that is, sites are sampled at fixed intervals at equal distances on all 3,056,332 sites (the sampling interval step size is equal to the total number of sites divided by the number of target sites), resulting in 31 SNP sites evenly distributed on the genome. The reason for setting up this control group is that, in the conventional practice of representing whole genome information with a small number of markers, fixed interval systematic sampling is the most representative traditional scheme, and using it as a control can objectively reflect the advantages of the method of the present invention.

[0101] 3. Experimental Methods For the three locus sets mentioned above, the genotypes of all 337 samples at the corresponding loci were extracted. The genetic similarity between the two unknown samples (Unknown1, Unknown2) and each sample in the known variety library was calculated (i.e., the proportion of genotypes that are identical at common comparable loci between the two samples), and the samples were sorted from high to low similarity. The known sample with the highest similarity indicates the most likely variety affiliation of the unknown sample.

[0102] 4. Experimental Results (1) All loci groups (baseline group, 3,056,332 loci): Unknown1 had the highest similarity to breed849 (0.9652), significantly higher than breed796 (0.9196), which ranked second, with a difference of 0.0457, indicating a clear and unique attribution; Unknown2 had the highest similarity to breed847 (0.9655), significantly higher than breed898 (0.9277), which ranked second, with a difference of 0.0379, indicating a clear and unique attribution. Based on this, the true species attribution of the two unknown samples was determined as follows: Unknown1 belongs to breed849, and Unknown2 belongs to breed847.

[0103] (2) The core locus group of this invention (28 loci): Unknown1 has the highest similarity to breed849 (0.9259), which is 0.1111 different from breed860 (0.8148), which ranks second; Unknown2 has the highest similarity to breed847 (0.8929), which is 0.0714 different from breed846 (0.8214), which ranks second. The identification conclusion of this group is completely consistent with the benchmark group (Unknown1 belongs to breed849, and Unknown2 belongs to breed847), and the similarity difference between the most similar variety and the second most similar variety (0.1111, 0.0714) is significantly greater than that of the benchmark group (0.0457, 0.0379), indicating that this invention can give a more clear and unambiguous correct identification conclusion consistent with the full number of loci using only 28 loci.

[0104] (3) Systematic sampling site group (31 sites): The highest similarity of Unknown1 was tied (breed809 and breed893 were both 0.9667), making it impossible to uniquely determine its species classification, while the true species classification breed849 was only ranked 11th (0.9333); the highest similarity sample of Unknown2 was breed867 (0.9655), while the true species classification breed847 was pushed to 63rd (0.8710). That is, even though the number of sites (31) was more than that of the present invention (28), the systematic sampling group gave incorrect species classifications for both unknown samples, and there were ties and difficult interpretations.

[0105] For a clear comparison, the results of the three sets are summarized in Table 1 below: Table 1 Comparison of the species classification results of the three sets of loci for unknown samples

[0106] 5. Experimental Conclusions (1) Based on the identification results of the most complete set of all loci (3,056,332 loci), the method of the present invention can accurately reproduce the identification conclusion of the full set of loci with only 28 core SNP loci, and correctly assign the two unknown samples to breed849 and breed847 respectively. The number of loci used is reduced compared to the full set of loci (from 3,056,332 to 28).

[0107] (2) When the number of loci is comparable to that of the present invention (31 loci), the existing conventional fixed interval sampling strategy cannot give the correct variety classification. There are cases where two unknown samples are incorrectly classified and difficult to interpret. This shows that simply relying on a small number of loci without optimizing the ability to distinguish between them is not enough to reliably identify varieties.

[0108] (3) The similarity difference between the most similar variety and the second most similar variety in the core locus group of the present invention (0.1111, 0.0714) is greater than that in the full locus group (0.0457, 0.0379), indicating that the loci selected by the present invention are targeted at the differences between varieties, and the clarity and reliability of the identification conclusions are improved while reducing the number of loci.

[0109] In summary, this embodiment demonstrates that the SNP screening method provided in this application can achieve variety differentiation comparable to that of the full number of loci, with a significantly fewer loci.

[0110] Example 4 This embodiment provides a computer device, including a memory, a processor, a computer program stored on the memory, a computer-readable storage medium, and a computer program product. When the processor executes the computer program, it implements the above-described method for screening variety SNP markers.

[0111] Figure 4 An example is a schematic diagram of the physical structure of an electronic device, such as... Figure 4 As shown, the electronic device may include a processor 22, a memory 21, a communication unit 23, and a communication bus 24, wherein the processor 22, the memory 21, and the communication unit 23 communicate with each other through the communication bus 24. The processor 22 can call logical instructions in the memory 21 to execute a method for screening variety SNP tags.

[0112] Furthermore, when the logical instructions in the aforementioned memory 21 can be implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium, including: USB flash drive, portable hard drive, read-only memory (ROM), random access memory (RAM), magnetic disk or optical disk, and other media that can store program code.

[0113] On the other hand, the present invention also provides a computer program product, which includes a computer program that, when executed by a processor, enables the computer to perform the above-described method for screening SNP markers of the execution variety.

[0114] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for screening SNP markers for a variety, characterized in that, It includes the following steps: S1: Perform adaptive filtering on variant sites: Manually or automatically calculate the following quality metrics of the variant sites: site quality value, missing rate and sequencing depth, determine the filtering threshold for each metric, filter the variant sites according to the specified quality metrics, and retain the sites that pass the thresholds for all specified metrics as a candidate site set; S2: Sort the candidate site set: Sort the candidate site set by using a genetic diversity indicator; the genetic diversity indicator is at least one selected from the group consisting of homozygote product, genotype diversity index, nucleotide polymorphism Pi, expected heterozygosity He, polymorphism information content PIC and Shannon's information index; S3: Chromosome polling sorting: Within the candidate site set with the same genetic diversity indicator, a chromosome polling mechanism is adopted for arrangement, so that the candidate sites are preferably from different chromosomes alternately, and the candidate sites enter an alternately arranged candidate queue for sorting; when the available sites on some chromosomes are exhausted in advance, "polling skip" is performed for the chromosomes exhausted in advance; and adjacent sites among the remaining available sites from the same chromosome are taken as candidate sites to enter the alternately arranged candidate queue for sorting; S4: Distinguishing ability increment evaluation with adaptive weight decrement and site filling: establishing a variety genotype fingerprint dictionary, where the fingerprint of each variety is an empty string in the initial state, evaluating each candidate site in sequence according to the sorted order obtained in step S3, for each candidate site, appending its genotype information to the fingerprint string of each variety, and calculating the number of varieties that can be uniquely distinguished after updating the variety genotype fingerprint dictionary; According to the site inclusion rule, sites that can increase the number of distinguishable varieties are included in the variety genotype fingerprint dictionary; If all varieties are not yet distinguished after traversing all candidate sites under the current weight coefficient, reduce the weight coefficient and re-traverse the candidate sites using the new weight coefficient; S5: Removing multi-copy or highly homologous regions.

2. The method for screening variety SNP markers according to claim 1, characterized in that, The site inclusion rule in step S4 includes: Let D be the number of varieties already distinguished in the current candidate site set, N be the total number of varieties, D_new be the number of varieties distinguishable after addition, W be the current weight coefficient, the dynamic threshold M = max(1, W × N), and the phase dividing lines S1 = 0.1×N, S2 = 0.9×N + 0.1×W×N; S1 refers to the distinguishability dividing line between the initial phase and the middle phase, and S2 refers to the distinguishability dividing line between the middle phase and the final phase; In the initial phase, where D<S1: when D is 0, if D_new - D ≥ 2, this site can be included in the variety genotype fingerprint dictionary; when D>0, any site that can increase at least one distinguishable variety is included in the variety genotype fingerprint dictionary; In the middle phase, where S1 ≤ D ≤ S2: if D_new - D>M, this site can be included in the variety genotype fingerprint dictionary; In the final phase, where D>S2: any site that can increase at least one distinguishable variety is included in the variety genotype fingerprint dictionary.

3. The method for screening variety SNP markers according to claim 2, characterized in that, W is initially set to 0.

7. If all candidate sites are still not distinguished under the current weight coefficient, the weight coefficient is reduced and a new weight coefficient is used to re-traverse the candidate sites. In the later stage, the inclusion criteria are gradually relaxed to supplement sites that are effective for a few difficult-to-distinguish varieties. Preferably, the weighting coefficient W gradually decreases from 0.7 to 0.

4. The method for screening variety SNP markers according to claim 2, characterized in that, If none of the candidate loci in the candidate locus set pass the incremental discrimination assessment, the added genotype information will be removed from the genotype fingerprint dictionary of all varieties.

5. The method for screening variety SNP markers according to claim 1, characterized in that, Step S5, removal of multiple copies or highly homologous regions, includes: extracting flanking sequences from the reference genome for each candidate site that has passed the discrimination ability assessment, and constructing local sequence fragments; comparing the local sequence fragments with the whole genome for sequence similarity, and removing candidate sites with multiple highly similar alignment positions; Preferably, the criterion for determining the existence of multiple highly similar alignment positions is as follows: if the alignment similarity of the second similar sequence is greater than or equal to a preset threshold and the alignment length is greater than or equal to a preset proportion of the total length of the flanking sequences, then it is determined that there are multiple highly similar alignment positions.

6. The method for screening variety SNP markers according to claim 1, characterized in that, After step S5 is completed, steps S4 and S5 are executed in multiple rounds. The results of each round of screening are automatically added to the "known sites" and "skipped sites" lists of the next round, so that subsequent rounds can expand from different candidate sites to obtain the core SNP site set.

7. The method for screening variety SNP markers according to claim 1, characterized in that, In step S2, the genetic diversity index is selected from the homozygous product and the genotype diversity index, and the two are combined and scored; the composite score is constructed with the homozygous product as the primary sorting key and the genotype diversity index as the secondary sorting key. The formula for calculating the homozygous product is: HP = f_ref × f_alt; f_ref refers to the frequency of the homozygous reference genotype in the sample, and f_alt refers to the frequency of the homozygous variant genotype in the sample; The formula for calculating the genotype diversity index GD is: GD = 1 - Σ(fi²); where fi is the frequency of the i-th genotype; The scoring method for the composite score includes: retaining 3-8 decimal places for the HP value as the high digit, and concatenating 4-8 decimal places of the GD value as the low digit to form a numerical composite score. First, the scores are sorted in descending order of HP, and if HP is the same, they are sorted in descending order of GD.

8. The method for screening variety SNP markers according to claim 7, characterized in that, In step S3, the chromosome polling and sorting method includes: if the score levels are the same, the candidate sites are first grouped according to their respective chromosomes, and the maximum number of candidate sites in each chromosome is used as the number of iterations. In each round, one site is selected from each chromosome that still has remaining candidate sites, forming a candidate queue that is preferentially distributed across chromosomes. At the same time, at the chromosome polling and sorting level, a single chromosome selection upper limit is set based on the number of currently selected sites on each chromosome. Before a candidate site enters the incremental evaluation of discrimination ability, it is first determined whether its chromosome has reached the current upper limit. If it has, the candidate site is skipped in the current upper limit round. After the current upper limit round has scanned the entire candidate site list, the single chromosome selection upper limit is increased, and the candidate site list is scanned again from high score to low score, so that the candidate sites that were previously skipped at this score level have another chance to enter the incremental evaluation of discrimination ability. Preferably, in step S1, the filtering threshold for the deletion rate is automatically determined by using the interquartile range statistical method; the filtering threshold for the sequencing depth is determined by the average value of the sequencing depth; and the filtering threshold for the site quality value is determined by taking half of the average quality value of all sites. Preferably, the average depth threshold is taken as 0.5 to 2 times the mean; the missing rate threshold is taken as the upper bound of IQR.

9. A computer device, comprising a memory, a processor, and a computer program stored in the memory, characterized in that, When the processor executes the computer program, it implements the method for screening variety SNP markers as described in any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the method for screening variety SNP markers as described in any one of claims 1 to 8.

Images