Quantitative Trait Multilocus Oscillation Search Genome-Wide Association Analysis System and Method

Abstract

This invention belongs to the field of bioinformatics, specifically relating to a quantitative trait multilocus oscillatory search genome-wide association analysis (GMO) system and method. This invention provides a method for performing GMO analysis using a quantitative trait multilocus oscillatory search genome-wide association analysis system. The method aims to identify SNP loci (PseQTNs) associated with specific phenotypes, and includes the following steps: S1 Data preparation: preparing SNP matrix files and phenotype matrix files; S2 Data preprocessing: removing duplicates from the entered SNP matrix and standardizing it column-wise, and standardizing the phenotype matrix column-wise; S3 Association analysis: using an oscillatory search strategy, screening out phenotype-associated SNP loci (PseQTNs) from the SNP matrix. This invention has significant and practical application value in the mining of genes related to important traits in plants and animals, and in the mining of genes related to complex human diseases.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics, specifically relating to a system and method for quantitative trait multi-site oscillatory search genome-wide association analysis. Background Technology

[0002] GWAS serves as a bridge between genotype and phenotype. With the development of phenomics and high-throughput genotyping technologies, GWAS has been widely applied in fields such as complex human diseases and animal genetics and breeding. Since the first publication of GWAS for sea beet in 2001, more than 50,000 SNPs (p < 5 × 10⁻⁶) have been discovered to date. -8 It is associated with many common diseases or traits, but only a very small number have been functionally validated and it has been subject to doubts about its "lack of heritability". [1] GWAS is essentially feature selection in machine learning: for a dataset (Y... i SNP ij Let i = 1, 2, ..., n; n be the number of samples. j = 1, 2, ..., m; m be the number of SNPs. For simulated data, SNPs = {SNP1, SNP2, ..., SNP...} m}, where QTNs is set to {QTN1, QTN2, ..., QTN} m '}, GWAS finds PseQTNs = {PseQTN1, PseQTN2, ..., PseQTN} m”}, m'≤m, m”≤m, QTN j ∈SNPs, PseQTN j ∈SNPs. The goal of GWAS is to make PseQTNs as similar as possible to QTNs. For real-world data where QTNs are unknown, the goal of GWAS is to ensure that the model built based on PseQTNs has the highest independent prediction accuracy and literature / experimental support for PseQTNs mapping gene results. Assuming m = 2 million SNPs and QTNs = 20, without considering epistatic effects, the difficulty of GWAS is to find exactly 20 QTNs from 2 million SNPs, with a probability of 1 / 10 for complete accuracy. 100 This is akin to finding a needle in a haystack; it highlights that among the three pillars of data, computing power, and algorithms, "algorithm is the core, algorithm is the soul, and algorithm is the key."

[0003] The quantitative trait GWAS has gone through four stages of development.

[0004] The first stage uses a simple model. R- or F-tests are used to examine the correlation between individual SNPs and phenotypes, excluding the population structure Q matrix and kinship (K matrix). A significance threshold of 0.05 is commonly used in statistical tests, meaning researchers must accept a false positive rate as high as 5%. This is an acceptable risk for single-item testing. However, when the number of SNPs reaches millions or even tens of millions, the probability inflation caused by multiple testing leads to a high false positive rate. Common multiple testing correction methods include Bonferroni correction, false discovery rate (FDR), and permutation tests. Bonferroni correction is too strict, resulting in an underestimation of statistical potential. The assumption of independence for each SNP, based on Bonferroni correction and FDR, does not hold in GWAS because SNPs in the least significant (LD) region of the genome will produce similar test statistics, leading to false positives. Permutation tests are considered the gold standard for determining the significance threshold, but they are computationally expensive and practically impractical.

[0005] The second stage is the general linear model (GLM). Age, sex, and population structure (Q) are added as fixed-effect covariates. Population stratification refers to the systematic differences in allele frequencies among subgroups, which can cause spurious associations between genotypes and phenotypes. Population structure is usually represented by principal components (PCs) calculated from the proportion of samples from subgroups or all genotypes, and is called the Q(PCs) matrix. PCs-based correction is currently the most commonly used method to control population stratification, but it can also lead to false negatives.

[0006] The third stage is the mixed linear model (MLM). Q and kinship (K matrix) are added as covariates, where Q represents fixed effects and K represents random effects. The random effects defined by K require estimation of parameters using maximum likelihood or constrained maximum likelihood methods, which is computationally intensive. A series of improvements include EMMA (2008). [2] EMMAX, 2010 [3] P3D, 2010 [4] CMLM, 2010 [4] FaST-LMM, 2011 [5] GEMMA, 2012 [6] GRAMMAR-Gamma, 2012 [7] SUPER, 2014 [8] ; BOLT-LMM, 2015 [9] FastGWA, 2019

[10] These improvements mainly focus on model solution optimization and K construction optimization, and their statistical potential is similar to that of conventional MLM.

[11] .

[0007] The first three stages all use a unit-point model. Unit-point models cannot escape the quagmire of multiple testing and significance thresholds. If the significance threshold is too low, due to multiple testing, a large number of false positive sites will inevitably be selected; if the significance threshold is too high, some true micro-effect QTNs will inevitably be missed, reducing statistical potential.

[0008] The fourth stage is the multi-site model (MuLM). Compared to the single-site model, the multi-site model significantly reduces false positives. There are currently three multi-site models. The first is MLMM.

[12] It is still based on MLM; forward selection continuously adds QTN as a fixed effect to the model, compressing the random effect to 0; backward selection removes redundant pseudo-QTN. The second is FarmCPU.

