Gene-phenotype correlation analysis model as well as establishment method and application thereof

Abstract

The invention discloses a gene-phenotype correlation analysis model as well as an establishment method and application thereof, and belongs to the technical field of biological medicines. The model establishment method comprises the following steps: S1, collecting known associated character-gene data to form a gene-character pair; s2, calculating the score of the rare mutation type of each gene by using a formula; s3, analyzing the correlation between the mutation score and the character through linear regression, calculating the weight of each mutation type, optimizing the weight combination, and taking the correlation R2 as an evaluation standard; s4, calculating a rare mutation load score of the sample gene according to a scoring formula and the optimized weight; and S5, analyzing the relevance between the mutation load fraction and the phenotype through a regression method, and constructing a gene-phenotype relevance analysis model. Compared with a traditional gene-base colouring method, the model has better reproducibility and complementarity, and can be used for discovering candidate risk genes of new characters or unknown diseases.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and in particular to a gene-phenotype association analysis model and an establishment method and application thereof. Background Art

[0002] Rare mutations are a key factor influencing diseases and traits. Genetic association analysis based on rare mutations can screen for candidate genes associated with diseases or traits, making it a common approach to identifying disease risk or causative genes. Assessing the pathogenicity or deleteriousness of rare mutations is crucial for identifying associated genes. While numerous prediction software programs are available that can predict different mutations and offer excellent assessment results, their effectiveness is unlikely to be completely equivalent to experimental verification. Furthermore, due to the sheer volume of data, often encompassing tens of millions of mutations, it is unrealistic to experimentally evaluate every locus.

[0003] Traditional gene-base collapsing methods based on rare mutations only consider different levels of mutation type and lack a systematic assessment of the deleteriousness of mutations. Other gene-base analysis methods, such as SAIGE-GENE and SKAT-o, while highly specific, lack sensitivity and may miss some candidate genes. Furthermore, these methods require integrating rare mutations into plink files. Due to the sparseness of rare mutations in the human population, these files are often large, increasing workload, time, and storage costs. Furthermore, current research has yet to clarify which combinations of mutation type weights can achieve similar analytical results, nor which weight combination ultimately yields the best results. Furthermore, comparative studies between various gene analysis methods are insufficient. For example, how newly developed methods compare with gene-based collapsing, which method is superior, and whether they complement each other are among the issues that need to be addressed. Summary of the Invention

[0004] The technical problem to be solved by the present invention is that the existing technology does not systematically assess the harmfulness of mutations and it is difficult to strike a balance between sensitivity and computational efficiency. A gene-phenotype association analysis model is provided. By analyzing the association data of known genes and traits and using machine learning technology and grid search methods to optimize parameters, a prediction model is finally constructed that can effectively evaluate the association between rare gene mutations and traits or diseases.

[0005] In one aspect, the present invention discloses a method for establishing a gene-phenotype association analysis model, comprising the following steps:

[0006] S1: Collect data on traits and genes with known associations to form gene-trait pairs, which are combinations of rare gene mutation types and trait associations;

[0007] S2: Use the following gene rare mutation scoring formula to calculate the score of each gene rare mutation type:

[0008] S=AC×log 10 (1 / AF)×Pr

[0009] Among them, S is the score of rare mutation type of gene, AC is the allele count, AF is the allele frequency, and Pr is the probability of deleteriousness.

[0010] S3: Based on the scores determined in step S2, the linear regression method is used to analyze the correlation between the scores of each rare mutation type in the gene-trait pair in S1 and the traits. The intercept of each rare mutation type in each trait is calculated as the weight reference value of each mutation type. The weight combination of each rare mutation type in different genes is evaluated by the hyperparameter optimization method to obtain the correlation R between the score and the trait. 2 As an evaluation criterion, when R 2 When the value exceeds the preset threshold, it is determined to be a rare mutation type weight combination, which is used as a parameter for the subsequent calculation of the rare mutation load score;

[0011] S4: Apply the scoring formula determined in step S2 and the weight combination determined in step S3 to calculate the rare mutation load score for each gene in each sample according to the following rare mutation load score formula:

[0012]

[0013] Wherein, GBS is the gene mutation burden score; i represents the rare gene mutation type; Beta is the weight of the rare gene mutation type i; S is the score of the rare mutation type i calculated according to the gene rare mutation score formula;

[0014] S5: Analyze the correlation between the rare mutation load score and the phenotype by linear regression or logistic regression method, thereby obtaining the gene-phenotype association analysis model based on the rare mutation score.

[0015] In some of these approaches, in step S1, the rare gene mutation types include at least one of: loss of function (LoF), non-synonymous, synonymous, and other protein coding region mutations (other mutations that may alter protein type and structure) that are not of the above mutation types. It is understood that "other mutations that may alter protein type and structure" refers to mutations that, in addition to the aforementioned loss of function, non-synonymous, and synonymous mutations, are mutations that are considered to have the risk of altering protein type and structure based on common knowledge in the art.

[0016] Furthermore, the loss-of-function mutation includes: at least one of a frameshift mutation, a splicing mutation and a truncation mutation; and / or, the other protein coding region mutations (other mutations that may change the type and structure of the protein) that are not the above mutation types include: at least one of a non-frameshift insertion / deletion mutation, a start codon loss mutation and a stop codon loss mutation.

[0017] In some of the schemes, in step S2, the Pr is used to evaluate non-synonymous mutations, and is obtained in the following manner: using at least three of SIFT, PolyPhen-2 (HDIVandHVAR), LRT, MutationTaster, MutationAssessor, FATHMM, PROVEAN, VEST4, MVP, MPC, PrimateAI3D, DEOGEN2, LIST.S2, CADD, DANN, fathmm-MKL, fathmm-XF, Eigen, GenoCanyon, fitCons and AlphaMissense data analysis software to perform pathogenicity assessment on each rare mutation, using the evaluation value of each data analysis system as input, and using a random forest model to predict the pathogenicity probability prediction value of the rare mutation, which is the Pr value of the non-synonymous mutation; and / or, Pr = 1 for other mutation types (loss-of-function mutations, synonymous mutations, and other protein coding region mutations other than the above mutation types).

