A method for constructing a multi-gene genetic risk assessment model using an ising model

Abstract

The application discloses a kind of construction methods of multi-gene genetic risk assessment model using Ising model, including reading SNP data of original sample, the correlation coefficient between each single nucleotide polymorphism is calculated;Obtain the public data of genome-wide association analysis of target disease, set the threshold of P value as the basis for subsequent screening SNP:According to the correlation coefficient of each SNP in original sample, construct Ising model, and add regular term;Using Ising machine to solve the solution of Ising model, obtain SNP for predicting target disease;Based on the screened SNP, the target disease is predicted using machine learning algorithm.The construction method of multi-gene genetic risk assessment model using Ising model disclosed by the application converts the process of screening SNP in traditional method into the solution of Ising model, and uses Ising machine to efficiently and accurately screen SNP for algorithm modeling PRS, thereby improving the prediction efficiency and accuracy of target disease.

Description

Technical Field

[0001] This invention belongs to the field of genetic risk assessment model construction technology, specifically relating to a method for constructing a polygenic genetic risk assessment model using the Ising model. Background Technology

[0002] Diseases are caused by a complex interplay of factors, which can be broadly categorized into genetic and environmental factors. To predict the genetic risk of a target disease from an individual's genetic data, many research institutions utilize large-scale patient samples to conduct association analyses between single nucleotide polymorphisms (SNPs) and phenotypic traits, and publicly release the results of genome-wide association studies (GWAS). For institutions with only small-scale sample data, they can construct polygenic risk score (PRS) models based on GWAS data published by large research institutions and validate these PRS models on smaller sample data.

[0003] However, even using the weighted effect values ​​of individual SNPs in publicly available GWAS data, it remains difficult to effectively predict phenotypic traits. One challenge is that most SNPs contribute very little to phenotypic traits, but their cumulative effect often has a relatively significant clinical impact. Therefore, in the actual construction of PRS, whether to select the most significant few SNPs or a large range of low-significance SNPs needs to be adjusted according to the specific circumstances. Another challenge is that SNPs that are relatively close together often have a very strong correlation (linkage disequilibrium, LD). Including strongly correlated SNPs in PRS modeling simultaneously will hinder the construction of the predictive model, leading to a decrease in predictive accuracy. Therefore, traditional PRS construction methods screen out strongly correlated SNPs before modeling.

[0004] Therefore, further improvements will be made to address the aforementioned issues. Summary of the Invention

[0005] The main objective of this invention is to provide a method for constructing a polygenic genetic risk assessment model using the Ising model. The method constructs the Ising model based on the correlation coefficient matrix of each SNP in the sample data, using 0 or 1 to represent the screening result of the SNP. The Ising model is solved efficiently and accurately using the Ising machine. Based on the SNPs screened by the Ising model, machine learning algorithms are used to model the target phenotypic traits, thereby improving the prediction efficiency and accuracy of PRS.

[0006] To achieve the above objectives, this invention provides a method for constructing a polygenic genetic risk assessment model using the Ising model, comprising the following steps:

[0007] Step S1: Read the genomic data of the sample and read the external reference data at the same time to perform SNP site quality control on the genomic data of the sample;

[0008] Step S2: Obtain GWAS statistics for the target phenotypic trait;

[0009] Step S3: Based on the acquired SNP locus data, calculate the correlation coefficient between each SNP locus;

[0010] Step S4: Based on the P-values ​​of each SNP in the GWAS statistics obtained in step S2, set a threshold for the P-values ​​used to filter SNPs.

[0011] Step S5: Generate the Ising model based on the correlation coefficient between the quality control and screened SNPs;

[0012] Step S6: Add a regularization term to the Ising model generated in step S5;

[0013] Step S7: Solve the Ising model after adding the regularization term in step S6 using the Ising machine;

[0014] Step S8: Based on the solution of the Ising model obtained in step S7, determine the SNP sites used for PRS modeling;

[0015] Step S9: Based on the selected SNP sites, use machine learning to model PRS to predict the target phenotypic trait.

[0016] As a further preferred technical solution to the above technical solution, for step S1:

[0017] The sample gene data read can be the result of whole genome sequencing (WGS) or whole exome sequencing (WES), or the result obtained based on Affy500k or Illumina330 (etc.) SNP chips (etc.).

