High-throughput genotype intelligent analysis method

Abstract

The invention discloses a high-throughput genotype intelligent analysis method. The method comprises the steps of data acquisition and preprocessing, molecular marker recognition and genetic relationship and hybridization advantage prediction. Performing standardization processing on the high-throughput genotype data, establishing triple quality evaluation standards of coverage rate, conversion / transversing rate and error rate, and correcting batch effect through a position effect index; then calculating a genetic similarity matrix among the samples based on the preprocessed data, and analyzing and identifying a tag SNP with high centrality as a molecular marker through a linkage imbalance network; the markers are used for estimating the genetic relationship of the population and constructing a phylogenetic tree, and meanwhile, the hybridization advantage is predicted based on the relationship index of the heterozygosity and the genetic distance. According to the method, the genotype data quality is remarkably improved, molecular marker screening is optimized, the population genetic relationship is accurately estimated, the hybridization advantage is accurately predicted, and an efficient bioinformatics solution is provided for modern breeding.

Description

Technical Field

[0001] The present invention belongs to the technical field of breeding analysis, and in particular relates to a high-throughput genotype intelligent analysis method. Background Art

[0002] With the rapid development of genomics and information technology, various high-throughput sequencing platforms such as Illumina, PacBio, and Oxford Nanopore have emerged continuously. The amount of data generated by a single sequencing run has rapidly increased from GB to TB, providing unprecedented amounts of genomic information for biological breeding. This large-scale genomic data provides strong support for breeding technologies such as whole-genome selection (Genomic Selection), molecular-assisted breeding (Marker-Assisted Selection), and genetic diversity assessment.

[0003] However, the genotype data formats from different sequencing platforms are not uniform, the data quality is uneven, and there is a lack of effective quality control standards and evaluation systems. As the number of sequencing batches increases, the systematic differences between batches (Batch Effect) become increasingly significant, posing severe challenges to data integration and analysis. Traditional SNP (single nucleotide polymorphism) screening methods, such as p-value-based screening or fixed-interval sampling, cannot fully consider the linkage disequilibrium (Linkage Disequilibrium) structure between SNP sites, resulting in information redundancy or loss of key variants. At the same time, the computational complexity of simultaneously analyzing millions of SNP sites in large-scale populations is too high, seriously affecting analysis efficiency. In terms of population kinship estimation, traditional methods such as the Van Raden method or the GCTA-GRM method often ignore the influence of population structure (Population Structure), resulting in biased kinship estimation. This bias is particularly significant in the presence of obvious population stratification, which is not conducive to breeders in line selection and hybrid design.

[0004] Therefore, an intelligent high-throughput genotyping method is urgently needed to overcome the above technical bottlenecks and improve breeding efficiency and breeding accuracy. Summary of the Invention

[0005] The present invention provides a high-throughput intelligent genotyping analysis method, which improves the quality of genotyping data, optimizes molecular marker screening, accurately estimates population relatedness, and accurately predicts hybrid vigor by introducing multiple quality assessment criteria, batch effect correction, linkage disequilibrium network analysis, and hybrid vigor relationship index modeling.

[0006] In order to achieve the above-mentioned purpose of the invention, the specific technical solutions are as follows: A high-throughput genotype intelligent analysis method, comprising the following steps: Step S1: collect high-throughput genotyping data, standardize the sample data format, perform quality assessment on the sample data and remove batch effects to obtain preprocessed genotyping data.

[0007] The method for quality assessment of sample data includes: calculating sample genotype coverage, calculating transition / transversion ratio in genotype data, and calculating genotype error rate.

[0008] Sample screening was performed based on quality assessment results, excluding genotype coverage. Sample.

[0009] Knockout transition / transversion ratio Abnormal deviation from the population A sample of and are the expectation and variance of the population transition / transversion ratio, respectively.

[0010] Elimination of genotyping error rate Sequencing batch data.

[0011] Step S2, using the pre-processed genotype data, calculate the genetic similarity matrix between samples and perform linkage disequilibrium analysis to identify tag SNPs as molecular markers.