[0018] In this study, Pr = 1 is used as a formula parameter rather than a biological pathogenicity assessment. The model only calculates the specific pathogenicity probability (Pr) for non-synonymous mutations using the random forest algorithm. Loss-of-function mutations, synonymous mutations, and other protein-coding region mutations not listed above are mathematically processed using Pr = 1, maintaining their original scores.

[0019] Furthermore, the random forest model was established by the following method: training was performed using non-synonymous mutations annotated as pathogenic, possibly pathogenic, benign, and possibly benign in the public database ClinVar20231112.

[0020] In some of the embodiments, the hyperparameter evaluation method in step S3 includes the following steps:

[0021] (1) Collect the best intercepts of different mutation types in each gene-trait pair calculated by the gene rare mutation score formula as the parameters for input linear regression analysis; preferably, the best intercept refers to the correlation R between the mutation type and the correlation R 2 The intercept at the highest point;

[0022] (2) within the input parameter range, parameter combinations are set with an interval of 0.1±0.02 to obtain multiple weight combination schemes;

[0023] (3) Calculate the correlation R of the known gene-trait pairs determined in step S1 for each weight combination 2 , calculate its average value, and use R 2 The weight combination with the highest average value is used as the weight combination applicable to all phenotypes; preferably, the R 2 The average value is greater than 0.1.

[0024] In some of the above-mentioned solutions, the hyperparameter evaluation method is manual parameter adjustment, grid search, random search or Bayesian optimization.

[0025] In some of the schemes, in step S4, the weight combination is: loss-of-function mutation is 0.10-1.00, non-synonymous mutation is 0.10-1.00, other mutations is 0.00-0.80 and synonymous mutation is 0.00-0.10.

[0026] Furthermore, the weight combination is: loss-of-function mutation is 0.70-0.80, non-synonymous mutation is 0.40-1.00, other mutations are 0.00-0.60 and synonymous mutations are 0.00-0.1.

[0027] Furthermore, the weight combination is loss-of-function mutation=0.8, non-synonymous mutation=0.7, other mutation=0.3, and synonymous mutation=0.

[0028] In some embodiments, the phenotype comprises a trait and / or a disease,

[0029] When the phenotype is a trait, in step S5, the correlation between the rare mutation load score and the trait is analyzed using a linear regression method;

[0030] When the phenotype is a disease, in step S5, the association between the rare mutation burden score and the disease is analyzed using a logistic regression method.

[0031] In some of the schemes, step S5 further includes: obtaining a P value through linear regression or logistic regression method, performing FDR correction on the P value to obtain an FDR value, and setting an FDR threshold to determine the correlation between the analyzed gene and the target trait or disease. When the FDR obtained by gene analysis is greater than the threshold, it is judged as a high-risk gene.

[0032] Furthermore, the FDR correction was performed using the Benjamini-Hochberg method.

[0033] Furthermore, the FDR threshold for traits was less than 0.001, and the FDR threshold for disease-related traits was less than 0.05.

[0034] In some of the solutions, step S5 is specifically as follows:

[0035] (1) Based on the sample ID, gene, and phenotype, an ID-gene dataset is constructed using rare mutation load scores. All genes are traversed, and the mutation load scores of each gene are analyzed with related traits or diseases through linear regression or logistic regression. The total number of people carrying the gene mutation must be greater than 100.

[0036] (2) Determine risk candidate genes based on the corrected FDR value.

[0037] On the other hand, the present invention also discloses a gene-phenotype association analysis model, which is obtained according to the above-mentioned establishment method.

[0038] On the other hand, the present invention also discloses a gene-phenotype association analysis system, comprising: a data storage module, a data analysis module and a result output module.

[0039] Data storage module: used to store the parameters and data for constructing the above-mentioned gene-phenotype association analysis model;

[0040] Data analysis module: used to obtain the whole genome sequencing or whole exome sequencing data of the sample to be evaluated, perform analysis and calculation according to the above-mentioned model establishment method, obtain the gene rare mutation load score of the sample to be evaluated, calculate the correlation between the gene rare mutation load score and the trait or disease, and determine the risk gene according to the preset FDR threshold;

[0041] Result output module: used to output and display the risk gene determination results.

[0042] On the other hand, the present invention also discloses a use of a risk gene in functional assessment, where the function is a trait, the risk gene and its corresponding function are selected from at least one of the following gene-trait pairs:

[0043] (1) Age at Menopause (AAM)-related genes: BANP, CHEK2, GNAZ, PNPLA8, or ZNF518A;

[0044] (2) Apolipoprotein A1 (APOEA)-related genes: ABCA1, ANGPTL3, APOA1, CETP, ICOSLG, LCAT, LIPC, LIPG, LPL, or SCARB1;

[0045] (3) Apolipoprotein B (APOEB)-related genes: ASGR1, LDLR, or PCSK9;

[0046] (4) Body Mass Index (BMI)-related genes: MC4R;

[0047] (5) Calcium metabolism (CAL)-related genes: ALB, ALPL, CASR, CNNM4 or FCGRT;

[0048] (6) Estimated Glomerular Filtration Rate (creatinine-based, EGCR)-related genes: ARG1, CSAG1, DKK1, IL10, MOSPD2, or PEG10;

[0049] (7) Estimated Glomerular Filtration Rate (cystatin-based, EGCY)-related genes: CLDN10, CST3, KCNN4, KRAS or TET2;

[0050] (8) Docosahexaenoic acid (DOA) metabolism-related genes: ANGPTL3, LDLR, or LIPC;

[0051] (9) Glycated hemoglobin A1C (HBA1C)-related genes: EGR3, GCK, JAZF1, PFKM, or RHAG;

