Quantitative character multi-site oscillation search whole genome association analysis system and method

Abstract

The invention belongs to the field of bioinformatics, and particularly relates to a quantitative character multi-site oscillation search whole genome association analysis system and method. The invention provides a method for performing whole genome association analysis by using a quantitative character multi-site oscillation search whole genome association analysis system, the method aims to find out SNP sites, namely PseQTNs, related to specific phenotypes, and the method comprises the following steps: S1, data preparation: preparing an SNP matrix file and a phenotype matrix file; s2, data preprocessing: de-duplicating the input SNP matrix, standardizing the matrix according to columns, and standardizing a phenotype matrix according to columns; and S3, correlation analysis: screening out phenotype-related SNP sites, namely PseQTNs, from the SNP matrix through an oscillation search strategy. The method has important and practical application value in the aspects of animal and plant important character related gene mining, human complex disease related gene mining and the like.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, and in particular relates to a system and method for quantitative trait multi-locus oscillation search and genome-wide association analysis. Background Art

[0002] GWAS is a bridge between genotype and phenotype. With the development of phenomics and high-throughput genotyping technology, GWAS has been widely used in the fields of complex human diseases and animal genetic breeding. Since the first GWAS of sea beet was published in 2001, more than 50,000 SNPs (p < 5 × 10 -8 ) are associated with a variety of common diseases or traits, but only a few have been functionally verified and are questioned for their “missing heritability” [1] GWAS is essentially feature selection in machine learning: i ,SNP ij ), i = 1, 2, ..., n; n is the number of samples. j = 1, 2, ..., m; m is the number of SNPs. For the simulated data, SNPs = {SNP1, SNP2, ..., SNP m}, set QTNs = {QTN1, QTN2, …, QTN m '}, PseQTNs found by GWAS = {PseQTN1, PseQTN2, …, PseQTN m”}, m'≤m, m"≤m, QTN j ∈SNPs,PseQTN j ∈SNPs. The goal of GWAS is to make PseQTNs and QTNs as similar as possible. For real data with unknown QTNs, the goal of GWAS is to make the model built based on PseQTNs have the highest independent prediction accuracy and literature / experimental support for the results of PseQTNs mapping genes. Assuming m = 2 million SNPs, QTNs = 20, and ignoring epistasis effects, the difficulty of GWAS is to find exactly 20 QTNs from 2 million SNPs, and the probability of being completely correct is 1 / 10 100 , which is like looking for a needle in a haystack; it highlights that among the three pillars of data, computing power, and algorithm, "algorithm is the core, algorithm is the soul, and algorithm is the key."

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

[0004] The first stage is a simple model. R or F-tests are used to test the correlation between a single SNP and the phenotype, excluding the population structure Q matrix and the kinship (K matrix). The significance threshold commonly used in statistical tests is 0.05, which means that researchers have to accept a false positive rate of up to 5%. This is an acceptable risk for a single test. However, when the number of SNPs is as high as millions or even tens of millions, the probability inflation caused by multiple tests leads to a high false positive rate. Commonly used multiple test correction methods include Bonferroni correction, false discovery rate FDR and permutation test. Bonferroni correction is too strict, resulting in low statistical power. The assumption that each SNP is independent based on the Bonferroni correction and FDR does not hold in GWAS, because SNPs in LD within the genomic region will produce similar test statistics, resulting in false positives. The permutation test is considered to be the gold standard method for determining the significance threshold, but it is computationally intensive and almost infeasible in practical applications.

[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 between subpopulations, which can cause false associations between genotypes and phenotypes. The principal components PCs calculated from the proportion of samples from subpopulations or all genotypes are usually used to represent the population structure, which 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 is a fixed effect and K is a random effect. The random effect defined by K requires the maximum likelihood method or constrained maximum likelihood method to estimate the parameters, 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 are mainly focused on model solution optimization and K construction optimization, and its statistical potential is similar to that of conventional MLM.

[11] .

[0007] The first three stages are all single-site models. Single-site models cannot escape the quagmire of multiple testing / significance thresholds. If the significance threshold is too low, a large number of false positive sites will inevitably be selected due to multiple testing; if the significance threshold is too high, some real micro-effect QTNs will inevitably be missed, reducing statistical power.

[0008] The fourth stage is the multi-locus model MuLM. Compared with the single-locus model, the multi-locus model can significantly reduce false positives. There are currently three multi-locus models. The first is MLMM

[12] , 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] , still based on MLM; using constrained maximum likelihood optimization to replace the standard of K explaining variance close to zero; using fixed effects and random effects alternately, independent random effects model to screen sites, independent fixed effects to verify sites, and using both alternately until no new candidate markers enter the model; because there are no fixed and random effects in the model at the same time, the confounding of different effects is avoided, and false positives and false negatives are controlled at the same time. The third one is BLINK

[14] , no longer uses MLM; uses the Bayesian information criterion to replace the constrained maximum likelihood estimation in the random effect, 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) The false positive rate is not low enough, the statistical power is not high enough, and the comprehensive evaluation index Matthews correlation coefficient MCC still has room for improvement. In short, there is no best, only better. 2) The accuracy rate (Top1 Acc%) of the first PseQTN selected is low. Generally speaking, the first PseQTN represents a gene with a stronger effect and is more important, and is the preferred object for functional verification. However, when the heritability is 0.8, 0.5, and 0.2 in random scenarios, Blink's Top1 Acc is only 89%, 76%, and 60%. 3) In complex scenarios such as TopKSNPs not containing QTNs and tight linkage between QTNs, Blink's MCC and Top1 Acc drop significantly. Affected by LD and low heritability, SNPS ranked high in simple models are not necessarily QTNs. For example, when the heritability is 0.2, R 2 The probability that the Top100SNPs do not contain QTN is 7.80%; at this time, Blink's MCC and Top1 Acc are as low as 0.0338 and 5%. 4) There is still room for improvement in the computational time complexity on large data sets.