[0012] Step S3, estimating population kinship and predicting hybrid vigor based on the obtained molecular markers.

[0013] Furthermore, the standardized sample data format is specifically as follows: The genotype data obtained from different sequencing platforms are uniformly converted into a standard format. Numerical representation, where 0 represents the homozygous major allele AA, 1 represents the heterozygous genotype Aa, and 2 represents the homozygous minor allele aa.

[0014] The missing data were filled using the multiple imputation method based on the local linkage disequilibrium LD structure. The calculation formula is: ,in Represents a sample At the site The genotype value to be imputed on Representation and location At the LD threshold The set of all non-deleted sites within the range, Indicates the site and site The weight between , For the site and site The LD coefficient between them.

[0015] Furthermore, the sample genotype coverage The calculation formula is: ,in, For samples The number of effective genotype sites in is the total number of sites.

[0016] Transition / transversion ratios in genotype data The calculation formula is: ,in, is the number of conversion mutations. The conversion mutations include: A G or C T, is the number of transversion mutations, including: A C, A T, G C or G T. The left and right represent the genes before and after mutation, respectively.

[0017] Genotyping error rate The calculation formula is: ,in, is the number of repeated sequencing samples, For repeated samples The number of inconsistent sites, For repeated samples The total number of sites in the CCP.

[0018] Furthermore, the method for removing batch effects includes: Constructing the position effect index ,in Indicates batch midpoint The average genotype of and Sites The mean and standard deviation across all batches, For batch Sequencing time, in days, and are the mean and standard deviation of sequencing time for all batches, respectively.

[0019] right The sites are batch corrected, and the correction formula is: ,in For the site The batch effect coefficient is estimated by minimizing the corrected between-batch variance: , Indicates the site A vector of genotype values ​​across all batches.

[0020] Furthermore, the calculation of the genetic similarity matrix between samples includes the following steps: Based on the preprocessed genotype data, for any two samples and , calculate the sample and samples At the site Correlation on : When the genotypes of two samples are exactly the same , when they share an allele , completely different at the same time .

[0021] Calculate the sum of the associations between two samples at all common non-missing sites: ,in is the number of common non-missing SNP sites.

[0022] Calculate genetic similarity weighted by allele frequencies: ; in, , is the site The weight of For the site minor allele frequency.

[0023] Furthermore, the linkage disequilibrium analysis is specifically as follows: using the pre-processed genotype data, the linkage disequilibrium coefficient between any two SNP sites is calculated. : ,in, , is the frequency of haplotype AB, and are the frequencies of the major alleles at sites A and B, respectively, and are the frequencies of the minor alleles at sites A and B, respectively.

[0024] Construct a linkage disequilibrium decay model: ,in, is the effective population size, and d is the genetic distance between two loci in Morgan units.

[0025] Linkage disequilibrium half-distance Evaluate the extent of linkage disequilibrium in a population, i.e. When , the linkage disequilibrium half-decline distance is: .

[0026] Furthermore, the identification tag SNP is used as a molecular marker specifically to construct a SNP site network based on linkage disequilibrium analysis , where the vertex set V represents all SNP sites, and the edge set E represents the linkage disequilibrium relationship. Greater than threshold When , there is an edge between these two sites.

[0027] Calculate the degree centrality of each SNP site: ,in is the indicator function.

[0028] For each connected component, the SNP site with the highest degree centrality is selected as the tag SNP of the connected component.

[0029] If the degree centrality of multiple sites is the same, the site with the minor allele frequency closest to 0.5 is selected as the tag SNP.

[0030] Calculate the capture rate of tag SNPs: , where T is the tag SNP set, is the total number of SNP sites.

[0031] Furthermore, the estimation of population kinship based on the obtained molecular markers includes: using the identified tag SNP as a molecular marker to calculate the common ancestor coefficient matrix between samples : ,in is the number of tag SNPs, Represents a sample In the The genotype value on the tag SNP, For the The minor allele frequency of each tag SNP.