[0052] (10) High-density lipoprotein cholesterol (HDL)-related genes: ABCA1, ANGPTL3, APOA1, APOA5, APOC3, CETP, FAM25G, LCAT, LIPC, LIPG, LPL, NR1H3, PPARG, or SCARB1;

[0053] (11) Intraocular pressure (IOP)-related genes: B3GNT7;

[0054] (12) Low-density lipoprotein cholesterol (LDL)-related genes: ANGPTL3, APOB, ASGR1, LDLR, or PCSK9;

[0055] (13) Omega-3 fatty acid metabolism (OTFA)-related genes: ANGPTL3 or LIPC;

[0056] (14) Omega-6 fatty acid metabolism (OSFA)-related genes: ABCA1, ANGPTL3, APOA5, APOB, LCAT, LDLR, LIPC, LIPG, PCSK9, or TIMD4;

[0057] (15) Phosphatidylcholine metabolism (PDCL)-related genes: ABCA1, ANGPTL3, APOA1, CETP, ICOSLG, LCAT, LIPC, LIPG, or SCARB1;

[0058] (16) Phosphoglyceride metabolism (PHG)-related genes: ABCA1, ANGPTL3, APOA1, ICOSLG, LCAT, LIPC, LIPG, or SCARB1;

[0059] (17) Polyunsaturated fatty acid metabolism (PFA)-related genes: ABCA1, ANGPTL3, APOA5, APOB, LCAT, LDLR, LIPC, LIPG, PCSK9 or TIMD4;

[0060] (18) Resting Heart Rate (RHR)-related genes: KCNJ5, PTEN, RGS6 or RRAD;

[0061] (19) Non-HDL cholesterol and non-LDL cholesterol remnant cholesterol (RMNC) related genes: ANGPTL3, APOB, LDLR or PCSK9;

[0062] (20) Sphingomyelin metabolism (SGM)-related genes: ABCA1, ANGPTL3, LCAT, LDLR, LIPC, LPL, PCSK9 or SCARB1;

[0063] (21) Height-related genes: GHRH, HMGA2, IHH, NPR2, NRK, PTPN11, SOCS2, or SPIN4;

[0064] (22) Total cholesterol metabolism (TCH)-related genes: ABCA1, ANGPTL3, APOB, LCAT, LDLR, or PCSK9;

[0065] (23) Total fatty acids metabolism (TFA)-related genes: ANGPTL3, ANGPTL4, APOA5, APOB, G6PC, LIPC, LIPG, PCSK9, or TIMD4;

[0066] (24) Total triglycerides metabolism (TTG)-related genes: ANGPTL3, ANGPTL4, APOA5, APOC3, G6PC, LPL, PPARG, or TIMD4;

[0067] When the function is a disease, the risk gene and its corresponding disease are selected from at least one of the following gene-disease pairs:

[0068] (1) Age-related macular degeneration (AMD)-related genes: PRPH2, ZNF664, or CFH;

[0069] (2) Alzheimer's disease (AD)-related genes: PSEN1, SORL1, or IZUMO1R;

[0070] (3) Asthma (AST)-related genes: FLG or ATP2A3;

[0071] (4) Atrial fibrillation (AF)-related genes: PLEC, PKP2, KDM5B, SLC9A9, LMNA, MUL1, or CTNNA3;

[0072] (5) Bipolar Disorder (BD)-related genes: ELFN1 or VSIG8;

[0073] (6) Colorectal cancer (CRC)-related genes: MLH1, MSH2 or MSH6;

[0074] (7) Breast Cancer (BC)-related genes: BRCA2, ATMIN, BRCA1, PALB2, ATM, or BARD1;

[0075] (8) Cardiovascular disease (CVD)-related genes: LDLR;

[0076] (9) Crohn's disease (CD)-related genes: KLHDC4;

[0077] (10) Epithelial Ovarian Cancer (EOC)-related genes: BRCA2, BRCA1, RAD51D or RAD51C;

[0078] (11) Hypertension (HT)-related genes: NOS3, REN or FES;

[0079] (12) Melanoma (MEL)-related genes: CDKN2A, POGK or DCLRE1B;

[0080] (13) Osteoporosis (OP)-related genes: LRP5 or WNT1;

[0081] (14) Parkinson's disease (PD)-related genes: MAP2K7;

[0082] (15) Primary open-angle glaucoma (POAG)-related genes: LTBP2;

[0083] (16) Prostate cancer (PC) related genes: ATM;

[0084] (17) Type 2 diabetes (T2D)-related genes: GCK, GIGYF1, C19orf57, HNF1A, MIB1, HMGXB4, IGF1R, or HNF4A;

[0085] (18) Venous Thromboembolism (VTE)-related genes: REC8, PROC, PROS1 or JAK2.

[0086] On the basis of conforming to the common sense in this field, the above-mentioned preferred conditions can be arbitrarily combined to obtain the preferred embodiments of the present invention.

[0087] The reagents and raw materials used in the present invention are commercially available.

[0088] The positive progress of the present invention is that it is a new method for screening candidate genes for diseases or traits by calculating the rare mutation load score. Through the effect examples, it can be seen that most of the risk genes screened out have been reported, which shows the reliability of this method. In addition, when comparing the disease association results with traditional methods, some new associated genes are obtained, and most of them have been reported. Therefore, it can be seen that this method can be used as an effective supplement to traditional methods in disease-related pathogenic or risk genes. Through this model, candidate risk genes can be obtained through association analysis by evaluating mutations through gene mutation types and prediction software, and new pathogenic genes or risk genes for diseases or traits can be discovered. Compared with the traditional method gene-base collapsing method, it has better reproducibility and complementarity. This model can play a complementary role in finding candidate risk genes associated with new traits or unknown diseases. BRIEF DESCRIPTION OF THE DRAWINGS

[0089] Figure 1 The process of determining the best formula.

[0090] Figure 2 The correlation and screening process of each formula for different gene-trait pairs.

[0091] Figure 3 Compute correlations for gene-trait pairs for hyperparameters.