[13] It is still based on MLM; it uses constrained maximum likelihood optimization instead of the K-explained variance near zero criterion; it alternates between fixed and random effects, using independent random effects models to screen loci and independent fixed effects models to validate loci, alternating between the two until no new candidate markers enter the model; because fixed and random effects do not appear simultaneously in the model, confounding effects are avoided, and false positives and false negatives are controlled simultaneously. The third is BLINK.

[14] It no longer uses MLM; it uses Bayesian information criteria to replace the constrained maximum likelihood estimation in random effects, and uses LD information to select candidate sites; it is faster than MLMM and FarmCPU.

[0009] The shortcomings of existing GWAS schemes, represented by Blink, include: 1) Insufficiently low false positive rate and statistical potential; the Matthews correlation coefficient (MCC) still has room for improvement. In short, there is no best, only better. 2) The accuracy (Top1 Acc%) of the first selected PseQTN is relatively low. Generally, the first PseQTN represents a gene with a strong and important effect, and is the first choice for functional validation. However, in random scenarios with heritability of 0.8, 0.5, and 0.2, Blink's Top1 Acc is only 89%, 76%, and 60%, respectively. 3) In complex scenarios such as TopKSNPs not containing QTNs or tight linkage between QTNs, Blink's MCC and Top1 Acc decrease significantly. Due to the influence of LD and low heritability, the top-ranking SNPS in simple models are not necessarily QTNs. For example, with a heritability of 0.2, R... 2 The probability that the Top 100 SNPs do not contain QTN is 7.80%; at this time, Blink's MCC and Top 1 Acc are as low as 0.0338 and 5%, respectively. 4) There is still room for improvement in computation time complexity on large datasets.

[0010] The reasons for the above shortcomings are multifaceted: 1) GWAS involves multiple disciplines such as bioinformatics, machine learning, and quantitative genetics, requiring a vast amount of background knowledge and a high barrier to entry; it requires comprehensive systemic consideration and overall improvement and optimization, balancing accuracy and computational efficiency; local improvements are often only suitable for specific scenarios and lack general applicability. Even this invention has not yet addressed the superordinate effect. 2) The influence of inertia. Due to computational efficiency constraints, GWAS before 2012 was a unit-point model; in the unit-point model, the introduction of Q and K is considered a milestone. [5] This has led to subsequent improvements consistently using Q and K, focusing solely on model solution optimization and K construction optimization. To date, the MLM model remains the mainstream. Even among the only three multi-site models, MLMM and FarmCPU are still based on MLM and initially included Q and K; it wasn't until 2019 that BLINK abandoned K for the first time (and currently, the only one). 3) It requires years of persistent effort (this invention took nearly four years to reach its basic form), the horizontal transfer of knowledge from multiple disciplines, close teamwork, keen awareness of feedback from numerous failed attempts and partial improvements to clarify the direction of future exploration, and the courage and confidence to abandon Q and K, which predecessors considered milestones.

[0011] References (also listed here for subsequent citations):

[0012] [1]Tam V,Patel N,Turcotte M,et al.Benefits and limitations of genome-wide association studies[J].Nature Reviews Genetics,2019,20(8):467-484.

[0013] [2]Kang HM,Zaitlen NA,Wade CM,et al.Efficient control of population structure in model organism association mapping[J].Genetics,2008,178(3):1709-1723.

[0014] [3]Kang H M,Sul J H,Service S K,et al.Variance component model toaccount for sample structure in genome-wide association studies[J].Naturegenetics,2010,42(4):348-354.

[0015] [4]Zhang Z,Ersoz E,Lai C Q,et al.Mixed linear model approach adaptedfor genome-wide association studies[J].Nature genetics,2010,42(4):355-360.

[0016] [5]Lippert C,Listgarten J,Liu Y,et al.FaST linear mixed models forgenome-wide association studies[J].Nature methods,2011,8(10):833-835.

[0017] [6]Zhou X,Stephens M.Genome-wide efficient mixed-model analysis forassociation studies[J].Nature genetics,2012,44(7):821-824.

[0018] [7]Svishcheva G R,Axenovich T I,Belonogova N M,et al.Rapid variancecomponents–based method for whole-genome association analysis[J].Naturegenetics,2012,44(10):1166-1170.

[0019] [8]Wang Q,Tian F,Pan Y,et al.A SUPER powerful method for genome wideassociation study[J].PloS one,2014,9(9):e107684.

[0020] [9] Loh PR, Tucker G, Bulik-Sullivan BK, et al. Efficient Bayesianmixed-model analysis increases association power in large cohorts[J]. Naturegenetics, 2015, 47(3): 284-290.

[0021]

[10] Jiang L, Zheng Z, Qi T, et al. A resource-efficient tool for mixedmodel association analysis of large-scaledata[J]. Nature genetics, 2019, 51(12):1749-1755.

[0022]

[11] Tan Lizhi, Zhao Yiqiang. Principles, optimization and application of mixture models in genome-wide association analysis [J]. Chinese Agricultural Science, 2023, 56(9): 1617-1632.

[0023]

[12] Segura V, Vilhjálmsson BJ, Platt A, et al. An efficient multi-locusmixed-model approach for genome-wide association studies in structured populations[J]. Nature genetics, 2012, 44(7):825-830.

[0024]

[13] Liu

[0025]

[14] Huang M,Liu X,Zhou Y,et al.BLINK:a package for the next level ofgenome-wide association studies with both individuals and markers in themillions[J].Gigascience,2019,8(2):giy154.

[0026]

[15] Neves H H R,Carvalheiro R,Queiroz S A.A comparison of statisticalmethods for genomic selection in a mice population[J].BMC genetics,2012,13:1-17.

[0027]