[0032] Computing the kinship matrix for population structure correction ,in is the group structure matrix, obtained through the group structure analysis software STRUCTURE, is the diagonal matrix of differentiation coefficients among different subpopulations.

[0033] Construct a phylogenetic tree based on the kinship matrix to construct an unrooted tree, and the branch length of the tree reflects the genetic distance between samples: .

[0034] Furthermore, the method for predicting hybrid vigor is specifically as follows: calculating the expected hybrid F1 heterozygosity of the parental genotypes : ,in, and Respectively represent parents and parents At the site The genotype value on .

[0035] Calculate the relationship index between the genetic distance between parents and hybrid vigor: ,in For parents and The genetic distance between them.

[0036] Constructing a hybrid vigor prediction model: ,in is the predicted hybrid vigor value, , , and is the regression coefficient, which is estimated by the least squares method based on the known hybrid combination phenotypic data and the corresponding molecular marker data.

[0037] Compared with the prior art, the present invention has the following beneficial effects: This invention significantly improves the reliability of genotype data by establishing triple quality assessment criteria for coverage, transition / transversion ratio, and error rate. It also introduces a position effect index to assess batch effects and effectively eliminates systematic biases between sequencing batches through an adaptive correction algorithm. The corrected inter-batch variance is reduced by 85%. Linkage disequilibrium network analysis identifies tag SNPs, reducing the number of molecular markers by 70% while maintaining a capture rate of over 90%, significantly reducing the computational complexity of subsequent analysis. The invention constructs a hybrid vigor relationship index based on heterozygosity and genetic distance, achieving a hybrid vigor prediction accuracy of over 85%, providing precise guidance for efficient breeding. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 This is a flow chart of a high-throughput genotype intelligent analysis method of the present invention. DETAILED DESCRIPTION

[0039] To make the objectives, technical solutions, and advantages of the present invention more clear, the technical solutions of the present invention are described clearly and completely below. Obviously, the embodiments described are only part of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0040] like Figure 1 FIG. 1 is a high-throughput genotyping intelligent analysis method of the present invention, comprising the following steps: Step S1: collect high-throughput genotyping data, standardize the sample data format, perform quality assessment on the sample data and remove batch effects to obtain preprocessed genotyping data.

[0041] High-throughput genotyping data can be derived from a variety of sequencing platforms, including but not limited to the Illumina HiSeq / NovaSeq series, the PacBio Sequel system, and the Oxford Nanopore PromethION. The raw data formats generated by each platform vary. For example, Illumina typically outputs data in FASTQ format, which may be in VCF or BCF format after genotyping analysis; PacBio may output data in BAM format; and Oxford Nanopore may output data in FAST5 format. This method first converts these different formats into a table containing sample ID, chromosome location, reference allele, variant allele, and genotype value; the standardized sample data format is specifically: The genotype data obtained from different sequencing platforms are uniformly converted into a standard format. Numerical representation, where 0 represents the homozygous major allele AA, 1 represents the heterozygous genotype Aa, and 2 represents the homozygous minor allele aa.

[0042] For example, for VCF files from the Illumina platform, tools such as VCFtools or bcftools can be used to extract SNP information and convert it to a standard format. For long-read data from PacBio or Oxford Nanopore, tools such as GATK HaplotypeCaller must first be used to perform SNP detection and then standardize the format. Normalized data typically includes columns such as chromosome number (Chr), physical position (Pos), reference allele (Ref), variant allele (Alt), and genotype values ​​(0 / 1 / 2) for each sample, facilitating data processing and analysis in environments such as Python and R.

[0043] The missing data were filled using the multiple imputation method based on the local linkage disequilibrium LD structure. The calculation formula is: ,in Represents a sample At the site The genotype value to be imputed on Representation and location At the LD threshold The set of all non-deleted sites within the range, Indicates the site and site The weight between , For the site and site The LD coefficient between them.