[0092] Figure 4 It is the optimal weight combination of the four major mutations.

[0093] Figure 5 Genes associated with continuous traits and their reporting status.

[0094] Figure 6 Reports on disease-related genes and their comparison with traditional collapsing methods. DETAILED DESCRIPTION

[0095] The present invention is further illustrated by way of examples below, but the present invention is not limited to the scope of the examples. Experimental methods in the following examples where specific conditions are not specified were performed according to conventional methods and conditions, or selected according to the product specifications.

[0096] Example 1 Construction of gene-trait association prediction method

[0097] (1) Traits and genes with known associations

[0098] To calculate the rare mutation score association formula for a trait with the most appropriate gene, we first need to find known gene-trait pairs with associations. Referring to the published article (Fiziev P, McRae J, Ulirsch JC, et al. Rare penetrant mutations confer severe risk of common diseases. Preprint.medRxiv.2023; 2023.05.01.23289356. Published 2023May 8. doi:10.1101 / 2023.05.01.23289356), the rare mutations and associated traits of the following genes were used as known gene rare mutation and trait association combinations, as shown in Table 1.

[0099] Table 1. Rare mutations in known associated genes and trait combinations

[0100]

[0101]

[0102] (2) Determination of the optimal formula

[0103] Rare gene mutations are divided into four categories: loss of function mutations (LoF, mainly including frameshift mutations, splicing mutations and stopgain mutations), non-synonymous mutations, synonymous mutations and other protein coding region mutations that are not of the above mutation types (including non-frameshift insertion and deletion, start codon loss and other mutations), see Table 2.

[0104] Table 2. Major and minor categories of mutation types

[0105]

[0106] Based on the type of rare mutation of an individual's gene, the number of alleles (AC), the allele mutation frequency (AF), and the harmful assessment score of non-synonymous mutations (Pr, Probability score for deleteriousness, ranging from 0 to 1), seven formulas were designed to evaluate the rare mutation scores of different mutation types. The optimal formula was obtained by calculating the RMSE (Root Mean Square Error) through linear regression.

[0107] Each mutation type is scored using the following seven formulas, and then the scores of each mutation type are added together to obtain the rare mutation score of the gene for that individual. Figure 1 The data were validated by screening 500,000 people from the UK biobank for data with the gene mutations and phenotypes listed in Table 1. The screening criteria were: for each specific trait-gene pair, the subject must carry the gene mutation and have a test record for the trait. Table 3 shows the amount of data used in this example.

[0108] Table 3. Total number of people and total number of mutations for each gene trait of the present invention

[0109]

[0110]

[0111] First is formula (1), the specific formula is as follows, only considering the influence of the number of allele mutations (AC). From Table 4, and Figure 2 From the A results, we can see that the use of AC can preliminarily show the correlation between gene mutations and diseases, but the effect is average in comparison.

[0112] S=AC (1)

[0113] Then, according to each gene-trait pair in Table 1, with the sample individual ID, four major mutations, and traits as the horizontal coordinates, and each person's AC as the main parameter of the four major mutations, linear regression (lm) was used to calculate the regression parameters of different gene-trait pairs with the minimum root mean square error RMSE, including Rsquared (R 2 , determination coefficient) and the intercepts of the four major categories of mutations (i.e., the weights of the preliminary assessment) are shown in Table 4.

[0114] Table 4. Regression intercept and correlation evaluation of formula (1)

[0115]

[0116]

[0117] Using the same method, the same parameters of formulas (2)-(6) were evaluated respectively. The specific formulas are as follows. All of them are evaluated based on the allele mutation frequency, but different conversion methods are used: reciprocal, reciprocal square root, ln, log2 and log10. The specific evaluation results are shown in Tables 5-9.

[0118] S=AC×1 / AF (2)

[0119]

[0120] S=AC×ln(1 / AF) (4)

[0121] S=AC×log2(1 / AF) (5)

[0122] S=AC×log 10 (1 / AF) (6)

[0123] Table 5. Regression intercept and correlation evaluation of formula (2)

[0124]

[0125]

[0126] Table 6. Regression intercept and correlation estimates for formula (3)

[0127]

[0128]

[0129] Table 7. Regression intercept and correlation estimates for formula (4)

[0130]

[0131]

[0132] Table 8. Regression intercept and correlation estimates for formula (5)

[0133]

[0134]

[0135] Table 9. Regression intercept and correlation estimates for formula (6)

[0136]

[0137]

[0138] By R 2, Figure 2 The results of A show that the reciprocal and the reciprocal square root are not as effective as the log permutation, and the effects of three different bases: 2, e, and 10 are similar. Since the allele frequency shows stable performance and good model fitting effect when the base is 10 as the statistical standard, and it has been verified by the inventors' experiments, log10 is used as the permutation function of AF.

[0139] Finally, the present invention further evaluates the impact of the probability of non-synonymous mutations on the entire gene-trait association. Since LoF has been proven to be able to directly associate rare gene mutations with traits in the form of AC, synonymous mutations have almost been proven to have weak correlations. Other mutations that may change protein type and structure have relatively low frequencies in the genome (see Table 2), and their individual contributions are difficult to accurately quantify. The present invention uses AC as its scoring benchmark value. For non-synonymous mutations, AC×Pr is used as its specific value in this embodiment. Regarding the evaluation of the harmful probability of non-synonymous mutations, reference was made to the machine learning-based method proposed by Shen Maoting et al. (2024) (Shen Maoting, Lin Junwei, Fan Xijie, et al. Prediction of pathogenicity of non-synonymous mutations based on machine learning and analysis of their characteristic importance [J]. Journal of Guangzhou Medical University, 2024, 52(5): 1-9. DOI: 10.3969 / j.issn.2095-9664.2024.05.01.). This method integrates over 20 different prediction methods to assess the pathogenicity of mutations on a probabilistic scale of 0-1. Based on this probability, this example assigns a weight of 10 times the probability of being harmful, i.e., 1-10. These 20 prediction methods include: SIFT, PolyPhen-2 (HDIV and HVAR), LRT, MutationTaster, MutationAssessor, FATHMM, PROVEAN, VEST4, MVP, MPC, PrimateAI3D, DEOGEN2, LIST.S2, CADD, DANN, fathmm-MKL, fathmm-XF, Eigen, GenoCanyon, fitCons, and AlphaMissense.

