Method for developing polymorphic molecular markers of leymus chinensis based on genome sequencing

Abstract

The application discloses a method for developing polymorphic molecular markers of Leymus mollis based on genome sequencing, and belongs to the field of genomics and molecular marker technology. The method comprises the following steps: obtaining genome sequencing raw data of a plurality of samples of Leymus mollis, extracting candidate regions with inconsistent coverage depth from the alignment results after quality control and alignment, and establishing a variation confidence score matrix based on base substitution frequency and insertion and deletion length distribution; performing collinearity clustering on all sample variation sites according to the score matrix to generate variation site clusters consistent across samples; performing local multiple alignment on the variation site clusters to identify simple repeat sequence regions and single nucleotide site regions; designing flanking primer pairs and competitive allele-specific probes respectively to form a primary molecular marker set; and using raw data to simulate amplification verification, screening markers with amplification success rate and polymorphism information content both exceeding a threshold value, and obtaining polymorphic molecular markers of Leymus mollis.

Description

Technical Field

[0001] This invention relates to the fields of genomics and molecular marker technology, specifically to a method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing. Background Technology

[0002] Leymus chinensis is an important forage and ecological grass species, and the development of polymorphic molecular markers is of great value for its genetic diversity assessment and molecular-assisted breeding. Techniques for developing molecular markers based on genome sequencing typically rely on detecting variations in single or small samples, then directly converting candidate sites into markers for experimental validation. When processing population sequencing data from multiple samples, this approach has significant shortcomings in terms of variation site accuracy and cross-sample consistency. On the one hand, relying solely on single features such as base substitutions or insertions / deletions to identify variations easily leads to misclassification of sequencing errors and alignment biases as true variations, resulting in a large number of false positive sites entering the downstream marker design stage. On the other hand, in a multi-sample context, the lack of testing for collinearity among variation sites in the population means that the output site set often contains a large number of "private" variations that exist only in individual samples or in regions of alignment errors. Markers designed from these sites exhibit unstable polymorphism in the population, resulting in a high amplification and validation failure rate. Furthermore, existing workflows often separate the detection of simple repetitive sequences and single nucleotide polymorphism sites, making it difficult to simultaneously and systematically obtain both types of markers from the same set of multi-sample sequencing data, leading to low marker conversion efficiency. The urgent problem to be solved is how to accurately identify cross-sample consistent polymorphic sites from population genome sequencing data and efficiently convert them into molecular markers with high polymorphism and high amplification success rate. Summary of the Invention

[0003] This paper provides a method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing, which solves the problems of low accuracy in identifying variant sites in multi-sample genome sequencing data, poor cross-sample consistency, and low efficiency in molecular marker transformation and verification in existing technologies. This method enables efficient and accurate development of highly polymorphic molecular markers from population sequencing data.

[0004] The objective of this invention can be achieved through the following technical solutions: This invention provides a method for developing polymorphic molecular markers in Leymus chinensis based on genome sequencing. By systematically identifying, clustering, and screening variants in the raw genome sequencing data of multiple Leymus chinensis samples, high-throughput and high-reliability polymorphic molecular marker development is achieved. The method first acquires raw genome sequencing data from multiple Leymus chinensis samples, performs quality filtering and error correction on each read, discarding low-quality reads and trimming low-quality bases at the ends. Preferably, during error correction, reads are split into k-mers, the frequency of each k-mer is counted, low-frequency k-mers are marked as erroneous k-mers, and reads containing erroneous k-mers are attempted to be corrected by replacing bases with high-quality ones. Correction is only accepted if the frequency of the corrected k-mers meets the requirements; otherwise, the read is discarded, thereby effectively eliminating the interference of sequencing errors on variant detection.

[0005] After obtaining high-quality sequencing data, the reads of each sample were initially aligned with the *Leymus chinensis* reference genome. The distribution of paired end distances was statistically analyzed, and reads with abnormal distances were removed. The remaining reads were then used to perform local re-alignment on regions that did not align or had low alignment quality in the initial alignment, generating a final alignment record sorted by coordinates, which served as the reference sequence alignment result for each sample. Based on this, the alignment results were scanned using a sliding window approach, and the average coverage depth and standard deviation of each window were calculated. Windows whose average coverage depth deviated from the average depth of the whole genome by more than a preset multiple were marked as depth aberration windows, and base substitution events and insertion / deletion events within them were recorded. Adjacent depth aberration windows were merged into continuous intervals. When both the interval length and base substitution frequency met the criteria, they were identified as candidate regions with inconsistent coverage depth. This strategy, based on a joint screening of depth deviation and variation events, can accurately identify genomic regions with genuine sequence variations and eliminate false positive regions caused by alignment errors or random fluctuations.

[0006] For each candidate region, the frequency of non-reference bases at each base position in all samples is counted as the base substitution frequency. The sequence length and occurrence frequency of each insertion / deletion event are recorded, and an insertion / deletion length distribution histogram is constructed. The initial substitution confidence score is obtained by multiplying the base substitution frequency by the inverse of the sequencing base quality value. The initial insertion / deletion confidence score is obtained by multiplying the mode length of the insertion / deletion length distribution by the occurrence frequency of the event. These are then arranged according to sample and location dimensions to form a variation confidence score matrix. Preferably, each element in the matrix is ​​a weighted sum of the initial substitution confidence score and the initial insertion / deletion confidence score, thus comprehensively reflecting the reliability of variations in different samples at different sites.

[0007] The variant sites shared by all samples are extracted from the variant confidence score matrix. Each site is represented as a vector with the confidence scores of each sample as elements, and the cosine similarity between any two site vectors is calculated. Variants with a cosine similarity greater than a similarity threshold are grouped into the same initial cluster. The similarity threshold can be set to 0.85, 0.90, or 0.95. For each initial cluster, it is checked whether the genomic coordinate order of the sites within the cluster remains consistent across all samples. If it is consistent, the initial cluster is retained; otherwise, it is split into multiple sub-clusters based on the degree of coordinate order conflict, ultimately resulting in a cross-sample consistent variant site cluster. Using vector similarity to construct collinear clusters can effectively aggregate variant sites that exhibit consistent linkage behavior across different samples, improving the generality and representativeness of the markers.