[16] Lozada D N,Ward B P,Carter A H.Gains through selection for grainyield in a winter wheat breeding program[J].PLoS One,2020,15(4):e0221603.

[0028]

[17] Zhao K,Tung C W,Eizenga G C,et al.Genome-wide association mappingreveals a rich genetic architecture of complex traits in Oryza sativa[J].Nature communications,2011,2(1):467.

[0029]

[18] Atwell S,Huang Y S,Vilhjálmsson B J,et al.Genome-wide associationstudy of 107phenotypes in Arabidopsis thaliana inbred lines[J].Nature,2010,465(7298):627-631.

[0030]

[19] Yin H,Zhou C,Shi S,et al.Weighted single-step genome-wideassociation study of semen traits in Holstein bulls of China[J].Frontiers ingenetics,2019,10:1053.

[0031]

[20] Romay M C,Millard M J,Glaubitz J C,et al.Comprehensive genotypingof the USA national maize inbred seed bank[J].Genome biology,2013,14:1-18.

[0032]

[21] Danilevskaya O N,Meng X,Selinger D A,et al.Involvement of theMADS-box gene ZMM4 in floral induction and inflorescence development in maize[J].Plant physiology,2008,147(4):2054-2069.

[0033]

[22] Buckler E S,Holland J B,Bradbury P J,et al.The geneticarchitecture of maize flowering time[J].Science,2009,325(5941):714-718.

[0034]

[23] Kumimoto R W,Zhang Y,Siefers N,et al.NF–YC3,NF–YC4 and NF–YC9 arerequired for CONSTANS-mediated,photoperiod-dependent flowering in Arabidopsisthaliana[J].The Plant Journal,2010,63(3):379-391.

[0035]

[24] Wang X,Wu L,Zhang S,et al.Robust expression and association ofZmCCA1 with circadian rhythms in maize[J].Plant Cell Reports,2011,30:1261-1272.

[0036]

[25] Jin M,Liu X,Jia W,et al.ZmCOL3,a CCT gene represses flowering inmaize by interfering with the circadian clock and activating expression ofZmCCT[J].Journal of integrative plant biology,2018,60(6):465-480.

[0037]

[26] Huang C,Sun H,Xu D,et al.ZmCCT9 enhances maize adaptation tohigher latitudes[J].Proceedings of the National Academy of Sciences,2018,115(2):E334-E341.

[0038]

[27] Liang Y,Liu Q,Wang X,et al.Zm MADS 69functions as a floweringactivator through the ZmRap2.7-ZCN 8regulatory module and contributes tomaize flowering time adaptation[J].New Phytologist,2019,221(4):2335-2347

[0039]

[28] Shi Y,Zhao X,Guo S,et al.ZmCCA1a on chromosome 10of maize delaysflowering of Arabidopsis thaliana[J].Frontiers in Plant Science,2020,11:78.

[0040]

[29] Li Q,Wu G,Zhao Y,et al.CRISPR / Cas9-mediated knockout and overexpression studies reveal a role of maize phytochrome C in regulating flowering time and plant height[J].Plant biotechnology journal,2020,18(12):2520-2532.

[0041]

[30] Su H,Liang J,Abou-Elwafa S F,et al.ZmCCT regulates photoperiod-dependent flowering and response to stresses in maize[J].BMC plant biology,2021,21:1-15.

[0042]

[31] Li Z,Gao F,Liu Y,et al.ZmGI2 regulates flowering time through multiple flower development pathways in maize[J].Plant Science,2023,332:111701.

[0043]

[32] Yang J,Wei H,Hou M,et al.ZmSPL13 and ZmSPL29 act together to promote vegetative and reproductive transition in maize[J].New Phytologist,2023,239(4):1505-1520.

[0044]

[33] Qi X,Wang Y,Zhang X,et al.Cloning,Characterization and Transformation of Methyltransferase 2a Gene(Zmet2a)in Maize(Zea mays L.)[J].Phyton(0031-9457),2024,93(7). Summary of the Invention

[0045] The MOST-GWAS proposed in this invention belongs to the fourth-stage multi-site MuLM model. It is the first to propose a forward feature selection oscillation search strategy and is completed using a new termination criterion.

[0046] This invention provides a method for genome-wide association analysis using quantitative trait multilocus oscillation search (PseQTNs), which aims to identify SNPs associated with specific phenotypes and includes the following steps:

[0047] S1 Data Preparation: Prepare SNP matrix files and phenotype matrix files; the SNP matrix files record genotype data, and the phenotype matrix files record quantitative trait data.

[0048] S2 data preprocessing: The entered SNP matrix is ​​deduplicated, standardized by column, and principal component analysis is performed. The phenotypic matrix is ​​standardized by column.

[0049] S3 association analysis: Using an oscillating search strategy, PseQTNs, which are phenotype-related SNP sites, are selected from the SNP matrix.

[0050] Specifically, the SNP matrix file supports .bed, .raw, and .csv files from the plink software, where the first line of the .csv file is the SNPID, and the remaining lines are the genotypes of each sample;

[0051] The phenotypic matrix file supports .csv files. The first row is the phenotypic name, the first column is the sample name, and the subsequent columns are the continuous phenotypic values. Multiple pictographic data can be analyzed sequentially.

[0052] In a specific implementation, the method for S3 correlation analysis is as follows:

[0053] 1) Calculate the R-values ​​of the SNP matrix S and a given phenotype Y in the phenotype matrix. 2 (Y, S) i And sort them in descending order by Rank, and denote the first SNP as . * S1; with |r( * S1,S i If | < 0.7, filter the Rank and record the second SNP as . * S2; and so on, until 10 relatively independent initial SNPs are obtained, denoted as S2; * S k k = 1, 2, ..., 10, serving as the starting points for 10 chains;

[0054] 2) Perform principal component analysis on the initial SNPs, and select the first few principal components (e.g., 2 or 3) according to the preset parameters, and denote them as PCs;