[0018] When performing SNP quality control, SNPs with frequencies less than a preset value can be removed, as well as SNPs that fail the Hardy-Weinberg Equilibrium (HWE) test. The reference data for SNP quality control can be external publicly available 1000 Genomes datasets or sample gene data.

[0019] As a further preferred technical solution to the above technical solution, the Ising model is constructed based on the correlation coefficient between SNPs to screen SNPs.

[0020] As a further preferred technical solution to the above technical solution, SNPs are screened by constructing the Ising model for predicting phenotypic traits.

[0021] As a further preferred technical solution to the above technical solution, in step S6, an L1 or L2 regularization term is added to the constructed Ising model.

[0022] As a further preferred technical solution to the above technical solution, the Ising machine is used to solve the Ising model to select the most suitable SNP for predicting phenotypic shape.

[0023] As a further preferred technical solution to the above technical solution, the Ising model is constructed based on the correlation coefficient between SNPs, and the Ising machine is used to solve it, including those based on optical chips, FPGAs, superconducting quantum mechanics, ion traps, CPUs, and GPUs.

[0024] As a further preferred technical solution to the above technical solution, after screening suitable SNPs for prediction based on the Ising model, machine learning algorithms including Lasso, Ridge regression, Elastic-Net, Bayesian neural network and Hidden Markov model are used to predict phenotypic traits.

[0025] After obtaining multiple solutions based on the Ising model, the optimal solution is determined as the output model based on the prediction accuracy of the training sample set or multiple cross-validation subsets.

[0026] As a further preferred technical solution to the above technical solution, the GWAS statistical results include information such as the chromosome number of each SNP, its position in the chromosome, P value, beta value or OR value, and SE (standard deviation).

[0027] To achieve the above objectives, the present invention provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the method for constructing a polygenic genetic risk assessment model using the Ising model.

[0028] To achieve the above objectives, the present invention provides a non-transitory computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the steps of the method for constructing a polygenic genetic risk assessment model using the Ising model.

[0029] The beneficial effects of this invention are as follows:

[0030] This invention transforms the process of selecting SNPs into solving the Ising model, and utilizes the high-speed solution capability of the Ising machine for combinatorial optimization problems to select the most suitable SNPs for PRS modeling more efficiently and accurately than traditional methods, thereby improving the efficiency and accuracy of PRS modeling. Attached Figure Description

[0031] Figure 1 This is a flowchart of a method for constructing a polygenic genetic risk assessment model using the Ising model according to the present invention.

[0032] Figure 2 This is an example Manhattan plot of GWAS statistics for obesity, calculated by taking the P-value as -log10(P), in a method for constructing a polygenic genetic risk assessment model using the Ising model according to the present invention.

[0033] Figure 3 This is a diagram illustrating the linkage disequilibrium effect among SNPs on chromosome 19 in a method for constructing a polygenic genetic risk assessment model using the Ising model according to the present invention. Detailed Implementation

[0034] The following description is intended to disclose the present invention and enable those skilled in the art to implement it. The preferred embodiments described below are merely examples, and other obvious variations will occur to those skilled in the art. The basic principles of the invention defined in the following description can be applied to other embodiments, modifications, improvements, equivalents, and other technical solutions that do not depart from the spirit and scope of the invention.

[0035] In a preferred embodiment of the present invention, those skilled in the art should note that the genomic data and the like involved in the present invention can be considered as prior art.

[0036] Preferred embodiment.

[0037] This invention discloses a method for constructing a polygenic genetic risk assessment model using the Ising model, comprising the following steps:

[0038] Step S1: Read the genomic data of the sample and read the external reference data at the same time to perform SNP site quality control on the genomic data of the sample;

[0039] Step S2: Obtain GWAS statistics for the target phenotypic trait;

[0040] Step S3: Based on the acquired SNP locus data, calculate the correlation coefficient between each SNP locus;

[0041] Step S4: Based on the P-values ​​of each SNP in the GWAS statistics obtained in step S2, set a threshold for the P-values ​​used to filter SNPs.

[0042] Step S5: Generate the Ising model based on the correlation coefficient between the quality control and screened SNPs;

[0043] Step S6: Add a regularization term to the Ising model generated in step S5;