[0008] For each cluster of variant sites, a predetermined flanking length is extended outward from its boundary to extract the original sequencing reads from each sample within this interval. These reads are then input into a progressive multiplexer to generate a local multiplex alignment matrix. Each column of the alignment matrix is ​​scanned. If a column contains three or more different base types, and the frequency of each type exceeds the lower frequency limit, the column is identified as a single nucleotide polymorphism (SNP) site. Consecutive SNP sites are then merged into a sequence polymorphism (SPP) region. Simultaneously, multiple consecutive columns are scanned. If an insertion / deletion pattern appears cyclically with a fixed repeat unit length, and the number of repeats varies across different samples, the interval is identified as a simple repeat sequence region. This step simultaneously identifies both length polymorphism and sequence polymorphism from the variant site clusters, providing clear target regions for subsequent marker design.

[0009] After obtaining the two types of sequence regions, corresponding molecular markers were designed. For simple repetitive sequence regions, conserved sequences without repetitive units on both sides were extracted as primer design templates. Flanking primer pairs were generated according to the constraints of consistent primer annealing temperature, no hairpin structure formation, and amplification product length within the target interval. For single nucleotide site regions, upstream and downstream flanking sequences were extracted. Two competitive allele-specific probes were designed based on the base types of the two alleles. The two probes were of the same length and the difference in annealing temperature did not exceed the temperature deviation threshold. Locked nucleic acid (LNI) modification was introduced on the terminal base of the probe near the single nucleotide site; preferably, the LNI modification was located at the penultimate or antepenultimate position of the 3′ end of the probe. All flanking primer pairs and their corresponding simple repetitive sequence regions, single nucleotide site regions and their corresponding competitive allele-specific probes were associated to form a primary molecular marker set. By differentially designing detection primers and probes of the two marker types and introducing LNI modification to enhance single-base discrimination ability, the accuracy and sensitivity of subsequent genotyping were significantly improved.

[0010] Before performing simulated amplification and verification of the primary molecular marker set using all original sample data, as a preferred technical solution of the present invention, the primary molecular marker set is redundantly merged. That is, the physical distance between any two markers on the reference genome is calculated, and marker pairs with a distance less than a threshold are marked as neighboring marker pairs. The primer or probe sequences of neighboring marker pairs are compared. If the sequence similarity exceeds the threshold, they are determined to be redundant marker pairs, and markers with low polymorphism information content are removed from the set. After sorting by genomic location, a redundancy-free primary molecular marker set is generated to avoid introducing redundant detection information due to markers covering the same or highly similar regions.

[0011] Finally, for each marker in the primary molecular marker set, reads from the corresponding region were extracted from the sequencing data of all samples. The amplification process was simulated, and the proportion of samples that successfully amplified products of the expected length was counted as the amplification success rate. Genotypes were calculated based on the abundance of different alleles in the amplified products, and the frequencies of each allele were statistically analyzed. These values ​​were then substituted into the polymorphism information content (PI) calculation formula to obtain the PI value for that marker. Markers with an amplification success rate greater than a success rate threshold and a PI value greater than an information content threshold were output as the final Leymus chinensis polymorphic molecular markers. This simulation validation strategy, without relying on additional wet experiments, utilizes existing sequencing data to perform dual screening of marker availability and polymorphism levels, ensuring that the developed molecular markers have good amplification stability and efficient polymorphism information capture capabilities in the population.

[0012] The beneficial effects of this invention are: By extracting candidate regions with inconsistent coverage depth and establishing a variant confidence scoring matrix, this approach changes the traditional reliance on a single variant signal for mutation detection. The confidence levels of substitution and insertion / deletion events are quantified by multiplying the base substitution frequency by the inverse of the sequencing base quality and the insertion / deletion length distribution mode by the number of occurrences, respectively. Scoring matrices are then constructed along sample and location dimensions. This multi-evidence fusion confidence scoring method can distinguish between true allelic variants and signals induced by random sequencing errors or low-quality alignments, significantly reducing the proportion of false positives entering downstream analysis.

[0013] Based on the variant confidence score matrix, collinear clustering is performed on variant sites in all samples, representing each variant site as a cross-sample confidence vector. Cosine similarity is used to calculate the overall variant pattern similarity between sites. Simultaneously, a genome coordinate order consistency test is used to retain or split the initial clusters, generating consistent variant site clusters across samples. This process ensures that the final output variant site clusters possess stable collinearity characteristics within the population, eliminating scattered false variants that only appear in a few samples or in local alignment anomalies. The obtained site clusters more accurately reflect the common polymorphic regions within the population.

[0014] After identifying simple repetitive sequence regions and single nucleotide site regions, differentiated marker design strategies were employed. For simple repetitive sequence regions, flanking primer pairs with consistent annealing temperatures, no hairpin structures, and controlled product lengths were designed. For single nucleotide site regions, competitive allele-specific probes modified with locked nucleic acids were designed, with both probes having the same length and similar annealing temperatures. This genotyping design effectively optimized the amplification and genotyping efficiency of both types of markers, with locked nucleic acid modification particularly enhancing the probes' ability to distinguish single-base mismatches.

[0015] Based on this, the primary molecular marker set was validated by simulated amplification using raw sequencing data from all samples. This involved extracting reads from target regions to simulate primer or probe binding and extension processes, statistically analyzing amplification success rates, and calculating polymorphism information content based on allele read abundance. Only markers with both indicators exceeding thresholds were retained. This validation step replaced batch wet screening with real-world population sequencing data, reducing experimental trial-and-error costs while ensuring that the output markers possess both high amplification success rates and rich polymorphic information. Attached Figure Description

[0016] The invention will now be further described with reference to the accompanying drawings.

[0017] Figure 1 This is a flowchart of the method for developing polymorphic molecular markers in Leymus chinensis based on genome sequencing; Figure 2 This is a flowchart of sequencing data quality filtering, error correction, and reference sequence alignment. Figure 3 This is a flowchart of the process for extracting candidate regions with inconsistent coverage depth and constructing the mutation confidence score matrix; Figure 4 This is a flowchart of the process for generating clusters of consistent variant sites across samples; Figure 5 This is a flowchart of molecular marker design based on local multiple alignment; Figure 6 This is a flowchart of the redundancy merging and simulated amplification verification process for Leymus chinensis polymorphic molecular markers. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] See Figure 1This invention provides a method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing, comprising the following steps: acquiring raw genome sequencing data from multiple Leymus chinensis samples; performing quality filtering and error correction on the raw data of each sample to construct a reference sequence alignment result for each sample; extracting candidate regions with inconsistent coverage depth from the alignment results of each sample, and establishing a variant confidence score matrix based on the base substitution frequency and insertion / deletion length distribution within the candidate regions; performing collinear clustering of variant sites in all samples according to the variant confidence score matrix to generate cross-sample consistent variant site clusters; performing local multiple alignment on the sequences within each variant site cluster to identify simple repetitive sequence regions with length polymorphism and single nucleotide site regions with sequence polymorphism; designing flanking primer pairs for simple repetitive sequence regions and competitive allele-specific probes for single nucleotide site regions to form a primary molecular marker set; and performing simulated amplification and verification of the primary molecular marker set using the raw genome sequencing data of all samples, and selecting markers with amplification success rate and polymorphism information content exceeding a threshold as Leymus chinensis polymorphic molecular markers.