[0055] 3) For each chain, use PCs+ * Sk As covariates, for Y and S i Find the partial correlation and denote the value of S corresponding to the largest absolute value of the partial correlation. i for k S1; referred to as chain extension step; then PCs+ * S k + k S1 is a covariate, relating Y and S. i Find the partial correlation and denote the value of S corresponding to the largest absolute value of the partial correlation. i for k S2; and so on, extending for a total of 10 steps;

[0056] 4) For 10 chains, each chain contains PCs and 11 SNPs, calculate the multiple correlation. The chain with the largest absolute value of the multiple correlation is the optimal chain.

[0057] 5) For the optimal chain, cut off * S+S1, PCs+S2+S3+…+S 10 As covariates, for Y and S i Finding partial correlations is called chain repair; let S be the value of the partial correlation with the largest absolute value. i For Q1; obtained the first PseQTN;

[0058] 6) Using PCs+Q1 as covariates, perform chain extension for 5 steps as in step 3), and chain repair as in step 5 to obtain Q2.

[0059] 7) Using PCs+Q1+Q2 as covariates, perform chain extension for 5 steps as in step 3), and chain repair as in step 5 to obtain Q3;

[0060] 8) Continue in this manner until the chain extension termination criterion is met: ΔR 2 <TV best / n, TV best =-20×h 2 +21;

[0061] 9) Let {Q1, Q2, ..., Q} be the integers of ... m Let PseQTNs be the multivariate inequality (Ri) of Y and PseQTNs + PCs. 2 PCs+QTNs The multiple correlation between Y and QTNs R 2 QTNs ;

[0062] 10) Repeat steps 1)-9) until all columns of the phenotypic matrix have been calculated, and output the results.

[0063] This invention provides a system for genome-wide association analysis using quantitative trait multi-site oscillation search, comprising the following modules: a data entry module, a data preprocessing module, an association analysis module, and a result output module;

[0064] The data preprocessing module is used to perform step S2 of the method; the correlation analysis module is used to perform step S3 of the method.

[0065] Specifically, the data entry module is used to output SNP matrix files and phenotype matrix files. More specifically, the SNP matrix files support .bed, .raw, and .csv files from the Plink software, where the first row of the .csv file is the SNPID, and the remaining rows are the genotypes of each sample; the phenotype matrix files support .csv files, where the first row is the phenotype name, the first column is the sample name, and subsequent columns are continuous phenotype values, supporting sequential analysis of multiple phenotypic tables.

[0066] In a specific implementation, the results output module outputs the results as follows: corresponding to the number of phenotypes in the input phenotype matrix, each phenotype is output as an Excel file containing two sheets. The first sheet contains SNPID, P-Value, and estimated effect size (only for PseQTNs), and the second sheet contains the heritability explained by PseQTNs.

[0067] The present invention also provides an apparatus containing the system, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, the processor executing the computer program, the computer program being encoded to implement the system.

[0068] Optionally, it also includes a results display device, such as a mobile phone, computer, or remote results display via the internet.

[0069] The method and system of this invention have significant advantages, including high accuracy (a high Matthews correlation coefficient (MCC) represents a balance between high statistical potential and low false positives), strong robustness (stable and effective in complex scenarios such as strong LD among QTNs and TopK SNPs without QTNs), and high computational efficiency (suitable for large-sample, high-density SNP data). Specifically:

[0070] 1) High precision and strong robustness. A pioneering forward feature selection oscillating search strategy is employed, including: a) A "borrowing sheep" strategy for the population structure Q matrix. Borrowing Q controls false positives, while returning Q increases statistical potential. Introducing Q as a covariate helps select QTNs, but once PseQTNs have been selected and the heritability explained by Q+PseQTNs is close to that of PseQTNs, Q should be phased out. b) Tree search. Instead of hastily identifying a SNP as a PseQTN, a few steps are taken and searched (chain extension), followed by a reassessment of the reliability of the previously selected SNPs (chain repair). c) Chain repair correction. Adding QTNs or strong LD-SNPs of QTNs as covariates can improve the accuracy of selecting other higher-effect QTNs. d) Under multiple covariate combinations, only SNPs are sorted rather than deleted, avoiding the quagmire of multiple testing and significance thresholds, ensuring that weak-effect QTNs have a chance to be selected at each step of chain extension, effectively improving statistical potential. e) Semi-correlation with LD candidate SNPs and partial correlation to determine PseQTNs. This avoids missing selections when there is strong LD among QTNs and ensures computational efficiency. f) Nested linkage with imbalanced LD candidate multiple initial points increases the diversity of initial points, making the model robust even in complex scenarios where TopKSNPs do not contain QTNs.

[0071] 2) Computationally efficient. a) K-free. Random effects defined by K require maximum likelihood or constrained maximum likelihood estimation of parameters, which is computationally intensive. At least for simulated data, the genetic model of phenotype Y itself does not contain K and Q, so the inclusion of K and Q in the final genetic model is not inevitable. Q (PCs) and K originate from all SNPs, and confounding factors between population structure Q, kinship K, and quantitative trait nucleotides (QTNs) can lead to false negatives. Similar attempts to remove K have been made: MLMM compresses K to 0, FarmCPU uses constrained maximum likelihood optimization to replace the standard of near-zero variance explained by K, and BLINK abandons K and no longer uses the MLM model. b) Uses incomplete tree search (width W = 10, depth D = 10, 10 chains). Complete tree search (width W = 10, depth D = 10, 10 chains). 10 (a chain) is very time-consuming. c) Perform semi-correlation and LD on all SNPs using matrix vectorization to determine a small number of candidate SNPs; perform partial correlation only on a small number of candidate SNPs to determine PseQTNs.