[0044] Step S7: Solve the Ising model after adding the regularization term in step S6 using the Ising machine;

[0045] Step S8: Based on the solution of the Ising model obtained in step S7, determine the SNP sites used for PRS modeling;

[0046] Step S9: Based on the selected SNP sites, use machine learning to model PRS to predict the target phenotypic trait.

[0047] Specifically, for step S1:

[0048] The sample gene data read can be the result of whole genome sequencing (WGS) or whole exome sequencing (WES), or the result obtained based on Affy500k or Illumina330 (etc.) SNP chips (etc.).

[0049] When performing SNP quality control, you can choose to remove SNPs with a frequency less than a preset value (set according to needs, such as 1%), as well as SNPs that fail the Hardy-Weinberg Equilibrium (HWE) test. When performing SNP quality control, the reference data is external publicly available 1000 Genomes dataset or sample gene data (etc.).

[0050] More specifically, the Ising model is constructed based on the correlation coefficients between SNPs to screen for SNPs.

[0051] Furthermore, SNPs were screened by constructing the Ising model for predicting phenotypic traits.

[0052] Furthermore, in step S6, an L1 or L2 regularization term is added to the constructed Ising model.

[0053] Preferably, the Ising machine is used to solve the Ising model to select the most suitable SNP for predicting phenotypic shape by utilizing its ability to handle combinatorial optimization problems.

[0054] Preferably, the Ising model is constructed based on the correlation coefficient between SNPs, and solved using the Ising machine, which includes those based on optical chips, FPGAs, superconducting quantum mechanics, ion traps, CPUs, and GPUs.

[0055] Preferably, after selecting suitable SNPs for prediction based on the Ising model, machine learning algorithms including Lasso, Ridge regression, Elastic-Net, Bayesian neural network and Hidden Markov model are used to predict phenotypic traits.

[0056] After obtaining multiple solutions based on the Ising model, the optimal solution is determined as the output model based on the prediction accuracy of the training sample set or multiple cross-validation subsets.

[0057] As a supplementary example, for step S1:

[0058] The SNP data of the samples were obtained and divided into training and validation sets. In practical applications, the samples can also be divided into multiple subsets for cross-validation. In this embodiment, the software tool GWAsimulator2 and the public dataset HAPMAP3 were used to simulate 40,000 SNPs located on chromosome 3, generating a training set of 2,000 participants and a validation set of 1,000 participants. A dataset of 5,000 participants was also generated for GWAS statistical analysis. Several SNPs (e.g., 400, 40, and 4) were randomly selected as pathogenic loci, and the true effect value β of the pathogenic SNP loci was randomly generated from a Gaussian distribution of ~N(0,1). i The individual's phenotypic traits are calculated using the following formula:

[0059] y j =e j +∑ i β i x i,j ;

[0060] In the above formula, x i,j For the data of the i-th SNP in the j-th individual sample, ∑ i β i x i,j Residual term e represents the genetic factors in the target phenotypic trait. j Represents environmental factors;

[0061] Personal environmental factors The phenotypic shape is randomly generated from a Gaussian distribution to ensure that the heritability h of the simulated phenotypic shape is 0.2.

[0062]

[0063] More specifically, in step S2, the GWAS statistical data includes the chromosome number of each SNP, its position in the chromosome, p-value, beta value or OR value, and SE (standard deviation) information. GWAS statistical results can be downloaded from external databases or analyzed from large sample datasets. In this embodiment, GWAS statistical analysis is performed based on simulation data from 5,000 individuals. In practical applications, GWAS statistical data for the target phenotypic trait can be obtained from institutions such as the UK BioBank and the China Kadoorie Biobank.

[0064] Furthermore, in step S3, when the number of SNP loci is too large, calculating the correlation coefficients between all SNPs would consume too much memory. Therefore, a moving window with a certain distance is set, and only SNPs within the moving window are calculated. In step S3, when calculating the correlation coefficients between each SNP, the moving window can be used to determine the calculation range. For example, on the same chromosome, a 1000×1000 correlation coefficient matrix can be calculated by grouping every 1000 SNPs.