[0044] Traditional methods such as mean or median imputation fail to account for linkage relationships between loci in the genome. This method, however, leverages the local LD ​​structure for interpolation, enabling more accurate estimation of missing values. In practice, the LD coefficient r² is first calculated between all locus pairs to construct an LD matrix. Then, for each missing locus, the set of neighboring loci within the LD threshold (r² > 0.6) is determined. The genotypic value of the missing locus is estimated based on the weighted average of the genotype values ​​and LD strength of these neighboring loci. This method is particularly suitable for species with slow LD decay, such as self-fertile crops (such as rice and wheat), and can achieve interpolation accuracy exceeding 90%.

[0045] The method for quality assessment of sample data includes: calculating sample genotype coverage, calculating transition / transversion ratio in genotype data, and calculating genotype error rate.

[0046] Sample screening was performed based on quality assessment results, excluding genotype coverage. Low coverage often means insufficient sequencing depth or poor DNA quality. The 85% coverage threshold set in this method is an empirical value derived from the analysis of a large amount of experimental data, and performs well in balancing data integrity and sample retention. For example, in a genotype analysis involving 1,000 rice varieties, setting an 85% coverage threshold eliminated approximately 5% of low-quality samples, which may introduce noise in subsequent analyses, while the remaining 95% of samples contained sufficient information for downstream analysis.

[0047] Knockout transition / transversion ratio Abnormal deviation from the population A sample of and are the expectation and variance of the population transition / transversion ratio, respectively.

[0048] Elimination of genotyping error rate sequencing batch data; setting a 5% error rate threshold takes into account the performance characteristics and cost-effectiveness of high-throughput sequencing platforms under the current technological level; taking crop molecular marker-assisted breeding as an example, when the error rate exceeds 5%, the reliability of genotype data for parent identification and kinship analysis will be significantly reduced, which may lead to breeding decision-making errors; each sequencing batch contains at least 2-3 replicate samples to accurately calculate the error rate.

[0049] Sample genotype coverage The calculation formula is: ,in, For samples The number of effective genotype sites in is the total number of sites.

[0050] Transition / transversion ratios in genotype data The calculation formula is: ,in, is the number of conversion mutations. The conversion mutations include: A G or C T, is the number of transversion mutations, including: A C, A T, G C or G T. The left and right represent the genes before and after mutation, respectively.

[0051] Genotyping error rate The calculation formula is: ,in, is the number of repeated sequencing samples, For repeated samples The number of inconsistent sites, For repeated samples The total number of sites in the CCP.

[0052] The method for removing batch effects includes: constructing a position effect index ,in Indicates batch midpoint The average genotype of and Sites The mean and standard deviation across all batches, For batch Sequencing time, in days, and are the mean and standard deviation of the sequencing time for all batches, respectively. The construction principle of this index is based on the observation that batch effects often show systematic changes with sequencing time; when the LI value is greater than 2.5, it indicates that the site has a significant systematic deviation in a specific batch. Taking a soybean genome sequencing project involving 20 batches and lasting two years as an example, about 15% of SNP sites showed obvious batch effects, mainly concentrated in areas with high GC content. After applying the batch correction of this method, the systematic differences between different batches were significantly reduced. The results of principal component analysis (PCA) before and after correction showed that the batch clustering phenomenon was basically eliminated.

[0053] right The sites are batch corrected, and the correction formula is: ,in For the site The batch effect coefficient is estimated by minimizing the corrected between-batch variance: , Indicates the site The vector of genotype values ​​in all batches. The batch correction coefficient is estimated by minimizing the corrected inter-batch variance, which is essentially a correction method based on linear regression. The coefficient reflects the trend strength of site j over sequencing time. The larger the value, the more significantly the site is affected by sequencing time.

[0054] For example, in a whole-genome sequencing analysis of rice, we found that approximately 8% of SNP sites showed significant linear time trends, and after correction, the inter-batch variance of these sites was reduced by an average of 78%.

[0055] Step S2, using the pre-processed genotype data, calculate the genetic similarity matrix between samples and perform linkage disequilibrium analysis to identify tag SNPs as molecular markers.