[0072] 3) New termination criteria. Multiple linear regression multiple correlation coefficient R0 2 Critical value R 2 (ɑ,ν,M) It is related to the degrees of freedom ν, the number of variables m, and the significance level α. Clearly, R... 2 The critical value is affected by both the number of samples and the number of variables. For forward feature selection, this invention, through extensive numerical experiments, has discovered a critical value ΔR.2 =R 2 (ɑ,ν,m+1) -R 2 (ɑ,ν,m) It is only affected by the sample size and α. Let ΔR 2 =

[0073] If TV / n, then TV is only related to ɑ. In high h... 2 At this time, more PseQTNs can be introduced to ensure high statistical potential, α should be large and TV should be small; at low h 2 At this time, it is necessary to more strictly control the introduction of PseQTNs to control false positives, with α being smaller and TV larger. The objective function is Max{MCC}, and h... 2 =1 hour TV best =1 is a constraint; the search yields different h values. 2 The empirical formula for the optimal value of TV is: TV best =-20×h 2 +21; h 2 When the number of characteristics increases from m to m+1 Characterization. When ΔR 2 <TV best The import process terminates when / n is encountered. Attached Figure Description

[0074] Figure 1 Manhattan plots and schematic diagrams showing changes in QTN ranking under different covariates. Note: 11,089 SNPs, 1,768 samples. Simulated phenotypic settings: QTNs = 10, h 2 =0.8, additive effect, effect size normally distributed. QTN index represents the effect, QTN1 has the largest effect. PCs = 3.

[0075] Figure 2 The MOST-GWAS flowchart of this invention. Detailed Implementation

[0076] The present invention will be described below through specific embodiments in order to better understand the present invention, but this does not constitute a limitation on the present invention.

[0077] 1. Related abbreviations, full names, and explanations:

[0078] Absolute Acc%: Absolute accuracy. True QTNs / Settled QTNs × 100%.

[0079] Blink: Bayesian information and linkage disequilibrium iteratively nested keyway. A multi-site GWAS algorithm.

[0080] FarmCPU: Fixed and random model circulating probability unification. A multi-site GWAS algorithm.

[0081] FDR: False discovery rate. A commonly used evaluation metric for the GWAS algorithm.

[0082] GLM: General linear model. A type of genetic model.

[0083] GP / GS: Genomic prediction / selection. It involves modeling a training population with known genotypes and phenotypes to predict the phenotype of a test population whose genotypes are only known.

[0084] GWAS: Genome-wide association studies. A statistical analysis method for identifying phenotype-associated QTNs or genes from a large number of genotypes (SNPs).

[0085] h 2 Heritability.

[0086] K: Kinship, a matrix of kinship relationships between individuals.

[0087] LD: Linkage disequilibrium. It refers to the non-random association between alleles at different loci, where the probability of two genes being inherited simultaneously in a population is greater than the probability of random inheritance.

[0088] MCC: Matthews' correlation coefficient. A commonly used evaluation metric for binary classification algorithms in machine learning.

[0089] MLM: Mixed linear model. A linear model that includes both fixed and random effects.

[0090] MLMM: Multi-locus mixed model, multi-locus mixed linear model, multi-locus GWAS algorithm.

[0091] MOST-GWAS: Multi-locus oscillatory search towards genome-wide association studies for quantitative traits, a multi-locus oscillatory search GWAS-regression algorithm for quantitative traits proposed in this patent application.

[0092] MuLM: Multi-locus model. The model contains multiple QTNs / SNPs.

[0093] PC: Principal Component Analysis. A linear dimensionality reduction algorithm.

[0094] PseQTNs: PseudoQTNs, candidate QTNs found by GWAS.

[0095] Q: Population structure, the population structure matrix, is usually represented by the first few principal components (PCs).

[0096] QTL: Quantitative trait loci. The location of genes controlling quantitative traits in the genome.

[0097] QTN: Quantitative trait nucleotide. The SNP within a QTL that truly contributes to quantitative trait variation.

[0098] R 2 Squared Pearson correlation coefficient: The square of the correlation coefficient. It measures the degree of linear correlation between two continuous variables and is a commonly used evaluation metric in regression algorithms in machine learning.

[0099] Relative Acc%: Relative accuracy. True QTNs / Selected QTNs × 100%.

[0100] SNP: Single nucleotide polymorphism. A DNA sequence polymorphism caused by a single nucleotide variation at the genomic level, occurring at a frequency greater than 1%.

[0101] TV: Threshold value.

[0102] Example

[0103] 1. Mathematical description of the problem

[0104] Existing dataset (Y) i ,S ij Let i = 1, 2, ..., n; n be the number of samples. j = 1, 2, ..., m; m be the number of SNPs. Y is the phenotype. Assume that Y is formed by only m' linear combinations of SNPs (QTNs) and that noise (heritability h) exists. 2 <1), How to find the m' SNPs (QTNs) as accurately as possible from m SNPs?

[0105] 2. The starting point of this invention