[0010] There are many reasons for the above shortcomings: 1) GWAS involves the intersection of multiple disciplines such as bioinformatics, machine learning, quantitative genetics, etc., and the required background knowledge is complex and the threshold is high; it requires comprehensive system consideration, overall improvement and optimization, and taking into account both accuracy and computational efficiency; local improvements are often only suitable for specific scenarios and have poor universality. Even the present invention does not involve epistatic effects for the time being. 2) The influence of inertial thinking. Due to the constraints of computational efficiency, GWAS before 2012 were all single-site models; in the single-site model, the introduction of Q and K was considered a milestone. [5] , which led to the subsequent improvements to continue using Q and K, focusing only on model solution optimization and K construction optimization. So far, the MLM model is still the mainstream. Even among the only three multi-site models, MLMM and FarmCPU are still based on MLM, and initially still include Q and K; it was not until 2019 that BLINK abandoned K for the first time (and currently the only one). 3) It requires years of persistent persistence (the invention took nearly four years to be basically formed), the horizontal transplantation of knowledge from multiple disciplines, the concerted efforts of the team, and the sensitivity to the feedback results of countless failed attempts and partial improvement feedback results in order to clarify the direction of the next exploration. It requires the courage and confidence to abandon Q and K, which predecessors considered to be milestones.

[0011] References (the literature that needs to be cited later is also listed here):