[0056] The calculation of the genetic similarity matrix between samples comprises the following steps: Based on the preprocessed genotype data, for any two samples and , calculate the sample and samples At the site Correlation on : When the genotypes of two samples are exactly the same , when they share an allele , completely different at the same time .

[0057] If two samples share a rare variant allele at a site with a major allele frequency of 0.95, this indicates that they are more closely related than if they share a site with an allele frequency of 0.5; the genetic similarity matrix calculated by weighted method shows an accuracy of about 15-20% higher than the simple matching coefficient method in distinguishing closely related samples and identifying parents.

[0058] The site-level association scoring scheme, based on the IBS (Identity By State) principle, intuitively reflects the degree of similarity between samples at the molecular level. For example, for a SNP locus, if the genotype of sample A is AA (coded as 0) and sample B is also AA, then S = 1; if sample B is AG (coded as 1), they share an A allele, S = 0.5; if sample B is GG (coded as 2), the two samples do not share any alleles, S = 0. This three-level scoring scheme is suitable for diploid species.

[0059] For polyploid species (such as wheat, cotton, etc.), it can be expanded to a more detailed multi-level score. In an analysis of hexaploid wheat, the correlation was subdivided into seven levels: 0, 0.17, 0.33, 0.5, 0.67, 0.83 and 1.0, which more accurately captured the strength of the kinship between different samples.

[0060] Calculate the sum of the associations between two samples at all common non-missing sites: ,in is the number of common non-missing SNP sites.

[0061] Calculate genetic similarity weighted by allele frequencies: ;in, , is the site The weight of For the site The minor allele frequency (pk) of a locus is considered to be high. Sites with minor allele frequencies (pk) close to 0 or 1 are weighted higher because variants at these sites are rare and less likely to be shared between samples, thus containing more information about kinship relationships. For example, in a population genomics study, a rare variant with a frequency of 0.01 may have several times greater discriminatory power than a common variant with a frequency of 0.5. To avoid excessive weights due to extreme frequencies (such as pk close to 0), a minimum frequency threshold (such as 0.01) can be set, or the empirical Bayesian method can be used to smooth the weights.

[0062] Constructing sample genetic similarity matrix , yes dimensional matrix, is the number of samples after preprocessing, The linkage disequilibrium analysis is specifically as follows: using the pre-processed genotype data, calculate the linkage disequilibrium coefficient between any two SNP sites : ,in, , is the frequency of haplotype AB, and are the frequencies of the major alleles at sites A and B, respectively, and are the frequencies of the minor alleles at site A and site B, respectively; r²=0 means the two sites are completely independent, and r²=1 means complete linkage.

[0063] Construct a linkage disequilibrium decay model: ,in, is the effective population size, and d is the genetic distance between two loci in Morgan units.

[0064] Linkage disequilibrium half-distance Evaluate the extent of linkage disequilibrium in a population, i.e. When , the linkage disequilibrium half-decline distance is: In a comparative study of multiple crops, the half-life distance of indica rice varieties was about 123 kb, and that of japonica rice was about 167 kb, reflecting that the japonica rice population experienced a stronger domestication bottleneck; the half-life distance of maize inbred line populations was about 10-30 kb, while the half-life distance of maize wild relatives was only 2-5 kb, reflecting the significant changes in effective population size during domestication and breeding.

[0065] Identifying tag SNPs as molecular markers involves: Constructing a SNP locus network based on linkage disequilibrium analysis , where the vertex set V represents all SNP sites, and the edge set E represents the linkage disequilibrium relationship. Greater than threshold When , there is an edge between these two sites.

[0066] Calculate the degree centrality of each SNP site: ,in is the indicator function.

[0067] For each connected component, the SNP site with the highest degree centrality is selected as the tag SNP of the connected component.

[0068] If the degree centrality of multiple sites is the same, the site with the minor allele frequency closest to 0.5 is selected as the tag SNP.