[0140] In this embodiment, the random forest model was established by the following method: training was performed using 1,344,221 non-synonymous mutations annotated as pathogenic, likely pathogenic, benign, and likely benign in the public database ClinVar20231112 (https: / / www.ncbi.nlm.nih.gov / clinvar / ).

[0141] The random forest model used was trained using conventional methods in the field, such as the "caret" package, with method = "rf", number = 10, and repeats = 5. The model was validated using a confusion matrix. The model achieved an internal validation accuracy of 0.9265 and an AUC of 0.974.

[0142] Therefore, according to formula (7), the same method is used to obtain its regression parameters. The specific formula is as follows, and the evaluation results are shown in Table 10.

[0143] S=AC×log 10 (1 / AF)×Pr (7)

[0144] Table 10. Regression intercept and correlation estimates for formula (7)

[0145]

[0146]

[0147] In summary, the present invention evaluates the effects of rare mutations AC, AF and non-synonymous mutation deleterious frequency scores (Pr) on gene-trait pairs, and finally determines formula (7) as the optimal formula, see Figure 2 A.

[0148] (3) Hyperparameter evaluation of the combination of mutation type weights and determination of the optimal weight

[0149] After determining the optimal formula for rare mutation scoring, linear regression was used to obtain the intercepts of different mutation types in each trait based on the correlation between the scores of different mutation types and different traits, which served as the initial weight reference for each mutation type. Since the intercepts of mutation types obtained for different traits were inconsistent, in order to find a universal mutation type weight, the optimal intercept (i.e., the highest R) of each gene-trait pair in different mutation types calculated according to the optimal formula was used. 2 The intercept of the value is used as one of the input parameters; then, within these ranges, with a spacing of 0.1, a grid search is used for evaluation, resulting in a total of 14,658 combinations. Then, the weighted combinations applicable to all traits are screened out, and the weighted scores and the correlation R of the traits are used to calculate the weighted scores. 2 as evaluation criteria.

[0150] The weighted score Sw can be calculated by the following formula:

[0151]

[0152] Finally, among all suitable combinations, the best weight is selected based on the mean and median.

[0153] By calculating 14658 combinations, the correlation R of 16 gene-trait pairs from Table 1 corresponding to each combination was obtained. 2 , and then calculate its mean and variance and other parameters, and according to R 2 The average value is greater than 0.1 to screen out the weight combinations with better performance (234 in total), see Figure 3 and Table 11.

[0154] Table 11R 2 Weight combinations with an average value greater than 0.1 (the top 50 weight combinations with better performance are shown as examples)

[0155]

[0156]

[0157] Finally, the best weight combination of the four major mutation types was screened out by the mean and median, namely LoF = 0.8, non-synonymous mutation = 0.70, other mutations = 0.30 and synonymous mutation = 0.00, see Figure 4 and Table 11.

[0158] Based on the above steps and analysis, the optimal formula for rare mutation scores and the optimal weights for the four major mutation types were finally determined. The following gene rare mutation load score formula was used to calculate the rare mutation load score for each person and each gene. The risk candidate genes for traits or diseases were screened through association analysis between phenotype and gene rare mutation load score, i.e., linear regression of traits and logistic regression of diseases. Figure 2 B.

[0159]

[0160] #Note: GBS = Gene Burden Score; i represents four mutation types, namely LoF mutation, nonsynonymous mutation, synonymous mutation, and other mutations; Beta is the optimal weight of mutation type (LoF = 0.8, nonsynonymous mutation = 0.7, other mutation = 0.3, synonymous mutation = 0); AC and AF are the number and frequency of allele mutations, respectively; Pr: probability of deleteriousness of nonsynonymous mutations, Pr = 1 for other mutation types.

[0161] (4) Correlating the deleterious mutation load score with the disease / trait through regression methods to screen for pathogenic or risk genes

[0162] The rare mutation load score of each gene in each sample is calculated according to the gene rare mutation load score formula, and then association analysis is performed through linear regression (for traits of continuous variables) or logistic regression (for diagnosis of disease) based on specific traits (such as low-density lipoprotein, etc.) or diseases (such as type 2 diabetes, etc.).

[0163] The specific steps are as follows: First, based on the sample ID, gene, and trait, a data set of ID and gene is formed using the rare mutation load score; then, by scanning all genes, the mutation load score of the gene is analyzed with the related trait or disease one by one through linear regression or logistic regression analysis, requiring that the total number of people carrying the gene mutation is greater than 100; the p-value obtained from the regression analysis is FDR corrected, for example, using the Benjamini-Hochberg method; finally, the risk candidate genes are determined based on the corrected FDR value, and the FDR threshold is set to screen the risk candidate genes for related traits or diseases, where the FDR for continuous traits is <0.001 and the FDR for disease-related traits is <0.05.

[0164] Example 2 A gene-phenotype association analysis model

[0165] This embodiment provides a specific implementation scheme of a gene-trait association analysis model based on rare mutation scores, which is implemented according to the following steps (1)-(5):

[0166] (1) Data collection and processing

[0167] Collect data on traits and genes with known associations as rare gene mutation types and trait association combinations, and preprocess these data, including quality control, outlier processing, and standardization, in preparation for subsequent analysis.

[0168] (2) Calculation of rare gene mutation scores

[0169] The collected mutations were annotated into four categories: loss-of-function mutations (including frameshift mutations, splicing mutations and truncation mutations), non-synonymous mutations, synonymous mutations and other protein coding region mutations not of the above mutation types (including non-frameshift insertion-deletion mutations, start codon loss mutations and stop codon loss mutations).