[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. Principle, optimization and application of mixed 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 the present invention belongs to the fourth-stage multi-site MuLM model, pioneers a forward feature selection oscillation search strategy, and is completed by adopting a new termination criterion.

[0046] The present invention provides a method for searching genome-wide association analysis using quantitative trait multi-locus oscillation, the goal of which is to find SNP sites associated with a specific phenotype, namely PseQTNs, comprising the following steps:

[0047] S1 Data preparation: Prepare SNP matrix file and phenotype matrix file; the SNP matrix file records genotype data, and the phenotype matrix file records quantitative trait data.

[0048] S2 Data preprocessing: remove duplicates from the input SNP matrix, standardize by column, and perform principal component analysis, and standardize the phenotype matrix by column;

[0049] S3 association analysis: Through the oscillation search strategy, SNP sites associated with the phenotype, namely PseQTNs, were screened out from the SNP matrix.

[0050] Specifically, the SNP matrix file supports .bed, .raw files, 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] Phenotype matrix files support .csv files. The first row is the phenotype name, the first column is the sample name, and the subsequent columns are continuous phenotype values. It supports sequential analysis of multiple columns of phenotypes.

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

[0053] 1) Calculate the R between the SNP matrix S and a phenotype Y in the phenotype matrix 2 (Y,S i ) and arrange them in descending order, and record the first SNP as * S1; with |r( * S1,S i )|<0.7Filter the Rank and record the second SNP as * S2; and so on, until 10 relatively independent initial SNPs are obtained, recorded as * S k , k = 1, 2, ..., 10, as the starting points of 10 chains;

[0054] 2) Performing principal component analysis on the initial SNPs, taking the first several principal components (e.g., 2 or 3) according to the preset parameters, and recording them as PCs;

[0055] 3) For each chain, PCs+ * Sk is a covariate, for Y and S i Find the partial correlation and record the S corresponding to the one with the largest absolute value of partial correlation i for k S1; called chain extension step; then PCs+ * S k + k S1 is a covariate, which affects Y and S i Find the partial correlation and record the S corresponding to the one with the largest absolute value of partial correlation i for k S2; and so on, a total of 10 steps;

[0056] 4) For 10 chains, each containing PCs and 11 SNPs, calculate the multiple correlation, and the chain corresponding to 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 is a covariate, for Y and S i Find the partial correlation, which is called chain repair; record the S corresponding to the one with the largest absolute value of partial correlation i For Q1; obtain the first PseQTN;

[0058] 6) Using PCs+Q1 as a covariate, perform chain extension for 5 steps as in step 3), and perform 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 perform chain repair as in step 5) to obtain Q3;

[0060] 8) And so on, 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 m} is PseQTNs, find the complex correlation R between Y and PseQTNs+PCs 2 PCs+QTNs , Y and the complex correlation of QTNs R 2 QTNs ;

[0062] 10) Repeat steps 1)-9) until all columns of the phenotype matrix are calculated in sequence, and output the results.

[0063] The present invention provides a system for searching genome-wide association analysis using quantitative trait multi-site oscillation, which comprises 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 execute step S2 of the method described above; and the association analysis module is used to execute step S3 of the method described above.

[0065] Specifically, the data entry module is used to output SNP matrix files and phenotype matrix files. More specifically, the SNP matrix file supports .bed, .raw files, 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 phenotype matrix file supports .csv files, the first line is the phenotype name, the first column is the sample name, and the subsequent columns are continuous phenotype values, supporting the analysis of multiple columns of phenotypes in sequence.

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

[0067] The present invention also provides a device containing the system, which includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program, and the computer program code implements the system.

[0068] Optionally, it also includes a result display device, such as a mobile phone, a computer or a network remote result display.

[0069] The method and system of the present invention have the significant advantages of high accuracy (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 between QTNs and TopKSNPs without QTNs), and high computational efficiency (applicable to large samples and high-density SNPs data). The details are as follows:

[0070] 1) High precision and strong robustness. The first forward feature selection oscillation search strategy, including: a) The "borrow sheep and divide the sheep" strategy is adopted for the population structure Q matrix. Borrow Q to control false positives, and return Q to improve statistical power. The introduction of Q as a covariate helps to select QTNs, but when PseQTNs are selected and the heritability explained by Q+PseQTNs is close to that explained by PseQTNs, Q should be withdrawn after completing its historical mission. b) Tree search. Do not hastily determine that a certain SNP is a PseQTN, but take a few steps and search a few steps (chain extension) first, and then look back to evaluate whether the originally selected SNP is credible (chain repair). c) Chain repair correction. QTNs or QTN's strong LD-SNPs are added as covariates to improve the accuracy of the selection of other parts of QTNs with higher effects. d) Under the combination of multiple covariates, only SNPs are sorted instead of deleted, avoiding the quagmire of multiple tests and significance thresholds, ensuring that weak-effect QTNs have the opportunity to be selected at each step of chain extension, effectively improving statistical power. e) Semi-correlation and LD candidate SNPs, partial correlation to determine PseQTNs. This avoids missed selection when there is strong LD between QTNs and ensures efficient calculation. f) Nested linkage disequilibrium LD candidate multiple initial points, increasing the diversity of initial points makes the model still robust in complex scenarios where TopKSNPs does not contain QTNs.

[0071] 2) Computationally efficient. a) Does not contain K. The random effects defined by K need to use the maximum likelihood method or constrained maximum likelihood method to estimate parameters, which is computationally intensive. At least for simulated data, the genetic model of phenotype Y itself does not contain K and Q, so it is not necessary to include K and Q in the final genetic model. Q (PCs) and K are derived from all SNPs, and the confounding between population structure Q, kinship K, and quantitative trait nucleotides (QTNs) will lead to false negatives. There have been similar attempts to remove K: MLMM compresses K to 0, FarmCPU uses constrained maximum likelihood optimization to replace the standard that the variance explained by K is close to zero, and BLINK abandons K and no longer uses the MLM model. b) Use incomplete tree search (width W = 10, depth D = 10, 10 chains). Complete tree search (width W = 10, depth D = 10, 10 10 c) Perform semi-correlation and LD on all SNPs using matrix vectorization operations to determine a small number of candidate SNPs; perform partial correlation on only a small number of candidate SNPs to determine PseQTNs.

[0072] 3) New termination criteria. Multiple linear regression multiple correlation coefficient R 2 Critical value R 2 (ɑ,ν,M) It is related to the degree of freedom ν, the number of variables m and the significance level ɑ. Obviously, R 2 The critical value is affected by both the number of samples and the number of variables. For forward feature selection, the present invention has found through a large number of numerical experiments that the critical value ΔR2 =R 2 (ɑ,ν,m+1) -R 2 (ɑ,ν,m) It is only affected by the number of samples and ɑ. Let ΔR 2 =

[0073] TV / n, then TV is only related to ɑ. 2 When more PseQTNs are introduced to ensure high statistical potential, ɑ should be large in width and TV should be small; 2 When PseQTNs are introduced, false positives should be controlled more strictly, ɑ should be small, and TV should be large. 2 =1 hour TV best = 1 is a constraint, search for different h 2 The empirical formula for the optimal value of TV is: TV best =-20×h 2 +21;h 2 When the number of features increases from m to m+1 When ΔR 2 <TV best / n terminates the import. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 Manhattan plot and QTNs ranking changes under different covariates. Note: 11,089 SNPs, 1,768 samples. Simulation phenotype settings: QTNs = 10, h 2 =0.8, additive effect, effect value is normally distributed. QTN number represents the effect, QTN1 has the largest effect. PCs=3.

[0075] Figure 2 MOST-GWAS flow chart of the present invention. DETAILED DESCRIPTION

[0076] The present invention is described below through specific implementation modes in order to better understand the present invention, but it does not constitute a limitation of the present invention.

[0077] 1. Abbreviations, full names and explanations of relevant terms:

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

[0079] Blink:Bayesian information and linkage disequilibrium iterativelynested keyway,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 GWAS algorithm evaluation indicator.

[0082] GLM: General linear model, a genetic model.

[0083] GP / GS: Genomic prediction / selection. Using a training population with known genotypes and phenotypes to model, predict the phenotype of a test population with only known genotypes.