[0106] SNPs are classified into four categories: The first category is TP-SNPs (True Positives), also known as QTNs. Note that TP-SNPs may exhibit varying degrees of redundancy due to linkage disequilibrium (LD). The second category is false-positive PS-SNPs caused by population stratification (PS). Population stratification results from systematic differences in allele frequencies between subpopulations within a population due to non-random mating between individuals, leading to a large number of false-positive loci in GWAS. Using Q(PCs) as a covariate is currently the most commonly used method to control for population stratification. The third category is false-positive LD-SNPs caused by linkage disequilibrium (LD). Linkage disequilibrium (LD) is a non-random association between alleles at different loci. A single TP-SNP often has multiple LD-SNPs, which is the main cause of false positives. The fourth category is irrelevant marker Ir-SNPs, which are not phenotypic. Due to low heritability, small sample size, and high density of SNPs, some Ir-SNPs may pass hypothesis testing.

[0107] This invention is based on three facts frequently observed in practice (such as...) Figure 1 As shown):

[0108] 1) Top-ranked SNPs in simple models S It doesn't necessarily have to be QTNs. Figure 1 Among the 10 QTNs in A, QTN9, which has the strongest effect, ranks 5th, while the 1st is the strong LD-labeled SNP of QTN9. 10118 .

[0109] 2) Including QTNs as a covariate can improve the ranking of other QTNs with higher effects. For example, using PCs+QTN9 as a covariate, from... Figure 1 B to Figure 1In the C category, QTN1's ranking rose from 44 to 19, and QTN2's ranking rose from 19 to 1.

[0110] 3) Including strong LD-SNPs of QTNs as covariates can also improve the ranking of some higher-effect QTNs. For example, from Figure 1 B to Figure 1 In the D category, the ranking of QTN1 also rose from 44 to 19, and the ranking of QTN2 also rose from 19 to 1 (genotype data are from the literature).

[15] ).

[0111] Therefore, a possible subset of PseQTNs (chain extension) can be obtained through forward search, such as PCs + SNPs. 10118 +QTN2+SNP 3657 +QTN1; then with {PCs, SNP 3657 QTN1} is the covariate used to correct for the first selected SNP. 10118 (Chain repair) At this point, QTN9 is highly likely to be selected correctly. This strategy is collectively known as oscillation search.

[0112] 3. Algorithm steps of this invention

[0113] The specific process of the MOST-GWAS algorithm of this invention is as follows: Figure 2 As shown, the specific steps are as follows:

[0114] 1) Matrix operations to calculate R 2 (Y, S) i And sort them in descending order by Rank, and denote the first SNP as . * S1; with |r( * S1,S i If | < 0.7, filter the Rank and record the second SNP as . * S2; and so on, until 10 relatively independent initial SNPs are obtained, denoted as S2; * S k k = 1, 2, ..., 10, serving as the starting points for 10 chains.

[0115] 2) Perform principal component analysis on SNPs and select a number of principal components (e.g., 2 or 3, default is 2) to represent the Q matrix, denoted as PCs.

[0116] 3) For each chain, use PCs+ * S k As covariates, for Y and S i Find the partial correlation and denote the value of S corresponding to the largest absolute value of the partial correlation. i for k S1; this is called a chain extension step. Then, PCs+ * Sk + k S1 is a covariate, relating Y and S. i Find the partial correlation and denote the value of S corresponding to the largest absolute value of the partial correlation. i for k S2; and so on, for a total of 10 steps.

[0117] 4) For 10 chains, each chain contains PCs and 11 SNPs, calculate the multiple correlation. The chain with the largest absolute value of the multiple correlation is the optimal chain.

[0118] 5) For the optimal chain, cut off * S+S1, PCs+S2+S3+…+S 10 As covariates, for Y and S i Finding partial correlations is called chain repair. Let S be the value of the partial correlation with the largest absolute value. i For Q1; obtain the first PseQTN.

[0119] 6) Using PCs+Q1 as covariate, extend the chain in 5 steps and repair the chain to obtain Q2.

[0120] 7) Using PCs+Q1+Q2 as covariates, extend the chain in 5 steps and repair the chain to obtain Q3.

[0121] 8) Continue in this manner until the chain extension termination criterion is met: ΔR 2 <TV best / n, TV best =-20×h 2 +21.

[0122] 9) Let {Q1 + Q2 + ... + Q} m’ Let PseQTNs be the total correlation coefficients (R) between Y and QTNs + PCs. 2 PCs+QTNs The multiple correlation between Y and QTNs R 2 QTNs .

[0123] 10) Output results: The number of phenotypes corresponding to the input phenotype matrix. Each phenotype is output as an Excel file containing two sheets. The first sheet contains SNPID, P-Value, and estimated effect size (only for PseQTNs). The second sheet contains the heritability explained by PseQTNs.

[0124] 4. Model Evaluation Criteria

[0125] Existing standards such as statistical power, Type I error rate, FDR (false positives / total positives), relative accuracy (Relative Acc%), and absolute accuracy (Absolute Acc%) all evaluate from a single perspective and are flawed. For example, if a model finds only one correct QTN (QTs) with a set QTN value of 20, the relative accuracy is very high (100%), but the statistical power and absolute accuracy are very low. Conversely, if a model finds 10,000 QTNs, including 20 true QTNs, the statistical power and absolute accuracy are very high, but the relative accuracy, Type I error rate, and FDR are very high.

[0126] The Matthews correlation coefficient (MCC), commonly used in machine learning, is a good comprehensive evaluation metric for binary classification. The confusion matrix for simulated data with known real QTNs is as follows:

[0127]

[0128] Accuracy is flawed when there is an imbalance between positive and negative samples: for example, with 10,000 SNPs and QTNs=10, even if a model fails to identify any (classifying all as negative), the accuracy is still as high as 99.90%. GWAS suffers from extreme imbalance between positive and negative samples, and MCC takes this imbalance into account.

[0129]

[0130] True Positive (TP): This SNP is a true QTN, and the model also classifies it as a QTN.