[0020] In specific implementation, please refer to Figure 2 The process of quality filtering and error correction for the raw data of each sample is as follows: Genome sequencing data from multiple samples of *Leymus chinensis* were obtained. Each sample's genome sequencing data contained multiple reads, each composed of multiple bases, with each base corresponding to a sequencing quality value. For each read in the genome sequencing data of each sample, the sequencing quality value was read base by base, and the arithmetic mean of the sequencing quality values ​​of all bases in the read was calculated. This arithmetic mean was used as the average quality value of the read. Reads with an average quality value lower than the average quality threshold were removed. The average quality threshold was set to Q20, meaning reads with an average quality value lower than Q20 were discarded. Simultaneously, the remaining reads underwent end trimming. The sequencing quality value was checked base by base from the 3′ end to the 5′ end of the read. When consecutive sequencing quality values ​​were lower than the end quality threshold, these consecutive bases were trimmed. The end quality threshold was set to Q15. After end trimming, the remaining portion of the read was used as the filtered read.

[0021] The filtered reads are split into all possible k-mers according to a preset k-mer length of 17 bases. Using a sliding window of 17 bases, starting from the 5' end of the filtered read, a k-mer of 17 bases is extracted for each base moved, until the window reaches the 3' end of the filtered read, thus obtaining all k-mers contained in the filtered read. The frequency of each k-mer in all filtered reads of the same sample is counted. K-mers with frequencies below a low frequency threshold of 2 are marked as erroneous k-mers. For reads containing erroneous k-mers, the base position corresponding to the erroneous k-mer is identified, and an attempt is made to replace the base at the corresponding position with the highest quality substitute base. The method for determining the substitution base with the highest quality value is as follows: Given an erroneous k-mer occurrence frequency of 1, each base in the erroneous k-mer is replaced with one of the other three base types, generating a corrected k-mer. The frequency of the corrected k-mer is then checked across all filtered reads of the same sample. The substitution base that maximizes the frequency of the corrected k-mer is selected as the substitution base with the highest quality value. A new corrected k-mer is generated. If the frequency of the corrected k-mer exceeds a high frequency threshold (set to 3), the correction is accepted; otherwise, the read containing the erroneous k-mer is discarded. After this quality filtering and error correction process, each sample yields a set of corrected sequencing reads.

[0022] The process of constructing the reference sequence alignment results for each sample is as follows: The quality-filtered and error-corrected sequencing reads of each sample are initially aligned with the *Leymus chinensis* reference genome using the BWA-MEM alignment tool and default settings. This yields the initial alignment position for each quality-filtered and error-corrected sequencing read. Based on the initial alignment position, the distribution of paired-end distances for each quality-filtered and error-corrected sequencing read is statistically analyzed. Paired-end distance refers to the fragment length between the two starting positions on the *Leymus chinensis* reference genome where a pair of reads aligns. The arithmetic mean and standard deviation of all paired-end distances are calculated. Abnormal reads with paired-end distances exceeding three standard deviations are removed; that is, read pairs whose absolute value of the difference between the paired-end distance and the arithmetic mean is greater than three times the standard deviation are discarded.

[0023] Remaining reads are used to perform local re-alignment on regions that did not align or whose alignment quality was below the quality threshold in the initial alignment. Remaining reads refer to the quality-filtered and error-corrected sequencing reads retained after removing abnormal reads. Regions with alignment quality below the quality threshold refer to the genomic regions corresponding to alignment records with alignment quality values ​​below 20. Local re-alignment uses GATK's IndelRealigner tool, with the *Leymus chinensis* reference genome as a reference, to extract reads covering regions that did not align or whose alignment quality was below the quality threshold from the remaining reads. Multiple alignments and relocations are performed in these local regions to generate the final alignment record for each sample. The final alignment record includes the alignment coordinates, alignment quality value, CIGAR string, and sequence information of the read. The final alignment records are sorted by alignment coordinates, and the reference sequence alignment results for each sample are output in BAM file format.

[0024] In specific implementation, please refer to Figure 3 The process of extracting candidate regions with inconsistent coverage depth from the alignment results of each sample is as follows: The reference sequence alignment results for each sample are obtained and stored in BAM file format. These results include the alignment coordinates of all quality-filtered and error-corrected sequencing reads in that sample on the *Leymus chinensis* reference genome, CIGAR strings, and sequencing quality values. The reference sequence alignment results for each sample are then scanned using a preset window length and step size. The window length is set to 500 base pairs, and the step size is set to 100 base pairs. Starting from the initial coordinates of the *Leymus chinensis* reference genome, a window of 500 base pairs is truncated every 100 base pairs until the entire chromosome of the *Leymus chinensis* reference genome is covered. For each window, the samtoolsdepth command is called to extract the number of sequencing reads covering each base position within the window. The arithmetic mean of the number of sequencing reads covering all base positions within the window is calculated, and the arithmetic mean is used as the average coverage depth within the window. At the same time, the standard deviation of the number of sequencing reads covering all base positions within the window is calculated, and the standard deviation is used as the standard deviation of the coverage depth within the window.

[0025] The average coverage depth of the whole genome is calculated by summing the average coverage depths of all windows and dividing by the total number of windows. A fold change threshold of 2.0 is set; a window is marked as a depth aberration window when its average coverage depth is greater than the average coverage depth of the whole genome multiplied by 2.0, or when its average coverage depth is less than the average coverage depth of the whole genome divided by 2.0. The criteria for determining a depth aberration window are that copy number variations or large fragment insertions / deletions cause the coverage depth to deviate from the expected coverage level of the whole genome. The fold change threshold of 2.0 corresponds to the signal intensity of heterozygous deletions or single copy duplications in a diploid genome. Base substitution events and insertion / deletion events are recorded within the depth aberration window. Base substitution events refer to the presence of base types in the aligned reads within the depth aberration window that are inconsistent with the bases in the *Leymus chinensis* reference genome. Insertion / deletion events refer to the presence of I or D operations in the CIGAR string of the aligned reads within the depth aberration window.