[0065] Furthermore, in step S4, based on the value of -log10(P), the preset interval (preferably [0,10], which can be adjusted according to needs) is divided into several (preferably 50) equal parts as the threshold of -log10(P), such as the thresholds of 1, 2, 3, 4, 5, 6, 7, 8, 9, 10.

[0066] Preferably, for step S5, constructing the Ising model based on the correlation coefficient matrix R, where:

[0067] The following formula is constructed based on the least squares method:

[0068]

[0069] In the above formula, w is a vector of elements taking values ​​of 0 or 1, which is the solution of the Ising model. In the above formula, k is a parameter that can be set according to the actual situation. In this embodiment, k takes the values ​​1, 2, 5, 10, 20, 50, and 100.

[0070] In step S6, an L1 or L2 regularization term is added to the Ising model constructed in step 5. In this embodiment, since w is a vector in which each element takes the value 0 or 1, L2 regularization is essentially the same as L1 regularization.

[0071]

[0072] In the above formula, λ is a parameter that controls the regularization term. In this embodiment, λ is taken as 0.01, 0.03, 0.1, 0.3.

[0073] In step S7, the Ising model obtained in step S6 is input to the Ising machine for solution. The Ising machine can solve combinatorial optimization problems with results of (0,1) at high speed. The Ising machine can be a computing device based on FPGA or optical chip, a computing device based on quantum systems such as low-temperature superconductors or ion traps, or an algorithm system based on CPU or GPU. In this embodiment, the simulated annealing algorithm is run on a CPU to solve the Ising model.

[0074] In step S8, based on the solution obtained from the Ising machine in step S7, the SNP used for predicting and modeling the target phenotypic trait is determined. In this embodiment, since 50 threshold values ​​are set for the P value, 7 values ​​are taken for the parameter k, and 4 values ​​are taken for the parameter λ, there are a total of 1400 solutions.

[0075] Preferably, in step S9:

[0076] In this embodiment, based on the solutions of 1400 Ising models, the beta value of the corresponding SNP in the GWAS statistical data is used to calculate the PRS prediction value according to the following formula:

[0077]

[0078] In the above formula, x i,j For the i-th SNP in the j-th individual sample, This refers to the beta value of the corresponding SNP in the GWAS statistics. This represents the predicted PRS value for the j-th individual sample. It should be noted that the target phenotypic trait in this embodiment is a numerical indicator; therefore, the above formula is a linear regression. In practical applications, if the target phenotypic trait is a case-control sample indicating whether or not the individual has the disease, logistic regression can be used using the following formula, where α is a constant in logistic regression:

[0079]

[0080] In step S9, based on the SNP selected in step S8, machine learning algorithms including LASSO, Ridge regression, random forest, Lassosum, LDpred, and neural networks are used to predict the target phenotypic trait; alternatively, the beta value or OR value of the corresponding SNP in the GWAS statistical data can be used directly. The target phenotypic trait can be a disease including diabetes or hypertension, or it can be an indicator value such as fasting blood glucose or systolic blood pressure.

[0081] The accuracy of 1400 predicted values ​​was calculated on both the training and validation sets. The set with the highest accuracy in the training set was taken as the final result, and its accuracy on the validation set was verified. In practical applications, the sample data can also be divided into multiple subsets for cross-validation. The standard for calculating accuracy can be the coefficient of determination between the predicted and true values ​​when the target phenotype is a numerical indicator, and the ROC-AUC value between the predicted and true disease status when the target phenotype is a case-control sample. In this embodiment, the number of pathogenic SNP sites was set to 400, 40, and 4 during the simulation data generation process. The validation accuracy under the three conditions is shown in Table 1 below.

[0082] Table 1

[0083] 400 0.0945 0.0985 40 0.0989 0.1086 4 0.1546 0.1553

[0084] It should be noted that in the traditional P+T method, which is used as a comparison object, 50 thresholds are also set based on the P-values ​​of GWAS statistics. However, when dealing with the cascading imbalance effects between SNPs, the squared correlation coefficient r is removed. 2 For SNPs with values ​​above 0.01, 0.05, 0.1, 0.2, 0.5, 0.8, and 0.95, the window size for calculating the correlation coefficient between SNPs is set to 50kb (Kilobase), 100kb, 200kb, and 500kb, respectively. The traditional P+T method yields a total of 1400 solutions, directly using the beta value of the corresponding SNP from the GWAS statistics, and selecting the solution with the highest accuracy in the training set as the final result. This embodiment leverages the Ising mechanism's ability to efficiently and accurately solve the Ising model when dealing with the chain imbalance effect between SNPs, achieving better prediction accuracy than the traditional P+T method under fair experimental conditions.