[0131] False Negative (FN): This SNP is a true QTN, but the model did not classify it as a QTN.

[0132] False Positive (FP): This SNP is a false QTN, and the model classifies it as QTN.

[0133] True Negative TN: This SNP is a false QTN, and the model does not classify it as a QTN.

[0134] In addition, the first QTN selected by the model often represents a gene with a strong and important effect, and is the first choice for functional validation; therefore, the accuracy of the first QTN (Top1ACC%) is also an important evaluation indicator.

[0135] Therefore, this invention selects MCC and Top1ACC% as the main evaluation indicators.

[0136] Verification Example 1: Simulated data, performance of various models under different heritability in random scenarios.

[0137] Using winter wheat dataset

[16] The number of SNPs was 11089, and the sample size was 1768. The random scenario refers to QTNs being randomly located on chromosomes without clustering. In each simulation, 20 SNPs were randomly selected as QTNs, meaning only 20 SNPs formed a linear combination of phenotype Y, and the additive effect values ​​of these 20 QTNs followed a normal distribution. Heritability h 2 Set three levels: 0.8, 0.5, and 0.2. 2 =0.8 indicates that 20 QTNs explain 80% of the phenotypic variation. The model is required to identify a subset of PseQTNs from 11089 SNPs, compare it with a defined subset of QTNs, and calculate the Matthews correlation coefficient (MCC), the top-1 accuracy (Accel%) of the first PseQTN, and the runtime. The simulation is repeated 100 times, and the results are expressed as the mean. Reference models include the unit-point model SUPER, and three existing multi-site models: MLMM, FarmCPU, and Blink.

[0138] Table 1. Performance of each model under different heritability in random scenarios.

[0139]

[0140]

[0141] Conclusion: Compared to the single-site model SUPER, all four multi-site models effectively controlled false positives. MOST-GWAS had the highest MCC, the highest overall Top1 Acc%, and the fastest speed.

[0142] Verification Example 2: Simulated data, performance of various models under different heritability in complex scenarios without QTN in the Top 100

[0143] SNPs with high ranking in simple models S It may not necessarily contain QTNs. Simulations showed that when heritability was 0.2, 0.5, and 0.8, R... 2 The probabilities of the Top 100 SNPs not containing QTNs were 7.80%, 3.90%, and 3.50%, respectively; R was low in heritability. 2 The absence of QTNs in the Top 100 is not a low-probability event. The first PseQTN selected by the model often represents a gene with a strong and important effect, making it the preferred target for functional validation. In R... 2 In complex scenarios where the Top100 SNPs do not include QTNs, evaluating the performance of each model, especially Top1Acc%, is particularly important.

[0144] Using winter wheat dataset

[16] The simulation process is the same as in Example 1, repeated 3000 times, from which 100 R values ​​are selected. 2 The dataset of Top 100 SNPs without QTNs is used in this example.

[0145] Table 2. Performance of Top 100 models with different genetic forces in complex scenarios without QTN

[0146]

[0147] Conclusion: MOST-GWAS has the highest MCC and Top1 Acc%. Compared with random scenarios (Table 1), MOST-GWAS has a significant advantage over the reference model in complex scenarios where the Top100 SNPs do not contain QTN.

[0148] Verification Example 3: Simulated Data, Performance of Different Models under Different Heredity in a Complex Scenario with Close Linkage of QTN

[0149] Using winter wheat dataset

[16] The simulation process is the same as in Example 1, but 5 pairs of QTNs among the selected 20 QTNs are closely linked, 0.7 <LD<0.95。

[0150] Table 3. Performance of different models under different heritability in complex scenarios with tight linkage of QTN

[0151]

[0152]

[0153] Conclusion: MOST-GWAS has the highest MCC and Top1 Acc%. Compared with random scenarios (Table 1), MOST-GWAS shows a more significant advantage over the reference model in complex scenarios with tightly linked QTN structures.

[0154] Verification Example 4: Simulated Data, Comparison of Multi-Species Performance of Blink and MOST-GWAS in Random Scenarios

[0155] To understand the generalizability of the models, the performance of Blink and MOST-GWAS models was further compared under five different species random scenarios. The winter wheat dataset contains 11,089 SNPs and 1,768 samples.

[16] The rice dataset contains 36,898 SNPs and 413 samples.

[17] (http: / / www.ricediversity.org / data / sets / 44kgwas / ), Arabidopsis thaliana dataset

[18] The mouse dataset contains 214,545 SNPs and 199 germplasm (https: / / archive.gramene.org / diversity / download_data.html).

[15] The Holstein cattle dataset contains 12,226 SNPs and 1,940 samples (1,000 males and 940 females).

[19] There are 42,551 SNPs and 5,024 samples.

[0156] For each species and for each heritability, 20 SNPs are randomly selected as QTNs in each simulation, and the simulation is repeated 100 times. The results are expressed as the mean.

[0157] Table 4. Comparison of multi-species performance of Blink and MOST-GWAS in random scenarios.

[0158]

[0159]

[0160] Conclusion: For five different species, MOST-GWAS showed better stability in MCC and Top1 Acc% than Blink.

[0161] Verification Example 5: Corn Flowering Time GWAS (Real Data)

[0162] The maize dataset contains 681,258 SNPs and 2,279 inbred lines.

[20] (http: / / cbsusrv04.tc.cornell.edu / users / panzea / filegateway.aspx?category=Genotypes), the phenotype is the flowering time, i.e. the number of days to spin silk (https: / / www.panzea.org / phenotypes).

[0163] Table 5. PseQTNs selected by Blink and MOSTR-GWAS

[0164]

[0165] Blink selected 49 PseQTNs, MOST-GWAS selected 42 PseQTNs, and together they selected 6 PseQTNs.