[0026] Adjacent depth anomaly windows are merged into contiguous intervals. Adjacent means the difference in the starting coordinates of two depth anomaly windows is less than or equal to the window length of 500 base pairs. After merging, one or more contiguous intervals are obtained, each consisting of at least one depth anomaly window. A minimum length limit is set at 1000 base pairs. When the length of a contiguous interval exceeds 1000 base pairs, the base substitution frequency within that interval is further calculated. The base substitution frequency is calculated by counting the total number of bases involved in all base substitution events within the contiguous interval, dividing this total by the total number of covered bases within the interval. A frequency threshold is set at 0.05; if the base substitution frequency is higher than 0.05, the contiguous interval is identified as a candidate region with inconsistent coverage depth. The frequency threshold of 0.05 corresponds to the lower limit of the expected frequency of true variant sites in the genome; substitution signals below this frequency may originate from sequencing errors rather than biological variations.

[0027] The process of establishing a variation confidence scoring matrix based on the base substitution frequency and insertion / deletion length distribution within candidate regions is as follows: For each base position within each candidate region with inconsistent coverage depth, the number of times a non-reference base appears at that base position in all samples is counted. The reference base is the base type of the *Leymus chinensis* reference genome at that position, and the non-reference base is a base type different from the reference base. The non-reference base frequency is calculated as the base substitution frequency, which is equal to the number of times a non-reference base appears at that base position divided by the total number of times that base position is covered in all samples. For each insertion / deletion event within each candidate region with inconsistent coverage depth, the length of the inserted or deleted sequence is recorded. The sequence length is extracted from the CIGAR string, with positive values ​​for the sequence length corresponding to insertion operations and negative values ​​for the sequence length corresponding to deletion operations. The number of occurrences of the insertion / deletion event in all samples is counted, and an insertion / deletion length distribution histogram is constructed. The x-axis of the insertion / deletion length distribution histogram represents the length of the inserted or deleted sequence, and the y-axis represents the number of occurrences.

[0028] The initial substitution confidence score is obtained by multiplying the substitution frequency at each base position by the reciprocal of the corresponding sequencing base quality value. The corresponding sequencing base quality value refers to the arithmetic mean of the sequencing quality values ​​corresponding to non-reference bases in all sequencing reads where a non-reference base appears at that position. The reciprocal of the sequencing base quality value converts the sequencing quality value into a sequencing error probability; the higher the quality value, the smaller the reciprocal. Multiplying the substitution frequency by the sequencing error probability results in higher confidence for substitution signals of high-quality bases. The initial insertion / deletion confidence score is obtained by multiplying the mode length of the histogram of the insertion / deletion length distribution for each insertion / deletion event by the number of times the event occurs. The mode length of the histogram of the insertion / deletion length distribution is the sequence length value that appears most frequently in the histogram, and the number of times the event occurs refers to the total number of times the insertion / deletion event occurs in all samples.

[0029] The preliminary confidence scores for substitutions and insertions / deletions within all candidate regions with inconsistent coverage depths are arranged by sample and location dimensions to form a variant confidence score matrix. Rows in the variant confidence score matrix correspond to samples, and columns correspond to the base positions or genomic coordinate intervals where insertion / deletion events occur within candidate regions. The value of each element in the matrix is ​​calculated using the following formula:

[0030] in, In the confidence score matrix of the variance, the first... The sample at the th Confidence scores at each mutation location For sample index number, The value range is from 1 to the total number of samples. This is the index number of the mutation location, which includes the location of base substitution and the interval where the insertion / deletion event occurs. The value ranges from 1 to the total number of variant locations in all candidate regions. Indicates the first The sample at the th The initial confidence score of the permutation at the i-th mutation position, if the i-th If each mutation position is a base substitution position, then Take the initial confidence score for the substitution at that base position; if the first If the mutation positions are the intervals where insertion and deletion events occur, then... Take 0. Indicates the first The sample at the th The initial confidence score of insertion / deletion at the mutation location, if the ... If the mutation positions are the intervals where insertion and deletion events occur, then... Take the initial insertion-missing confidence score for this insertion-missing event; if the first... If each mutation position is a base substitution position, then Take 0. The confidence weighting coefficient for substitution is set to 0.6, based on the fact that the detection of base substitution events is less affected by sequencing errors, thus assigning a higher weight. The insertion missing confidence weight coefficient is set to 0.4, based on the fact that the detection of insertion missing events is greatly affected by alignment errors, and therefore a low weight is assigned.

[0031] In specific implementation, please refer to Figure 4 The process of generating consistent cross-sample variant site clusters by performing collinear clustering on variant sites of all samples based on the variant confidence score matrix is ​​as follows: Variants shared by all samples are extracted from the variant confidence score matrix. The variant confidence score matrix is ​​a two-dimensional numerical matrix. Each row of the variant confidence score matrix corresponds to a sample, and the number of rows is the total number of samples. Each column of the variant confidence score matrix corresponds to a variant position, and the number of columns is the total number of variant positions detected in all candidate regions. The types of variant positions include base substitution positions and intervals containing insertion / deletion events. The value in the u-th row and v-th column of the variant confidence score matrix represents the confidence score of the u-th sample at the v-th variant position. Variants shared by all samples refer to the variant sites in the variant confidence score matrix where, for a given variant position v, all samples have a confidence score greater than zero at that variant position v. After extracting the variant sites shared by all samples, a vector is constructed for each extracted variant site. The number of elements in the vector is equal to the total number of samples, and the value of the a-th element in the vector is the confidence score of the a-th sample at that extracted variant site.

[0032] Calculate the cosine similarity between any two extracted variant site vectors. The cosine similarity is calculated using the following formula:

[0033] in, Let represent the cosine similarity between the p-th extracted variant site vector and the q-th extracted variant site vector. and This refers to the index number of the extracted variant sites. and The values ​​range from 1 to the total number of extracted variant sites, and . This represents the total number of samples, i.e., the number of elements in each extracted mutation site vector. This represents the confidence score of the i-th sample in the p-th extracted variant site vector. For sample index number, The value range is 1 to . Let represent the confidence score of the i-th sample in the q-th extracted variant site vector. The numerator is the dot product of the two variant site vectors, and the denominator is the product of the magnitudes of the two variant site vectors.