[0084] GWAS: Genome-wide association studies, a statistical analysis method to find QTNs or genes associated with phenotypes from a large number of genotypes (SNPs).

[0085] h 2 :Heritability, heritability.

[0086] K: Kinship, the kinship matrix between individuals.

[0087] LD: Linkage disequilibrium. The non-random association between alleles at different loci, manifested in the probability of two genes being inherited simultaneously in a population being greater than the probability of their random combination being inherited.

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

[0089] MLM: Mixed linear model, a linear model with fixed effects 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, quantitative trait multi-locus oscillatory search GWAS-regression. This patent application proposes a quantitative trait multi-locus GWAS algorithm.

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

[0093] PC:Principal component (analysis), a linear dimensionality reduction algorithm.

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

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

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

[0097] QTN: Quantitative trait nucleotide. The SNP within the QTL that actually contributes to the variation of the quantitative trait.

[0098] R 2 :Squared Pearson correlation coefficient, square of correlation coefficient. Measures the linear correlation between two continuous variables and is a commonly used evaluation indicator for regression algorithms in machine learning.

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

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

[0101] TV: Threshold value, threshold.

[0102] Example

[0103] 1. Mathematical description of the problem

[0104] Existing datasets (Y i ,S ij ), i = 1, 2, ..., n; n is the number of samples. j = 1, 2, ..., m; m is the number of SNPs. Y is the phenotype. Assume that there are only m' linear combinations of SNPs (QTNs) that constitute Y and there is noise (heritability h 2 <1), how to find these m' SNPs (QTNs) as accurately as possible from m SNPs?

[0105] 2. The invention's entry point

[0106] SNPs are divided into four categories: the first category is TP-SNPs (True Positive), namely QTNs. Note that there may be different degrees of redundancy between TP-SNPs due to linkage disequilibrium LD. The second category is false positive PS-SNPs caused by population stratification PS; population stratification is the systematic difference in allele frequencies between subpopulations in a population due to non-random mating between individuals, which leads to the selection of a large number of false positive sites in GWAS; using Q (PCs) as a covariate is currently the most commonly used control population stratification method. 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 TP-SNP often has multiple LD-SNPs, which is the main cause of false positives. The fourth category is irrelevant markers Ir-SNPs, which are not related to the phenotype themselves. Due to low heritability, small samples, high density SNPs and other reasons, some Ir-SNPs may pass the hypothesis test.

[0107] The present invention is based on three facts often observed in practice (such as Figure 1 shown):

[0108] 1) Top-ranked SNPs in the simple model S Not necessarily QTNs. Figure 1 Among the 10 QTNs in A, QTN9, which has the strongest effect, ranks 5, and the 1st one is the strong LD marker SNP of QTN9. 10118 .

[0109] 2) Adding QTNs as covariates can improve the ranking of other QTNs with higher effects. Figure 1 B to Figure 1The ranking of QTN1 increased from 44 to 19, and the ranking of QTN2 increased from 19 to 1.

[0110] 3) The strong LD-SNPs of QTNs were added as covariates, which could also improve the ranking of some high-effect QTNs. Figure 1 B to Figure 1 The ranking of QTN1 also increased from 44 to 19, and the ranking of QTN2 also increased from 19 to 1 (the genotype data are from the literature

[15] ).

[0111] Therefore, a possible subset of PseQTNs (chain extension) can be obtained by forward search, for example: PCs+SNP 10118 +QTN2+SNP 3657 +QTN1; then {PCs, SNP 3657 , QTN1} is the first selected SNP for covariate correction 10118 (Chain repair), at this time QTN9 is likely to be selected correctly. This strategy is collectively called oscillation search.

[0112] 3. Algorithm steps of the present invention

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

[0114] 1) Matrix operation to find R 2 (Y,S i ) and arrange them in descending order, and record the first SNP as * S1; with |r( * S1,S i )|<0.7Filter the Rank and record the second SNP as * S2; and so on, until 10 relatively independent initial SNPs are obtained, recorded as * S k , k = 1, 2, ..., 10, as the starting point of 10 Chains.