[0069] Calculate the capture rate of tag SNPs: , where T is the tag SNP set, is the total number of SNP sites.

[0070] The capture rate (CR) is a key metric for evaluating the representativeness of tag SNPs. It quantifies the extent to which a collection of tag SNPs represents the genetic variation at all SNP loci. A CR value closer to 1 indicates a better representativeness of the tag SNPs. In practical applications, there is often a trade-off between the number of tag SNPs (cost) and the capture rate (information content). Based on experience, when the threshold τ is set to 0.8, a CR of 0.85-0.95 is typically achieved, while reducing the number of markers by 70-85%.

[0071] Visualizing the relationship between capture rate and number of tag SNPs can help determine the optimal balance. For example, in a soybean breeding project, an initial pool of 420,000 SNPs was screened using a threshold of τ = 0.8, resulting in 36,000 tag SNPs and a capture rate of 0.91. Further relaxing the threshold to τ = 0.7 increased the number of tag SNPs to 48,000, raising the capture rate to 0.95, but at a 33% cost increase. Ultimately, the τ = 0.8 solution was chosen as the optimal balance.

[0072] Step S3, estimating population kinship and predicting hybrid vigor based on the obtained molecular markers.

[0073] The method of estimating population kinship based on the obtained molecular markers includes: using the identified tag SNP as a molecular marker to calculate the common ancestor coefficient matrix between samples : ,in is the number of tag SNPs, Represents a sample In the The genotype value on the tag SNP, For the The minor allele frequency of each tag SNP.

[0074] Computing the kinship matrix for population structure correction ,in is the group structure matrix, obtained through the group structure analysis software STRUCTURE, is the diagonal matrix of differentiation coefficients among different subpopulations.

[0075] Construct a phylogenetic tree based on the kinship matrix to construct an unrooted tree, and the branch length of the tree reflects the genetic distance between samples: .

[0076] The method for predicting hybrid vigor is specifically as follows: calculating the expected hybrid F1 heterozygosity of the parental genotypes : ,in, and Respectively represent parents and parents At the site The genotype value on .

[0077] Calculate the relationship index between the genetic distance between parents and hybrid vigor: ,in For parents and The genetic distance between them.

[0078] Constructing a hybrid vigor prediction model: ,in is the predicted hybrid vigor value, , , and is the regression coefficient, which is estimated by the least squares method based on the known hybrid combination phenotypic data and the corresponding molecular marker data.

[0079] According to biological knowledge, the quadratic term It is usually a negative value, reflecting the potential decrease in fitness caused by excessive heterozygosity. Regarding model fitting, it is recommended to use cross-validation methods, such as 10-fold cross-validation or leave-one-out cross-validation, to evaluate predictive performance. For example, in a study of maize heterosis, this model was applied to predict the yield of 380 hybrid combinations. The cross-validated prediction accuracy reached 0.72, significantly higher than the 0.58-0.65 of traditional prediction methods. The model also performed stably under different environmental conditions, verifying its practicality and reliability.

[0080] An example of the application of the present method in rice breeding is as follows: Genome resequencing data for 520 rice varieties (including indica, japonica, and intermediate varieties) from various regions were collected. Sequencing depths ranged from 15-30X, with coverage exceeding 96%, and the total amount of raw sequencing data reached 35TB. The data were generated using 28 sequencing batches across three different sequencing platforms (Illumina HiSeq 2500, NovaSeq 6000, and MGI DNBSEQ-T7) over a period of three years.

[0081] Applying the data preprocessing method of the present invention, data from different platforms were first converted to a standardized format, resulting in the detection of approximately 4.3 million high-quality SNPs. Quality assessment revealed that 18 samples had genotype coverage below 85%, 12 samples had transition / transversion ratios (average value of 2.14 ± 0.18) that deviated from the population by more than 3 standard deviations, and three sequencing batches had genotype error rates exceeding 5%. These low-quality data were removed, leaving 487 valid samples.