[0034] Extracted variant sites with cosine similarity greater than a similarity threshold are grouped into the same initial cluster. The similarity threshold is set to one of 0.85, 0.90, or 0.95. The selection criteria are as follows: when the similarity threshold is set to 0.85, the angle between the corresponding cosine similarity values ​​is approximately 31.8 degrees, allowing for a moderate degree of confidence score variation; when the similarity threshold is set to 0.90, the angle between the vectors is approximately 25.8 degrees, requiring high consistency in vector direction; when the similarity threshold is set to 0.95, the angle between the vectors is approximately 18.2 degrees, clustering only variant sites with highly similar confidence score patterns into one cluster. The specific value of the similarity threshold can be selected from 0.85, 0.90, or 0.95 based on the degree of conservation of variation among the samples. The clustering process uses a connected component strategy: the extracted variant sites are used as nodes. When the cosine similarity between two extracted variant sites is greater than the selected similarity threshold, an undirected edge is established between the two extracted variant sites. After all undirected edges are established, all connected components in the graph are extracted by breadth-first search or depth-first search. The extracted variant sites corresponding to the nodes in each connected component constitute an initial cluster.

[0035] For each initial cluster, the order of genomic coordinates of the variant sites within the initial cluster is extracted. The order of genomic coordinates refers to the linear arrangement of all variant sites within the initial cluster, arranged in ascending order according to their chromosome number on the *Leymus chinensis* reference genome, and then arranged in ascending order of their initial coordinates within the same chromosome. For any two variant sites A and B within the initial cluster arranged by genomic coordinates, if the coordinate of variant site A on the *Leymus chinensis* reference genome is smaller than the coordinate of variant site B, then variant site A is placed before variant site B in the order. The order of any two variant sites within the initial cluster is checked one by one to ensure consistency across all samples. Consistency is determined by either all samples having a coordinate smaller than or greater than the coordinate of variant site B. If the order of coordinates for all pairs of variant sites within the initial cluster remains consistent across all samples, the initial cluster is retained.

[0036] If, within an initial cluster, the coordinate order of at least one pair of variant sites is inconsistent across different samples (i.e., in some samples, the coordinates of variant site A are less than the coordinates of variant site B, while in other samples, the coordinates of variant site A are greater than the coordinates of variant site B), then the initial cluster is determined to have a coordinate order conflict. The initial cluster with coordinate order conflicts is split into multiple sub-clusters according to the degree of conflict. The splitting process is as follows: a directed order graph of variant sites within the initial cluster is constructed, where each node represents a variant site within the initial cluster. For each pair of variant sites, the number of samples where the coordinates of variant site A are less than the coordinates of variant site B is counted. If this number exceeds half the total number of samples, a directed edge from variant site A to variant site B is added to the directed order graph. Directed cycles are detected in the directed order graph, with each directed cycle corresponding to a set of variant sites with coordinate order conflicts. For each directed cycle, the absolute value of the difference between the number of supporting samples and the number of opposing samples for each directed edge in the cycle is calculated, and the directed edge with the smallest absolute difference is removed to break the directed cycle. Repeat the process of detecting and breaking directed cycles until the directed order graph becomes a directed acyclic graph. For each node in the connected maximal subgraph of the resulting directed acyclic graph, a subcluster is formed. Both the initial cluster and the resulting subclusters are labeled as cross-sample consistent clusters of variant sites.

[0037] In specific implementation, please refer to Figure 5The process of performing local multiple alignments on sequences within each cross-sample consistent variant site cluster is as follows: A predetermined flanking length of 500 base pairs is extended outwards from the boundary of each cross-sample consistent variant site cluster. This length is determined based on the principle that extending 500 base pairs upstream and downstream of the variant site cluster completely covers the local genomic region where the variant site is located, while reserving sufficient flanking sequences for subsequent primer and probe design. The original sequencing read set within the extended interval for each sample is extracted. The original sequencing read set refers to all sequencing reads that overlap with the extended interval, extracted from the raw genome sequencing data of each sample according to the chromosome number and start / end coordinates of the extended interval using the samtoolsview command.

[0038] The raw sequencing read sets of all samples are merged and input into the progressive multiplex aligner. The progressive multiplex aligner employs an iterative progressive alignment algorithm, and its core architecture includes a distance matrix construction module, a guide tree construction module, and a progressive alignment module. The distance matrix construction module performs pairwise local alignments of all input raw sequencing reads, using a fast distance estimation method based on k-mer matching counts to obtain the sequence similarity distance between any two raw sequencing reads, generating a distance matrix. The guide tree construction module receives the distance matrix, calculates the branch length using an unweighted group averaging method, and constructs a binary guide tree, with leaf nodes corresponding to each raw sequencing read. The progressive alignment module performs sequence alignment operations from bottom to top along the guide tree. At each internal node of the guide tree, it performs a global alignment of the alignment groups corresponding to the two child nodes. The global alignment uses the Needleman-Wunsch dynamic programming algorithm, with a gap opening penalty of 10 and a gap extension penalty of 0.5. Each cross-sample consistent variant site cluster generates a corresponding local multiple alignment matrix. Each column of the local multiple alignment matrix corresponds to an alignment column position, which refers to the character column at the same alignment position of all reads in the multiple alignment. Each row of the local multiple alignment matrix corresponds to a read of one sample. The element in the r-th row and c-th column of the matrix represents the base character of the r-th read at the c-th alignment column position. The space character indicates that the read has a missing character at the alignment column position.

[0039] The process of identifying single nucleotide polymorphism (SNP) sites in a local multiple alignment matrix is ​​as follows: The local multiple alignment matrix is ​​scanned column by column, and the frequency of different base types in each column is counted. The base types include adenine, thymine, cytosine, and guanine, denoted as A, T, C, and G, respectively. If more than three different base types exist in a column, and the frequency of each base type exceeds the lower frequency limit (set to 0.1), then that column is identified as a SNP site. The lower frequency limit is set to 0.1 because base types with frequencies below 0.1 are usually introduced by PCR amplification errors or sequencing errors; setting the lower frequency limit to 0.1 effectively distinguishes genuine SNP signals from low-level sequencing noise. Consecutive single nucleotide polymorphism (SNP) sites are merged into sequence polymorphism (SPP) regions. The merging criteria are: two adjacent SNP sites are aligned in the local multiple alignment matrix, and the physical distance between the two SNP sites on the Leymus chinensis reference genome is less than 50 base pairs. Multiple SNP sites that meet the merging criteria are grouped into the same sequence polymorphism region.