[0115] 2) Perform principal component analysis on SNPs, and select several (e.g., 2 or 3, 2 by default) principal components to represent the Q matrix according to preset parameters, denoted as PCs.

[0116] 3) For each chain, PCs+ * S k is a covariate, for Y and S i Find the partial correlation and record the S corresponding to the one with the largest absolute value of partial correlation i for k S1; called chain extension step. Then PCs+ * Sk + k S1 is a covariate, which affects Y and S i Find the partial correlation and record the S corresponding to the one with the largest absolute value of partial correlation i for k S2; and so on, extending a total of 10 steps.

[0117] 4) For 10 chains, each containing PCs and 11 SNPs, calculate the multiple correlation, and the chain corresponding to 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 is a covariate, for Y and S i Find the partial correlation, which is called chain repair. The S corresponding to the one with the largest absolute value of partial correlation is recorded as i For Q1; obtain the first PseQTN.

[0119] 6) PCs+Q1 was used as a covariate, the chain was extended for 5 steps, and the chain was repaired to obtain Q2.

[0120] 7) PCs+Q1+Q2 were used as covariates, the chain was extended for 5 steps, and the chain was repaired to obtain Q3.

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

[0122] 9) Note {Q1+Q2+…+Q m’} is PseQTNs, find the complex correlation R between Y and QTNs+PCs 2 PCs+QTNs , Y and the complex correlation of QTNs R 2 QTNs .

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

[0124] 4. Model evaluation criteria

[0125] Existing standards such as statistical power, type I error, FDR (false positive number / total positive number), relative accuracy (Relative Acc%), and absolute accuracy (Absolute Acc%) are all evaluated from a single perspective and are flawed. For example, if QTNs = 20, and a model only finds one and it is correct, then the relative accuracy is 100%, which is very high, but the statistical power and absolute accuracy are very low; if a model finds 10,000 and contains 20 true QTNs, then the statistical power and absolute accuracy are very high, but the relative accuracy is very low, the type I error rate is very high, and the FDR is very high.

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

[0127]

[0128] When the positive and negative samples are unbalanced, the accuracy is flawed: for example, if there are 10,000 SNPs and QTNs=10, a model does not find any (all are judged as negative), and the accuracy is as high as 99.90%. In GWAS, the positive and negative samples are extremely unbalanced, and MCC takes into account the sample imbalance.

[0129]

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

[0131] False Negative FN (False Negative): The SNP is a true QTN, but the model does not judge it as QTN.

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

[0133] True Negative TN: This SNP is a false QTN and is not judged as a QTN by the model.

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

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

[0136] Verification Example 1: Simulated data, performance of each model under different heritability in random scenarios

[0137] Using the winter wheat dataset

[16] , 11089 SNPs, 1768 samples. The random scenario means that the QTNs are randomly located on the chromosomes without clustering. Each simulation randomly selects 20 SNPs as QTNs, that is, only 20 SNPs linear combinations constitute phenotype Y, and the additive effect values ​​of the 20 QTNs conform to the normal distribution. Heritability h 2 Set three levels: 0.8, 0.5, and 0.2. 2 =0.8 means that 20 QTNs explain 80% of the variation of the phenotype. The model is required to find a subset of PseQTNs from 11089 SNPs, compare it with the set QTNs subset, and calculate the Matthews correlation coefficient MCC, the first PseQTN accuracy Top1 Acc%, and the running time. Repeat the simulation 100 times, and the results are expressed as the mean. The reference models include the single-site model representative SUPER, and the existing three 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 with the single-site model SUPER, the four multi-site models all 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, Top 100 excluding QTN, performance of each model under different heritability in complex scenarios