[0170] The following gene rare mutation scoring formula is used to calculate the score of each mutation:

[0171] S=AC×log 10 (1 / AF)×Pr

[0172] Among them: S is the score of rare mutation type of gene, AC is the number of alleles, AF is the allele mutation frequency, and Pr is the harmful probability score, ranging from 0 to 1.

[0173] In actual calculations: For non-synonymous mutations, the system uses a random forest model to predict the Pr value, which uses the evaluation results of software such as SIFT, PolyPhen-2, and MutationTaster as input (see Example 1 for details).

[0174] For other mutations, considering their small effects, the Pr value is set to 1

[0175] (3) Hyperparameter optimization to determine weight combination

[0176] The hyperparameter method (manual parameter adjustment) is used to evaluate the weight combination of different mutation types. The optimization process is as follows:

[0177] Collect the optimal intercepts of each gene-trait pair under different mutation types calculated by the gene rare mutation scoring formula as input parameters; within these parameter ranges, set parameter combinations with an interval of 0.1 to obtain multiple weight combination schemes; calculate the correlation R of the known gene-trait pair determined in step (1) corresponding to each weight combination 2 , calculate its mean and variance and other parameters, according to R 2 The average value was greater than 0.1, which was used to select the weight combination that was suitable for all traits. Based on the average and median, the optimal weight combination was determined to be: LoF = 0.8, non-synonymous mutation = 0.7, other mutations = 0.3, and synonymous mutation = 0.

[0178] (4) Calculation of rare mutation load score

[0179] Using the scoring formula in step (2) and the optimal weight combination determined in step (3), the rare mutation load score of each gene in each sample is calculated according to the following gene rare mutation load score formula:

[0180]

[0181] Where: GBS is the gene mutation burden score, i represents the rare gene mutation type, Beta_ is the weight of the rare gene mutation type i (the optimal weight determined in step (3)), and S is the score of the rare mutation type i calculated according to the scoring formula.

[0182] (5) Establishment of correlation analysis model

[0183] Based on the relationship between gene mutation load score and traits / diseases, an association analysis model was established: for specific traits (such as low-density lipoprotein, etc.), a linear regression method was used with GBS as the independent variable and the trait value as the dependent variable; for disease status (such as type 2 diabetes, etc.), a logistic regression method was used with GBS as the independent variable and the disease status (0 / 1) as the dependent variable.

[0184] Based on the sample ID, gene, and trait, an ID-gene dataset is constructed using rare mutation load scores. All genes are traversed, and the mutation load scores of each gene are associated with related traits or diseases. The total number of people carrying the gene mutation must be greater than 100.

[0185] The p-values ​​obtained from the regression analysis were FDR-corrected, and risk candidate genes were identified based on the corrected FDR values. The FDR threshold for continuous traits was less than 0.001, and the FDR threshold for disease-related traits was less than 0.05.

[0186] Through the above steps (1)-(5), the system successfully established a gene-trait association analysis model based on rare mutation scores.

[0187] Effect Example 1 Validation of a gene-trait association analysis model based on rare mutation scores

[0188] (1) Screening of trait-associated risk genes using deleterious mutation load scores

[0189] In order to verify that the harmful mutation load score in the gene-trait association analysis model based on rare mutation score in Example 2 can effectively screen risk genes associated with traits, data from 500,000 people in the UK biobank with gene mutations in Table 12 and 25 different trait phenotypes were screened for verification. The screening criteria were: for each specific trait-gene pair, the subject must carry the mutation of the gene and have a detection record for the trait. Linear regression analysis was performed using the model in Example 2. The mutation load score constructed using the most suitable mutation type weight and mutation factors (AF, log10, 1 / AF and Pr) was used to obtain the credibility of the associated gene by associating with the trait. 149 genes associated with 24 traits were screened using FDR<0.001, of which 128 gene-traits were reported, accounting for 85.91%, indicating that the model has a certain degree of reliability ( Figure 5 , Table 12).

[0190] Table 12 Association results of rare mutation scores of trait-related genes

[0191]

[0192]

[0193]

[0194]

[0195]

[0196]

[0197]

[0198]

[0199] (2) Screening of disease-related genes using deleterious mutation load scores

[0200] In order to verify that the deleterious mutation load score in the gene-trait association analysis model based on rare mutation score in Example 2 can effectively screen risk genes associated with diseases, data from 500,000 people in the UK biobank with gene mutations in Table 13 and 28 different disease phenotypes were screened for validation. The screening criteria were: for each specific trait-gene pair, the subject must carry the mutation in the gene and have a test record for the disease. Logistic regression analysis was performed using the model in Example 2. Using an FDR < 0.05, 55 genes associated with 18 diseases were screened, of which 45 disease-gene pairs were reported, accounting for 81.82% ( Figure 6 A, Table 13). This result also illustrates the credibility of the association analysis based on the gene mutation load score. In addition, in order to further verify its effectiveness, the traditional gene-based collapsing method (Fiziev P, McRae J, Ulirsch JC, et al. Rare penetrant mutations confer severe risk of common diseases. Preprint.medRxiv.2023; 2023.05.01.23289356. Published2023May 8. doi: 10.1101 / 2023.05.01.23289356) was used to verify and compare the same group of data screened from the above-mentioned UK biobank 500,000 people. Collapsing screened a total of 109 associated genes for 22 diseases, of which 83 disease-gene pairs have been reported, accounting for 76.15% ( Figure 6 A, Table 14), the reported rate is similar to that of the method, but the total number is larger. It was found that there are 33 groups of disease-gene pairs that are overlapping in the two different methods, all of which are reported disease-gene pairs ( Figure 6 B). Compared with the traditional collapsing method, the method of

[15] should be an effective supplement to this method and be useful in finding more new disease-associated genes.

[0201] Table 13 Association results of rare mutation scores of disease-related genes

[0202]

[0203]