[0040] The process of identifying simple repeating sequence regions in a local multiple alignment matrix is ​​as follows: A sliding window of 20 alignment columns is used to scan the local multiple alignment matrix, moving one alignment column at a time. Within each sliding window, the presence of insertion / deletion patterns that cyclically occur with a fixed repeat unit length is detected. The fixed repeat unit length ranges from 1 to 6 base pairs. The detection method is as follows: For each read sequence covered by the sliding window, the sequence length of each read within the window is statistically analyzed. A pattern is identified where the difference between the peak length of the occurrence frequency and its adjacent lengths shows a periodic variation. When the periodic variation period length is fixed and there are at least 5 complete periods, an insertion / deletion pattern with a fixed repeat unit length cyclically occurs within the sliding window. If the pattern shows a difference in the number of repetitions across different samples (i.e., the sequence length distribution of reads from different sample sources shows two or more separate peaks within this interval), then the interval corresponding to the sliding window is identified as a simple repeating sequence region. The difference in the number of repetitions is determined by comparing the peak positions of the sequence length distribution in different sample reads. The difference between the peak positions is divided by the length of the repetition unit to obtain the difference in the number of repetitions.

[0041] When designing molecular markers, for each simple repeat sequence region, conserved sequences without repeating units within a 200-base-pair length range flanking the simple repeat sequence region are extracted. The method for determining conserved sequences without repeating units is as follows: a sequence matching search is performed on the known repeating unit sequences within the 200-base-pair range flanking the simple repeat sequence region, excluding segments containing complete copies of the repeating unit; the remaining segments are considered conserved sequences. Primer design templates are selected from the conserved sequences, and the conserved sequences are input into the primer design software Primer3, with constraints set. The constraint for consistent primer annealing temperature is that the annealing temperature of all flanking primers is set to 60°C, with an allowable deviation range of ±1°C. The constraint for primers not forming hairpin structures is that the 3′ end complementarity score of the primers is less than 30; the complementarity score is the base pairing score within the primer molecule. The constraint for the length of the primer pair amplification product within the target range is that the target amplification product length range is set to 100 to 300 base pairs. Primer3 automatically searches for forward and reverse primer sequences that meet the above constraints on the primer design template, generating flanking primer pairs.

[0042] For each single nucleotide polymorphism (SNP) region, flanking sequences of 150 base pairs upstream and downstream of the SNP site were extracted as probe design templates. The extraction method involved extending 150 base pairs upstream and downstream from the center coordinates of the SNP site on the *Leymus chinensis* reference genome, and then extracting the sequence of that segment. Two competitive allele-specific probes were designed based on the base types of the two alleles at the SNP site. Both competitive allele-specific probes were of the same length, initially set to 20 nucleotides. The annealing temperature of the probes was calculated using the following formula:

[0043] in, The annealing temperature of the probe is in degrees Celsius. The value is the molar concentration of sodium ions in the reaction system, which is 0.05 mol / L. This value is based on the typical concentration of sodium ions in the standard PCR amplification buffer. This is the percentage of the sum of the number of guanine and cytosine bases in the probe sequence relative to the total number of bases in the probe sequence, and is a real number ranging from 0 to 100. The value is 20, representing the probe length in nucleotides. The annealing temperatures of the two competing allele-specific probes are calculated, ensuring the temperature difference between them does not exceed a temperature deviation threshold of 2°C. If the temperature difference exceeds 2°C, the probe length is adjusted by gradually increasing or decreasing the number of 5′ terminal bases within the range of 18 to 25 nucleotides, and the annealing temperature is recalculated until the temperature difference is less than or equal to 2°C. Locked nucleic acid (LCA) modification is applied to the probe ends near the single nucleotide polymorphism (SNP) site. The LCA modification is located at the penultimate or antepenultimate base of the 3′ end of the probe, meaning the probe base complementary to the allele base at the SNP site is modified with LCA. The ribose ring of the LCA-modified base is locked by a 2′-O,4′-C methylene bridge.

[0044] Each flanking primer pair was associated with its corresponding simple repeat sequence region, and each single nucleotide site region was associated with its corresponding two competing allele-specific probes. The association method involved mapping the genomic coordinates of the simple repeat sequence region and the single nucleotide site region to the flanking primer pair sequence and the two competing allele-specific probe sequences, respectively, using the genomic coordinates of both regions as keys in the data structure. All association results were summarized into a primary molecular marker set, which included the genomic location information, marker type, primer or probe sequence, and associated variant region information for each marker.

[0045] In specific implementation, please refer to Figure 6 Before performing simulated amplification validation of the primary molecular marker set using the raw genome sequencing data of all samples, the process of redundancy merging of the primary molecular marker set is as follows: The physical distance between any two markers in the primary molecular marker set on the Leymus chinensis reference genome is calculated. The physical distance is defined as the absolute value of the difference between the starting coordinates of the first marker and the second marker on the Leymus chinensis reference genome. Marker pairs with a physical distance less than a distance threshold are marked as neighboring marker pairs. The distance threshold is set to 1000 base pairs. The rationale for setting the distance threshold to 1000 base pairs is that when the genomic physical distance between two molecular markers is less than 1000 base pairs, the two molecular markers are usually located within the same linkage block in genetic linkage analysis, and the polymorphic information they provide is highly overlapping.

[0046] For each neighboring marker pair, the flanking primer sequences or probe sequences of the two markers are compared separately. When both markers in a neighboring marker pair are simple repeat markers, the forward and reverse primer sequences of each marker's flanking primer pair are extracted, and the sequence similarity between the forward primer sequence of the first marker and the forward primer sequence of the second marker is calculated. Simultaneously, the sequence similarity between the reverse primer sequence of the first marker and the reverse primer sequence of the second marker is also calculated. When both markers in a neighboring marker pair are single nucleotide site markers, two competitive allele-specific probe sequences of each marker are extracted, and the sequence similarity between the two probe sequences of the first marker and the two probe sequences of the second marker is calculated. When one marker in a neighboring marker pair is a simple repeat marker and the other is a single nucleotide site marker, only the sequence similarity between the flanking primer sequences and probe sequences is calculated. Sequence similarity was calculated using a global alignment algorithm. Two sequences to be aligned were input into a Needleman-Wunsch dynamic programming algorithm, with a matching score of 1, a mismatch penalty of -1, a gap opening penalty of -2, and a gap extension penalty of -0.5. The optimal alignment score was obtained, and the sequence similarity was calculated by dividing the optimal alignment score by the length of the shorter sequence. If the sequence similarity of the flanking primer or probe sequences of two markers exceeded a similarity threshold, the two markers were considered a redundant marker pair. The similarity threshold was set to 0.8. The rationale for setting the similarity threshold to 0.8 is that when the sequence similarity exceeds 0.8, the flanking primers or probes of the two markers are prone to cross-hybridization during amplification, producing non-specific amplification products, making it difficult to independently distinguish the two markers in the experiment.