[0082] Applying a batch effect correction method, we identified significant batch effects (LI > 2.5) for approximately 12.5% ​​of the 538,000 SNPs, primarily concentrated between batches spanning a long sequencing time span. After correction, the coefficient of variation (CV) for these sites across different batches decreased from an average of 0.28 to 0.04, and batch clustering was largely eliminated in PCA analysis.

[0083] Linkage disequilibrium analysis revealed that the LD half-life distance of the rice population was approximately 150 kb, with the indica subpopulation (approximately 123 kb) being shorter than the japonica subpopulation (approximately 167 kb), reflecting historical differences between the subpopulations. A SNP locus network was constructed using an r² threshold of >0.8, resulting in approximately 465,000 connected components, each containing an average of 9.2 SNPs. The most representative loci within each connected component were selected based on degree centrality, ultimately identifying 48,723 tagging SNPs with a capture rate (CR) of 0.93, reducing the number of molecular markers by 88.7% while preserving the majority of genetic variation.

[0084] Based on these tag SNPs, the genetic relationships between samples were calculated. After correction for population structure, a phylogenetic tree was constructed, which accurately reflected the genetic relationships between rice varieties, with a 92.5% consistency with known pedigree information. In a simulated hybridization experiment, the expected heterozygosity and hybrid vigor relationship index were calculated for all possible hybrid combinations (approximately 118,000 combinations). Combined with existing yield data for 325 hybrid combinations, a hybrid vigor prediction model was established: The model achieved a prediction accuracy of 85.3% in a 10-fold cross-validation, significantly exceeding the approximately 70% accuracy of traditional prediction methods based on parental phenotypes. Based on this model, 50 hybrid combinations with high heterosis potential were recommended, 12 of which demonstrated yields exceeding the parental average by more than 45% in field trials, demonstrating the practical application of the method in breeding.

[0085] Through this embodiment, the application of the method of the present invention in rice breeding significantly improves the quality of genotype data, reduces the number of redundant markers, accurately estimates the relationship between varieties and accurately predicts hybrid vigor, provides strong technical support for rice molecular design breeding, and can be extended to the molecular breeding process of other crop varieties.

[0086] The specific implementation methods described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific implementation method of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A high-throughput genotyping intelligent analysis method, characterized in that: The method comprises the following steps: Step S1, collecting high-throughput genotyping data, standardizing the sample data format, and performing quality assessment on the sample data and removing batch effects to obtain preprocessed genotyping data; The method for quality assessment of sample data includes: calculating sample genotype coverage, calculating the conversion / transversion ratio in the genotype data, and calculating the genotype error rate; screening samples based on the quality assessment results, eliminating samples with genotype coverage of samples; removing transition / transversion ratios Abnormal deviation from the population A sample of and are the expectation and variance of the population transition / transversion ratio, respectively; the genotyping error rate is excluded Sequencing batch data; Step S2, using the preprocessed genotype data, calculate the genetic similarity matrix between samples and perform linkage disequilibrium analysis to identify tag SNPs as molecular markers; Step S3, estimating population kinship and predicting hybrid vigor based on the obtained molecular markers.

2. The high-throughput genotyping intelligent analysis method according to claim 1, characterized in that: The standardized sample data format is specifically: The genotype data obtained from different sequencing platforms are uniformly converted into a standard format; the genotype coding is Numerical representation, where 0 indicates homozygous major allele AA, 1 indicates heterozygous genotype Aa, and 2 indicates homozygous minor allele aa; The missing data were filled using the multiple imputation method based on the local linkage disequilibrium LD structure. The calculation formula is: ,in Representation sample At the site The genotype value to be imputed on Representation and location At the LD threshold The set of all non-deleted sites within the range, Indicates site and site The weight between , For the site and site The LD coefficient between them.

3. The high-throughput genotyping intelligent analysis method according to claim 2, characterized in that: Sample genotype coverage The calculation formula is: ,in, For samples The number of effective genotype sites in is the total number of sites; Transition / transversion ratios in genotype data The calculation formula is: ,in, is the number of conversion mutations; the conversion mutations include: A G or C T, is the number of transversion mutations, including: A C, A T, G C or G T; The left and right represent the genes before and after mutation, respectively. Genotyping error rate The calculation formula is: ,in, is the number of repeated sequencing samples, For repeated samples The number of inconsistent sites, For repeated samples The total number of sites in the CCP.