[0204]

[0205] Table 14 Collapsing analysis results

[0206]

[0207]

[0208]

[0209] Although the above describes specific embodiments of the present invention, it should be understood by those skilled in the art that these are merely illustrative and that various changes or modifications may be made to these embodiments without departing from the principles and essence of the present invention. Therefore, the scope of protection of the present invention is defined by the appended claims.

Claims

1. A method for establishing a gene-phenotype association analysis model, characterized in that: The following steps are involved: S1: Collect data on traits and genes with known associations to form gene-trait pairs, which are combinations of rare gene mutation types and trait associations; S2: Use the following gene rare mutation scoring formula to calculate the score of each gene rare mutation type: S=AC×log 10 (1 / AF)×Pr Among them, S is the score of rare mutation type of gene, AC is the number of alleles, AF is the allele mutation frequency, and Pr is the probability of harmfulness; S3: Based on the scores determined in step S2, the linear regression method is used to analyze the correlation between the scores of each rare mutation type in the gene-trait pair in S1 and the traits. The intercept of each rare mutation type in each trait is calculated as the weight reference value of each mutation type. The weight combination of each rare mutation type in different genes is evaluated by the hyperparameter optimization method to obtain the correlation R between the score and the trait. 2 As an evaluation criterion, when R 2 When the value exceeds the preset threshold, it is determined to be a rare mutation type weight combination, which is used as a parameter for the subsequent calculation of the rare mutation load score; S4: Apply the scoring formula determined in step S2 and the weight combination determined in step S3 to calculate the rare mutation load score for each gene in each sample according to the following rare mutation load score formula: Wherein, GBS is the gene mutation burden score; i represents the rare gene mutation type; Beta is the weight of the rare gene mutation type i; S is the score of the rare mutation type i calculated according to the gene rare mutation score formula; S5: Analyze the correlation between the rare mutation load score and the phenotype by linear regression or logistic regression method, thereby obtaining the gene-phenotype association analysis model based on the rare mutation score.

2. The establishment method according to claim 1, characterized in that In step S1, the rare gene mutation types include: at least one of loss-of-function mutations, non-synonymous mutations, synonymous mutations, and other protein coding region mutations that are not of the above mutation types; Preferably, the loss-of-function mutation includes at least one of a frameshift mutation, a splicing mutation, and a truncation mutation; and / or, the other protein coding region mutations other than the above mutation types include at least one of a non-frameshift insertion / deletion mutation, a start codon loss mutation, and a stop codon loss mutation.

3. The establishment method according to claim 2, characterized in that: In step S2, the Pr is used to evaluate non-synonymous mutations and is obtained by: using at least three of SIFT, PolyPhen-2 (HDIVandHVAR), LRT, MutationTaster, MutationAssessor, FATHMM, PROVEAN, VEST4, MVP, MPC, PrimateAI3D, DEOGEN2, LIST.S2, CADD, DANN, fathmm-MKL, fathmm-XF, Eigen, GenoCanyon, fitCons and AlphaMissense data analysis software to perform pathogenicity assessment on each rare mutation, using the evaluation value of each data analysis system as input, and using a random forest model to predict the pathogenicity probability prediction value of the rare mutation, which is the Pr value of the non-synonymous mutation; and / or, Pr = 1 for other mutation types; Preferably, the random forest model is established by the following method: training is performed using non-synonymous mutations annotated as pathogenic, possibly pathogenic, benign, and possibly benign in the public database ClinVar20231112.

4. The establishment method according to claim 1, wherein: The hyperparameter evaluation method in step S3 comprises the following steps: (1) Collect the best intercepts of different mutation types in each gene-trait pair calculated by the gene rare mutation score formula as the parameters for input linear regression analysis; preferably, the best intercept refers to the correlation R between the mutation type and the correlation R 2 The intercept at the highest point; (2) within the input parameter range, parameter combinations are set with an interval of 0.1±0.02 to obtain multiple weight combination schemes; (3) Calculate the correlation R of the known gene-trait pairs determined in step S1 for each weight combination 2 , calculate its average value, and use R 2 The weight combination with the highest average value is used as the weight combination applicable to all phenotypes; preferably, the R 2 The average value is greater than 0.

1.

5. The establishment method according to claim 4, characterized in that: In step S4, the weight combination is: loss-of-function mutation is 0.10-1.00, non-synonymous mutation is 0.10-1.00, other mutations are 0.00-0.80, and synonymous mutations are 0.00-0.10; Preferably, the weight combination is: loss-of-function mutation is 0.70-0.80, non-synonymous mutation is 0.40-1.00, other mutations are 0.00-0.60 and synonymous mutations are 0.00-0.1; More preferably, the weight combination is loss-of-function mutation=0.8, non-synonymous mutation=0.7, other mutation=0.3, and synonymous mutation=0.

6. The establishment method according to claim 1, wherein: The phenotype includes traits and / or diseases, When the phenotype is a trait, in step S5, the correlation between the rare mutation load score and the trait is analyzed using a linear regression method; When the phenotype is a disease, in step S5, the association between the rare mutation burden score and the disease is analyzed using a logistic regression method.

7. The establishment method according to claim 6, characterized in that: The step S5 further comprises: obtaining a P value by linear regression or logistic regression method, performing FDR correction on the P value to obtain an FDR value, and setting an FDR threshold to determine the correlation between the analyzed gene and the target trait or disease. When the FDR obtained by gene analysis is greater than the threshold, it is determined to be a high-risk gene; Preferably, the FDR correction is the Benjamini-Hochberg method; More preferably, the FDR threshold for traits is less than 0.001, and the FDR threshold for disease-related traits is less than 0.

05.

8. The establishment method according to claim 7, characterized in that: The step S5 is specifically as follows: (1) Based on the sample ID, gene, and phenotype, an ID-gene dataset is constructed using rare mutation load scores. All genes are traversed, and the mutation load scores of each gene are analyzed with related traits or diseases through linear regression or logistic regression. The total number of people carrying the gene mutation must be greater than 100. (2) Determine risk candidate genes based on the corrected FDR value.