[0047] Markers with low polymorphism information content (PIP) in all redundant marker pairs are removed from the primary molecular marker set. The calculation method for PIP is described in the subsequent simulated amplification validation step. For redundant marker pairs without PIP data (i.e., the markers in the redundant marker pairs have not yet been validated by simulated amplification), a simplified PIP estimation is performed on each of the two markers in the redundant marker pair. The simplified PIP estimation method is as follows: quickly extract the base composition variation of each region of the two markers from the raw genome sequencing data of all samples, count the read support number of the reference allele and the substitute allele, calculate the allele frequency based on the read support number, substitute the allele frequency into the PIP calculation formula to obtain an approximate PIP value, and use the approximate PIP value as the basis for removal. The markers with low polymorphism information content in each redundant marker pair are removed from the primary molecular marker set. The remaining markers are then sorted in the primary molecular marker set according to their genomic location. The sorting method is as follows: first, they are sorted in ascending order by chromosome number, and within the same chromosome, they are sorted in ascending order by the starting coordinate of the marker. The redundant primary molecular marker set is then regenerated.

[0048] The process of simulating amplification and validation of the deredundant primary molecular marker set using raw genome sequencing data from all samples is as follows: For each marker in the deredundant primary molecular marker set, reads of the region containing the marker are extracted from the raw genome sequencing data of all samples, extending 200 base pairs upstream and downstream from the marker's genomic coordinates. The simulation calculation is performed using a read-matching-based computational framework. The input to the simulation amplification algorithm is the marker's flanking primer pair sequences or competitive allele-specific probe sequences, the read set of the region containing the marker extracted from all samples, and the reference sequence of the region containing the marker on the Leymus chinensis reference genome. The output of the simulation amplification algorithm is the amplified read set for each sample and the genotype determination result for each sample.

[0049] When the marker type is a simple repeat sequence marker, the simulated amplification process is as follows: The forward and reverse primer sequences of the flanking primer pair are locally aligned with each read using the Smith-Waterman algorithm, allowing a maximum of two mismatches. When the forward primer sequence aligns to the 5' end of the read within the allowed mismatch range, and the reverse complementary sequence of the reverse primer aligns to the 3' end of the same read, the sequence fragment between the 3' end of the forward primer and the 3' end of the reverse complementary sequence is extracted as the simulated amplification product. The sequence length of the simulated amplification product is calculated. If the sequence length falls within the target range (100 to 300 base pairs), the sample is considered to have successfully amplified the product of the expected length. When the label type is single nucleotide site labeling, the simulated amplification process is as follows: two competitive allele-specific probes are locally aligned with each read. The alignment uses the Smith-Waterman algorithm, which allows a maximum of one mismatch and does not allow mismatches to occur within the last 5 bases of the 3′ end of the probe. When one of the two competitive allele-specific probes can align with a read, the read is identified as the amplification product of the allele corresponding to that probe, and the read abundance of the amplification product is recorded.

[0050] The number of samples that successfully amplified the product of the expected length was counted. A successfully amplified sample was defined as one containing at least one read from the simulated amplification product. The ratio of the number of successfully amplified samples to the total number of samples was used as the amplification success rate for the marker. For each marker, the genotype of each sample was calculated based on the abundance of different allele reads in the simulated amplification product. When the marker was a simple repeat sequence marker, the allele type was determined by the copy number of the repeat unit in the simple repeat sequence region. The number of reads corresponding to different repeat unit copy numbers in the simulated amplification product was counted, and the repeat unit copy number with the highest read number was taken as the genotype of the sample. When the marker was a single nucleotide polymorphism (SNP) marker, the allele type was determined by the base type of the SNP site. The number of reads corresponding to two competing allele-specific probes was counted. A homozygous genotype was defined when the ratio of the number of reads corresponding to the two probes was greater than 3, and a heterozygous genotype was defined when the ratio was between 0.33 and 3. The types of alleles observed in all samples and the frequency of each allele were counted. The number of allele types is the total number of different alleles appearing in all samples, and the allele frequency is the number of samples containing a particular allele divided by the total number of samples. The polymorphism information content (PIM) values ​​were obtained by substituting the allele frequencies into the polymorphism information content formula. The polymorphism information content calculation formula is as follows:

[0051] in, The numerical value representing the polymorphic information content of the tag. The value of is a real number between 0 and 0.5. A higher value indicates a higher degree of polymorphism in the tag. This represents the total number of allele types observed in all samples. It is an integer greater than or equal to 2. Indicates the first Frequency of alleles equal to the The number of samples in which each allele appeared divided by the total number of samples. The value of is a real number between 0 and 1. Indicates the first Frequency of alleles This is the index number of the allele. The value range is 1 to integers, and Different values, equal to the The number of samples in which each allele appears is divided by the total number of samples. The first term of the formula is the summation term. The second term represents the sum of squares of all allele frequencies. It represents twice the sum of the squared products of the frequencies of all pairs of alleles.

[0052] Markers with an amplification success rate greater than a success rate threshold of 0.8 are considered valid for amplification validation. This threshold is based on the principle that a marker possesses sufficient experimental stability and sample applicability when it can be successfully amplified in over 80% of samples. Markers with a polymorphism information content (PI) value greater than a PI threshold of 0.25 are considered valid for polymorphism. This threshold is based on the principle that PI values ​​above 0.25 indicate moderate to high polymorphism, providing effective information for genetic diversity analysis, population structure analysis, and linkage map construction. Markers with both an amplification success rate and a PI value greater than the success rate threshold are considered valid for validation. All valid markers are output as *Leymus chinensis* polymorphic molecular markers, including their genomic coordinates, marker type, flanking primer pair sequence or competitive allele-specific probe sequence, amplification success rate, and PI value.