4. The high-throughput genotyping intelligent analysis method according to claim 1, characterized in that: The method for removing batch effect includes: Constructing the position effect index ,in Indicates batch midpoint The average genotype of and Sites The mean and standard deviation across all batches, For batch Sequencing time, in days, and are the mean and standard deviation of sequencing time for all batches, respectively; right The sites are batch corrected, and the correction formula is: ,in For the site The batch effect coefficient is estimated by minimizing the corrected between-batch variance: , Indicates site A vector of genotype values ​​across all batches.

5. The high-throughput genotyping intelligent analysis method according to claim 4, characterized in that: The genetic similarity matrix between the samples is calculated and comprises the following steps: Based on the preprocessed genotype data, for any two samples and , calculate the sample and samples At the site Correlation on : When the genotypes of two samples are exactly the same , when they share an allele , completely different at the same time ; Calculate genetic similarity weighted by allele frequencies: ; in, , is the site The weight of For the site minor allele frequency.

6. The high-throughput genotyping intelligent analysis method according to claim 5, characterized in that: The linkage disequilibrium analysis is specifically as follows: using the pre-processed genotype data, the linkage disequilibrium coefficient between any two SNP sites is calculated. : ,in, , is the frequency of haplotype AB, and are the frequencies of the major alleles at sites A and B, respectively, and are the frequencies of the minor alleles at site A and site B, respectively; Construct a linkage disequilibrium decay model: ,in, is the effective population size, d is the genetic distance between two loci, in Morgans; Linkage disequilibrium half-distance Evaluate the extent of linkage disequilibrium in a population, i.e. When , the linkage disequilibrium half-decline distance is: .

7. The high-throughput genotyping intelligent analysis method according to claim 6, characterized in that: The identification tag SNP is used as a molecular marker specifically to construct a SNP site network based on linkage disequilibrium analysis. , where the vertex set V represents all SNP sites, and the edge set E represents the linkage disequilibrium relationship. Greater than threshold When , there is an edge between these two sites; Calculate the degree centrality of each SNP site: ,in is the characteristic function; For each connected component, the SNP site with the highest degree centrality is selected as the tag SNP of the connected component; If the degree centrality of multiple sites is the same, the site with the minor allele frequency closest to 0.5 is selected as the tag SNP; Calculate the capture rate of tag SNPs: , where T is the tag SNP set, is the total number of SNP sites.

8. The high-throughput genotyping intelligent analysis method according to claim 1, characterized in that: The method of estimating population kinship based on the obtained molecular markers includes: using the identified tag SNP as a molecular marker to calculate the common ancestor coefficient matrix between samples : ,in is the number of tag SNPs, Representation sample In the The genotype value on the tag SNP, For the The minor allele frequency of each tag SNP; Computing the kinship matrix for population structure correction ,in is the group structure matrix, obtained through the group structure analysis software STRUCTURE, is the diagonal matrix of differentiation coefficients among different subpopulations; Construct a phylogenetic tree based on the kinship matrix to construct an unrooted tree, and the branch length of the tree reflects the genetic distance between samples: .

9. The high-throughput genotyping intelligent analysis method according to claim 8, characterized in that: The method for predicting hybrid vigor is specifically as follows: calculating the expected hybrid F1 heterozygosity of the parental genotypes : ,in, and Respectively represent parents and parents At the site Genotype values ​​on ; Calculate the relationship index between the genetic distance between parents and hybrid vigor: ,in For parents and The genetic distance between Constructing a hybrid vigor prediction model: ,in is the predicted hybrid vigor value, , , and is the regression coefficient, which is estimated by the least squares method based on the known hybrid combination phenotypic data and the corresponding molecular marker data.

Images