9. A gene-phenotype association analysis model, characterized in that: Obtained by the establishment method according to any one of claims 1 to 8.

10. A gene-phenotype association analysis system, characterized in that: include: Data storage module, data analysis module and result output module, Data storage module: used for storing parameters and data for constructing the gene-phenotype association analysis model according to claim 9; A data analysis module is used to obtain whole genome sequencing or whole exome sequencing data of the sample to be evaluated, perform analysis and calculation according to the model establishment method described in any one of claims 1 to 8, obtain the gene rare mutation load score of the sample to be evaluated, calculate the correlation between the gene rare mutation load score and the trait or disease, and determine the risk gene according to a preset FDR threshold; Result output module: used to output and display the risk gene determination results.

11. A use of a risk gene in functional assessment, characterized in that: When the function is a trait, the risk gene and its corresponding function are selected from at least one of the following gene-trait pairs: (1) Menopausal age-related genes: BANP, CHEK2, GNAZ, PNPLA8, or ZNF518A; (2) apolipoprotein A1-related genes: ABCA1, ANGPTL3, APOA1, CETP, ICOSLG, LCAT, LIPC, LIPG, LPL, or SCARB1; (3) Apolipoprotein B-related genes: ASGR1, LDLR, or PCSK9; (4) BMI-related genes: MC4R; (5) calcium metabolism-related genes: ALB, ALPL, CASR, CNNM4, or FCGRT; (6) Genes related to glomerular filtration rate estimation based on creatinine: ARG1, CSAG1, DKK1, IL10, MOSPD2, or PEG10; (7) Genes related to glomerular filtration rate estimation based on cystatin: CLDN10, CST3, KCNN4, KRAS or TET2; (8) Docosahexaenoic acid metabolism-related genes: ANGPTL3, LDLR, or LIPC; (9) Glycated hemoglobin-related genes: EGR3, GCK, JAZF1, PFKM or RHAG; (10) high-density lipoprotein cholesterol-related genes: ABCA1, ANGPTL3, APOA1, APOA5, APOC3, CETP, FAM25G, LCAT, LIPC, LIPG, LPL, NR1H3, PPARG, or SCARB1; (11) Intraocular pressure-related genes: B3GNT7; (12) LDL cholesterol-related genes: ANGPTL3, APOB, ASGR1, LDLR, or PCSK9; (13) ω-3 fatty acid metabolism-related genes: ANGPTL3 or LIPC; (14) ω-6 fatty acid metabolism-related genes: ABCA1, ANGPTL3, APOA5, APOB, LCAT, LDLR, LIPC, LIPG, PCSK9, or TIMD4; (15) Phosphatidylcholine metabolism-related genes: ABCA1, ANGPTL3, APOA1, CETP, ICOSLG, LCAT, LIPC, LIPG, or SCARB1; (16) Genes related to phosphoglyceride metabolism: ABCA1, ANGPTL3, APOA1, ICOSLG, LCAT, LIPC, LIPG, or SCARB1; (17) Polyunsaturated fatty acid metabolism-related genes: ABCA1, ANGPTL3, APOA5, APOB, LCAT, LDLR, LIPC, LIPG, PCSK9, or TIMD4; (18) Resting heart rate-related genes: KCNJ5, PTEN, RGS6, or RRAD; (19) Non-HDL cholesterol and non-LDL cholesterol remnant cholesterol-related genes: ANGPTL3, APOB, LDLR, or PCSK9; (20) Sphingomyelin metabolism-related genes: ABCA1, ANGPTL3, LCAT, LDLR, LIPC, LPL, PCSK9, or SCARB1; (21) Height-related genes: GHRH, HMGA2, IHH, NPR2, NRK, PTPN11, SOCS2, or SPIN4; (22) Total cholesterol metabolism-related genes: ABCA1, ANGPTL3, APOB, LCAT, LDLR, or PCSK9; (23) total fatty acid metabolism-related genes: ANGPTL3, ANGPTL4, APOA5, APOB, G6PC, LIPC, LIPG, PCSK9, or TIMD4; (24) Total triglyceride metabolism-related genes: ANGPTL3, ANGPTL4, APOA5, APOC3, G6PC, LPL, PPARG, or TIMD4; when the function is a disease, the risk gene and its corresponding disease are selected from at least one of the following gene-disease pairs: (1) Age-related macular degeneration-related genes: PRPH2, ZNF664, or CFH; (2) Alzheimer's disease-related genes: PSEN1, SORL1, or IZUMO1R; (3) Asthma-related genes: FLG or ATP2A3; (4) atrial fibrillation-related genes: PLEC, PKP2, KDM5B, SLC9A9, LMNA, MUL1, or CTNNA3; (5) Bipolar disorder-related genes: ELFN1 or VSIG8; (6) Colorectal cancer-related genes: MLH1, MSH2, or MSH6; (7) Breast cancer-related genes: BRCA2, ATMIN, BRCA1, PALB2, ATM, or BARD1; (8) Cardiovascular disease-related genes: LDLR; (9) Crohn's disease-related genes: KLHDC4; (10) Epithelial ovarian cancer-related genes: BRCA2, BRCA1, RAD51D or RAD51C; (11) Hypertension-related genes: NOS3, REN or FES; (12) Melanoma-related genes: CDKN2A, POGK, or DCLRE1B; (13) Osteoporosis-related genes: LRP5 or WNT1; (14) Parkinson's disease-related genes: MAP2K7; (15) Primary open-angle glaucoma-related genes: LTBP2; (16) Prostate cancer-related genes: ATM; (17) Type 2 diabetes-related genes: GCK, GIGYF1, C19orf57, HNF1A, MIB1, HMGXB4, IGF1R, or HNF4A; (18) Venous thromboembolic disease-related genes: REC8, PROC, PROS1 or JAK2.

Images