[0166] Table 6. Alignment results of PseQTNs selected by MOST-GWAS with known flowering genes (PCs = within 2,900 Kb)

[0167]

[0168] 42 PseQTNs explained 71.97% of the heritability among 17 reported maize flowering genes. [21-33] Six known flowering genes were identified in the study; the most recent one was 528 bp.

[0169] Table 7. Alignment results of PseQTNs selected by Blink with known flowering genes (PCs = within 2,900 Kb)

[0170]

[0171] 49 PseQTNs explained 73.79% of the heritability among 17 reported maize flowering genes. [21-33] Four known flowering genes were matched in the study, and two QTNs matched the same gene three times; the most recent match was 6500 bp.

[0172] MOST-GWAS matched 6 known genes with 42 PseQTNs, the closest being 528 bp; while Blink matched only 4 known genes with 49 PseQTNs, the closest being 6500 bp. MOST-GWAS is significantly superior to Blink.

Claims

1. A method for genome-wide association analysis using quantitative trait multi-site oscillation search, characterized in that, The method aims to identify SNP sites, or PseQTNs, associated with a specific phenotype, and includes the following steps: S1 Data Preparation: Prepare SNP matrix files and phenotype matrix files; the SNP matrix files record genotype data, and the phenotype matrix files record quantitative trait data. S2 data preprocessing: The entered SNP matrix is ​​deduplicated, standardized by column, and principal component analysis is performed. The phenotypic matrix is ​​standardized by column. S3 association analysis: Using an oscillating search strategy, PseQTNs (SNP sites) related to the phenotype are screened from the SNP matrix; wherein the SNP matrix file supports .bed, .raw, and .csv files from the plink software, where the first line of the .csv file is the SNPID, and the remaining lines are the genotypes of each sample; The phenotypic matrix file supports .csv files. The first row is the phenotypic name, the first column is the sample name, and the subsequent columns are continuous phenotypic values. Multiple table types can be analyzed sequentially. The method for S3 correlation analysis is as follows: 1) Calculate the R²(Y, Si) of the SNP matrix S and a certain phenotype Y in the phenotype matrix, and sort them in descending order. Let the first SNP be denoted as . S1; with |r( Filter the Rank if S1,Si)|<0.7 and denote the second SNP as S2; and so on, until 10 relatively independent initial SNPs are obtained, denoted as *Sk, k=1, 2, ..., 10, which serve as the starting points of the 10 chains; 2) Perform principal component analysis on the initial SNPs, and select the first few principal components according to the preset parameters, denoted as PCs; 3) For each chain, use PCs + Sk is a covariate. The partial correlation between Y and Si is calculated, and the Si corresponding to the largest absolute value of the partial correlation is denoted as kS1; this is called chain extension one step; then PCs + Sk + kS1 are covariates. Calculate the partial correlation between Y and Si, and denote the Si corresponding to the largest absolute value of the partial correlation as kS2. Continue in this manner for a total of 10 steps. 4) For 10 chains, each chain contains PCs and 11 SNPs, calculate the multiple correlation. The chain with the largest absolute value of the multiple correlation is the optimal chain. 5) For the optimal chain, cut off *S +S1, and use PCs +S2 +S3 + … +S10 as covariates to calculate the partial correlation between Y and Si, which is called chain repair; denote the Si corresponding to the one with the largest absolute value of partial correlation as Q1; obtain the first PseQTN; 6) Using PCs + Q1 as covariates, perform chain extension for 5 steps as described in step 3), and chain repair as described in step 5 to obtain Q2; 7) Using PCs + Q1 + Q2 as covariates, perform chain extension for 5 steps as in step 3), and chain repair as in step 5 to obtain Q3; 8) Continue in this manner until the chain extension termination criterion is met: ΔR2 <TVbest / n, , where n is the number of samples; 9) Let {Q1, Q2, ..., Qm} be PseQTNs. Find the multiple correlation R2PCs+QTNs between Y and PseQTNs+PCs, and the multiple correlation R2QTNs between Y and QTNs. 10) Repeat steps 1)-9) until all columns of the phenotypic matrix have been calculated, and output the results.

2. A system for genome-wide association analysis using quantitative trait multi-site oscillation search, characterized in that, It includes the following modules: data entry module, data preprocessing module, correlation analysis module, and result output module; The data preprocessing module is used to perform step S2 of the method as described in claim 1; the correlation analysis module is used to perform step S3 of the method as described in claim 1.

3. The system as described in claim 2, characterized in that, The data entry module is used to input SNP matrix files and phenotypic matrix files.

4. The system as described in claim 3, characterized in that, The SNP matrix file supports .bed, .raw, and .csv files from the plink software. The first line of the .csv file is the SNPID, and the remaining lines are the genotypes of each sample. The phenotypic matrix file supports .csv files. The first row is the phenotypic name, the first column is the sample name, and the subsequent columns are the continuous phenotypic values. Multiple pictographic data can be analyzed sequentially.

5. The system as described in claim 2, characterized in that, The results output module outputs the results as follows: corresponding to the number of phenotypes in the input phenotype matrix, each phenotype is output as an Excel file containing two sheets. The first sheet contains SNPID, P-Value, and estimated effect value, and the second sheet contains the heritability explained by PseQTNs.

6. An apparatus comprising the system of any one of claims 2 to 5, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, the processor executing the computer program, the computer program being encoded to implement the system of any one of claims 2 to 5.

7. The device as described in claim 6, characterized in that, It also includes equipment for displaying results.

8. The device as described in claim 7, characterized in that, The results are displayed via mobile phone, computer, or remote network.

Images