[0053] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing, characterized in that, Includes the following steps: We obtained raw genome sequencing data from multiple samples of Leymus chinensis. After quality filtering and error correction of the raw data for each sample, we constructed reference sequence alignment results for each sample. Candidate regions with inconsistent coverage depths are extracted from the alignment results of each sample, and a variation confidence score matrix is ​​established based on the base substitution frequency and insertion / deletion length distribution within the candidate regions. Based on the variance confidence score matrix, collinear clustering is performed on all variance sites in all samples to generate consistent variance site clusters across samples; Local multiple alignments were performed on sequences within each cluster of variant sites to identify simple repeat sequence regions with length polymorphism and single nucleotide site regions with sequence polymorphism. Flanking primer pairs were designed for simple repetitive sequence regions, and competitive allele-specific probes were designed for single nucleotide site regions, forming a set of primary molecular markers; The primary molecular marker set was amplified and validated using the raw genome sequencing data of all samples. Markers with amplification success rate and polymorphism information content exceeding the threshold were selected as polymorphic molecular markers of Leymus chinensis.

2. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 1, characterized in that, The steps for constructing the reference sequence alignment results for each sample include: The quality-filtered and error-corrected sequencing reads of each sample were initially aligned with the Leymus chinensis reference genome to obtain the initial alignment position of each read. Based on the initial alignment position, the distribution of the pairing end distance of each read segment is statistically analyzed, and abnormal read segments with pairing end distances exceeding three standard deviations are removed; The remaining read segments are used to perform local re-alignment on regions that were not aligned or whose alignment quality was below the quality threshold in the initial alignment, generating the final alignment record for each sample. The final alignment record is then sorted by coordinates and output as the reference sequence alignment result for each sample.

3. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 1, characterized in that, The step of extracting candidate regions with inconsistent coverage depth from the alignment results of each sample includes: The reference sequence alignment results of each sample are scanned with a preset window length and step size, and the average coverage depth and standard deviation of coverage depth within each window are calculated. Windowings whose average coverage depth deviates from the average coverage depth of the whole genome by more than a fold deviation threshold are designated as depth aberration windows, and base substitution events and insertion / deletion events within the depth aberration windows are recorded; Adjacent depth anomaly windows are merged into continuous intervals. When the length of a continuous interval is greater than the lower limit of the length and the base substitution frequency within the interval is higher than the frequency threshold, the continuous interval is identified as a candidate region with inconsistent coverage depth.

4. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 1, characterized in that, The step of establishing a variation confidence score matrix based on the base substitution frequency and insertion / deletion length distribution within the candidate region includes: For each base position within each candidate region, count the number of times a non-reference base appears at that position in all samples, and calculate the non-reference base frequency as the base substitution frequency. For each insertion or missing event within each candidate region, record the length of the inserted or missing sequence, and count the number of occurrences of the insertion or missing event in all samples to construct a histogram of the insertion or missing length distribution. The initial substitution confidence score is obtained by multiplying the substitution frequency at each base position by the reciprocal of the corresponding sequencing base quality value. The product of the mode length of the histogram of the insertion missing length distribution for each insertion missing event and the number of times the event occurs is used as the initial score for the insertion missing confidence score. Arrange the preliminary confidence scores for permutations and insertions / missings within all candidate regions according to the sample dimension and the location dimension to form a variation confidence score matrix.

5. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 4, characterized in that, The value of each element in the variation confidence score matrix is ​​a weighted sum of the initial confidence score for permutation and the initial confidence score for insertion / deletion.

6. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 1, characterized in that, The step of performing collinear clustering of all sample variant sites based on the variant confidence score matrix to generate consistent variant site clusters across samples includes: Extract the shared variant sites from the variant confidence score matrix, and represent each variant site as a vector, with each element of the vector corresponding to the confidence score of the sample at that site; Calculate the cosine similarity between any two variant site vectors, and group variant sites with a cosine similarity greater than the similarity threshold into the same initial cluster; For each initial cluster, the order of the genomic coordinates of the variant sites within the cluster is extracted. If the coordinate order is consistent across all samples, the initial cluster is retained; otherwise, the initial cluster is split into multiple sub-clusters according to the degree of conflict in the coordinate order. Both the retained initial clusters and the sub-clusters obtained from the split are labeled as cross-sample consistent variant site clusters.

7. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 6, characterized in that, The similarity threshold for the cosine similarity is set to 0.85, 0.90, or 0.

95.

8. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 1, characterized in that, The step of performing local multiple alignments of sequences within each cluster of variant sites to identify simple repeat sequence regions with length polymorphism and single nucleotide site regions with sequence polymorphism includes: Extend a predetermined flanking length outward from the boundary of each variant site cluster, and extract the set of original sequencing reads for each sample within the extended interval; The raw sequencing reads of all samples are input into the progressive multiple alignment unit to generate a local multiple alignment matrix corresponding to each cluster of variant sites. Each column of the matrix corresponds to an alignment column position, and each row corresponds to a read of a sample. Scan each column of the local multiple alignment matrix. If there are three or more different base types in a column and the frequency of each type exceeds the lower frequency limit, then the column is identified as a single nucleotide polymorphism site. Consecutive single nucleotide polymorphism sites are merged into a sequence polymorphism region. Scanning multiple consecutive columns of the local multiple alignment matrix, if there is an insertion / deletion pattern that appears cyclically with a fixed repeating unit length, and the number of repetitions of this pattern varies in different samples, then the interval corresponding to the multiple consecutive columns is determined to be a simple repeating sequence region.

9. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 1, characterized in that, The steps of designing flanking primer pairs for simple repetitive sequence regions and competitive allele-specific probes for single nucleotide site regions to form a primary molecular marker set include: For each simple repetitive sequence region, conserved sequences without repetitive units are extracted from both sides of the region within a preset length range. Primer design templates are selected from the conserved sequences, and flanking primer pairs are generated according to the constraints of consistent primer annealing temperature, primers not forming hairpin structures, and primer amplification product length within the target range. For each single nucleotide site region, flanking sequences of predetermined lengths upstream and downstream of the site are extracted as probe design templates. Two competitive allele-specific probes are designed according to the base types of the two alleles at the site. The two probes are of the same length and the annealing temperature difference does not exceed the temperature deviation threshold. Locked nucleic acid modification is used at the end of the probe near the single nucleotide site. Each flanking primer pair was associated with its corresponding simple repeat sequence region, and each single nucleotide site region was associated with its corresponding two competing allele-specific probes. All association results were summarized into a primary molecular marker set.

10. The method for developing polymorphic molecular markers of Leymus chinensis based on genome sequencing according to claim 9, characterized in that, The locked nucleic acid modification is located at the penultimate or tertiary base of the 3′ end of the probe.

Images