[0085] The present invention also discloses an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the method for constructing a polygenic genetic risk assessment model using the Ising model.

[0086] The present invention also discloses a non-transitory computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the steps of the method for constructing a polygenic genetic risk assessment model using the Ising model.

[0087] In this invention, the Ising machine is a solver specifically designed for combinatorial optimization problems, which are difficult to solve in efficient time on classical computers. The challenge in applying the Ising machine to practical problems lies in how to effectively transform the target problem into an Ising model that the Ising machine can handle.

[0088] It is worth mentioning that the technical features such as genomic data involved in this patent application should be regarded as prior art. The specific structure, working principle, and possible control methods and spatial arrangement of these technical features can be adopted using conventional choices in the field, and should not be regarded as the inventive point of this patent. This patent will not be further elaborated in detail.

[0089] For those skilled in the art, modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this invention should be included within the protection scope of this invention.

Claims

1. A method for constructing a polygenic genetic risk assessment model using the Ising model, characterized in that, Includes the following steps: Step S1: Read the genomic data of the sample and read the external reference data at the same time to perform SNP site quality control on the genomic data of the sample; Step S2: Obtain GWAS statistics for the target phenotypic trait; Step S3: Based on the acquired SNP locus data, calculate the correlation coefficient between each SNP locus; Step S4: Based on the P-values ​​of each SNP in the GWAS statistics obtained in step S2, set a threshold for the P-values ​​used to filter SNPs. Step S5: Generate the Ising model based on the correlation coefficients between the quality control and screened SNPs; Step S5 constructs the Ising model based on the correlation coefficient matrix R, where: The following formula is constructed based on the least squares method: ； In the above formula The vector represents the values ​​of each element as 0 or 1, which is the solution of the Ising model. In the above formula, k is a parameter set according to the actual situation. Step S6: Add a regularization term to the Ising model generated in step S5; Step S7: Solve the Ising model after adding the regularization term in step S6 using the Ising machine; Step S8: Based on the solution of the Ising model obtained in step S7, determine the SNP sites used for PRS modeling; Step S9: Based on the selected SNP sites, use machine learning to model PRS to predict the target phenotypic trait; After selecting suitable SNPs for prediction based on the Ising model, machine learning algorithms including Lasso, Ridge regression, Elastic-Net, Bayesian neural network and Hidden Markov model are used to predict phenotypic traits. After obtaining multiple solutions based on the Ising model, the optimal solution is determined as the output model based on the prediction accuracy of the training sample set or multiple cross-validation subsets.

2. The method for constructing a polygenic genetic risk assessment model using the Ising model according to claim 1, characterized in that, The Ising model is constructed based on the correlation coefficients between SNPs to screen for SNPs.

3. The method for constructing a polygenic genetic risk assessment model using the Ising model according to claim 1, characterized in that, SNPs were screened by constructing the Ising model for predicting phenotypic traits.

4. The method for constructing a polygenic genetic risk assessment model using the Ising model according to claim 1, characterized in that, In step S6, an L1 or L2 regularization term is added to the constructed Ising model.

5. The method for constructing a polygenic genetic risk assessment model using the Ising model according to claim 1, characterized in that, By utilizing the Ising machine's ability to handle combinatorial optimization problems, the Ising model is solved, and the most suitable SNP for predicting phenotypic shape is selected.

6. The method for constructing a polygenic genetic risk assessment model using the Ising model according to claim 1, characterized in that, The Ising model is constructed based on the correlation coefficient between SNPs, and solved using Ising machines in various forms, including those based on optical chips, FPGAs, superconducting quantum mechanics, ion traps, CPUs, and GPUs.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the method for constructing a polygenic genetic risk assessment model using the Ising model as described in any one of claims 1 to 6.

8. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the computer program implements the steps of the method for constructing a polygenic genetic risk assessment model using the Ising model as described in any one of claims 1 to 6.

Images