[0143] Top ranked SNPs by simple model S It does not necessarily contain QTNs. Simulation results show that when the heritability is 0.2, 0.5, and 0.8, R 2 The probability that the top 100 SNPs do not contain QTNs is 7.80%, 3.90%, and 3.50%, respectively; when the heritability is low, R 2 It is not uncommon for the top 100 to not contain QTNs. The first PseQTN selected by the model often represents a gene with a stronger effect and is more important, and is the first choice for functional verification. 2 In the complex scenario where the Top100 SNPs do not contain QTNs, it is particularly important to evaluate the performance of each model, especially Top1Acc%.

[0144] Using the winter wheat dataset

[16] The simulation process is the same as in Example 1, and the simulation is repeated 3000 times, from which 100 R 2 The Top100 SNPs without QTNs dataset was used in this example.

[0145] Table 2. Performance of each model under different heritability in the Top 100 complex scenario without QTN

[0146]

[0147] Conclusion: MOST-GWAS has the highest MCC and Top1 Acc%. Compared with the random scenario (Table 1), the advantages of MOST-GWAS over the reference model in the complex scenario where the Top100 SNPs do not contain QTN are extremely obvious.

[0148] Verification Example 3: Simulated data, performance of each model under different heritability in a complex scenario with tight linkage of QTN

[0149] Using the winter wheat dataset

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

[0150] Table 3. Performance of each model under different heritability in the complex scenario of QTN tight linkage

[0151]

[0152]

[0153] Conclusion: MOST-GWAS has the highest MCC and Top1 Acc%. Compared with the random scenario (Table 1), the advantage of MOST-GWAS over the reference model is more obvious in the complex scenario of QTN tight linkage.

[0154] Verification Example 4: Simulated data, comparison of Blink and MOST-GWAS multi-species performance in random scenarios

[0155] To understand the universality of the model, we further compared the performance of the Blink and MOST-GWAS models in five different species random scenarios. The winter wheat dataset has 11,089 SNPs and 1,768 samples.

[16] , the rice dataset has 36,898 SNPs and 413 samples

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

[18] There are 214,545 SNPs and 199 accessions (https: / / archive.gramene.org / diversity / download_data.html), mouse dataset

[15] There are 12,226 SNPs and 1,940 samples (1,000 male samples and 940 female samples), Holstein cattle dataset

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

[0156] Twenty SNPs were randomly selected as QTNs in each simulation for each species and each heritability, and the simulation was repeated 100 times, and the results were expressed as the mean.

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

[0158]

[0159]

[0160] Conclusion: For 5 different species, the MCC and Top1 Acc% of MOST-GWAS are stably better than those of Blink.

[0161] Verification Example 5: Maize Flowering Time GWAS (real data)

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

[20] (http: / / cbsusrv04.tc.cornell.edu / users / panzea / filegateway.aspx?category=Genotypes), and the phenotype is flowering time, i.e., silking days (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. Comparison 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 aligned in the sequence; the closest one was 528 bp.

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

[0170]

[0171] 49 PseQTNs explained 73.79% of the heritability, and 17 reported maize flowering genes [21-33] Four known flowering genes were mapped in the PCR, and three times two QTNs were mapped to the same gene; the closest one was 6500 bp.

[0172] MOST-GWAS mapped 42 PseQTNs to 6 known genes, the nearest 528 bp, while Blink mapped 49 PseQTNs to only 4 known genes, the nearest 6500 bp. MOST-GWAS is significantly better than Blink.

Claims

1. A method for searching genome-wide association analysis using quantitative trait multi-locus oscillation, characterized in that: The method aims to find SNP sites associated with a specific phenotype, namely PseQTNs, and comprises the following steps: S1 Data preparation: Prepare SNP matrix files and phenotype matrix files; the SNP matrix file records genotype data, and the phenotype matrix file records quantitative trait data; S2 Data preprocessing: remove duplicates from the input SNP matrix, standardize by column, and perform principal component analysis, and standardize the phenotype matrix by column; S3 association analysis: Through the oscillation search strategy, SNP sites associated with the phenotype, namely PseQTNs, were screened out from the SNP matrix.

2. The method according to claim 1, characterized in that The SNP matrix file supports .bed, .raw files, 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; Phenotype matrix files support .csv files. The first row is the phenotype name, the first column is the sample name, and the subsequent columns are continuous phenotype values. It supports sequential analysis of multiple columns of phenotypes.

3. The method according to claim 1, characterized in that The method of S3 association analysis is as follows: 1) Calculate the R between the SNP matrix S and a phenotype Y in the phenotype matrix 2 (Y,S i ) and arrange them in descending order, and record the first SNP as * S1; with |r( * S1,S i )|<0.7Filter the Rank and record the second SNP as * S2; and so on, until 10 relatively independent initial SNPs are obtained, recorded as * S k , k = 1, 2, ..., 10, as the starting points of 10 chains; 2) Performing principal component analysis on the initial SNPs, taking the first several principal components (e.g., 2 or 3) according to the preset parameters, and recording them as PCs; 3) For each chain, PCs+ * S k is a covariate, for Y and S i Find the partial correlation and record the S corresponding to the one with the largest absolute value of partial correlation i for k S1; called chain extension step; then PCs+ * S k + k S1 is a covariate, which affects Y and S i Find the partial correlation and record the S corresponding to the one with the largest absolute value of partial correlation i for k S2; and so on, a total of 10 steps; 4) For 10 chains, each containing PCs and 11 SNPs, calculate the multiple correlation, and the chain corresponding to the largest absolute value of the multiple correlation is the optimal chain; 5) For the optimal chain, cut off * S+S1, PCs+S2+S3+…+S 10 is a covariate, for Y and S i Find the partial correlation, which is called chain repair; record the S corresponding to the one with the largest absolute value of partial correlation i is Q1; Obtained the first PseQTN; 6) Using PCs+Q1 as covariate, perform chain extension for 5 steps as in step 3), and perform chain repair as in step 5) to obtain Q2; 7) Using PCs+Q1+Q2 as covariates, perform chain extension for 5 steps as in step 3), and perform chain repair as in step 5) to obtain Q3; 8) And so on, until the chain extension termination criterion is met: ΔR 2 <TV best / n, TV best =-20×h 2 +21; 9) Let {Q1, Q2, ..., Q m } is PseQTNs, find the complex correlation R between Y and PseQTNs+PCs 2 PCs+QTNs , Y and the complex correlation of QTNs R 2 QTNs ; 10) Repeat steps 1)-9) until all columns of the phenotype matrix are calculated in sequence, and output the results.

4. A system for searching genome-wide association analysis using quantitative trait multi-locus oscillation, 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 execute step S2 of the method described in any one of claims 1 to 3; and the association analysis module is used to execute step S3 of the method described in any one of claims 1 to 3.

5. The system according to claim 4, characterized in that The data entry module is used to output SNP matrix files and phenotype matrix files.

6. The system according to claim 5, characterized in that The SNP matrix file supports .bed, .raw files, 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; Phenotype matrix files support .csv files. The first row is the phenotype name, the first column is the sample name, and the subsequent columns are continuous phenotype values. It supports sequential analysis of multiple columns of phenotypes.

7. The system according to claim 4, characterized in that The result output module outputs the results as follows: corresponding to the number of phenotypes in the input phenotype matrix, an Excel file is output for each phenotype, which contains two sheets. The first sheet contains SNPID, P-Value, and estimated effect value (only for PseQTNs), and the second sheet contains the heritability explained by PseQTNs.

8. A device containing the system according to any one of claims 4 to 7, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program, and the computer program code implements the system according to any one of claims 4 to 7.

9. The device according to claim 8, characterized in that Also includes result display equipment.

10. The device according to claim 9, characterized in that The result display device is a mobile phone, a computer or a network remote result display.

Images