A gene-phenotype association analysis model and a method for establishing and applying the same

Abstract

The application discloses a gene-phenotype correlation analysis model and a method and application thereof, and belongs to the technical field of biological medicine. The model establishment method comprises the following steps: S1, collecting known trait-gene data to form a gene-trait pair; S2, calculating the rare mutation type score of each gene by using a formula; S3, analyzing the correlation between the mutation score and the trait by linear regression, calculating the weight of each mutation type, and optimizing the weight combination, so that the correlation R 2 is taken as the evaluation standard; S4, calculating the rare mutation load score of a sample gene according to the scoring formula and the optimized weight; S5, analyzing the correlation between the mutation load score and the phenotype by a regression method, and constructing a gene-phenotype correlation analysis model. Compared with a traditional gene-base collapsing method, the model has good reproducibility and complementarity, and can be used for discovering candidate risk genes of new traits or unknown diseases.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, and in particular to a gene-phenotype association analysis model, its establishment method, and its application. Background Technology

[0002] Rare mutations are key factors influencing diseases and traits. Gene association analysis based on rare mutations can screen for candidate genes associated with diseases or traits, and is a common method for identifying disease risk or pathogenic genes. Assessing the pathogenicity or harmfulness of rare mutations is crucial for finding associated genes. Although numerous prediction software programs exist that can predict different mutations and have good assessment results, these results cannot always perfectly match experimental validation. Furthermore, due to the sheer volume of data—often tens of millions of mutations—it is impractical to conduct experimental evaluations at every locus.

[0003] Traditional gene-based collapsing methods, based on rare mutations, only consider different levels of mutation type and lack a systematic assessment of the harmfulness of mutations. Other types of gene-based analysis methods, such as SAIGE-GENE and SKAT-o, while possessing very high specificity, lack sensitivity and may miss some candidate genes. Furthermore, these methods require integrating rare mutations into plink files, which, due to the sparsity of rare mutations in the population, often result in large files, increasing workload, time, and storage costs. In addition, current research has not clearly defined which combinations of mutation type weights can achieve similar analytical results, or which weight combination ultimately achieves the best results. Moreover, comparative studies between various gene analysis methods are insufficient; for example, how to compare newly developed methods with gene-based collapsing, which method is better, and whether they are complementary assessments are also issues that need to be addressed. Summary of the Invention

[0004] The technical problem to be solved by this invention is that the assessment of the harmfulness of mutations is not systematic in the existing technology and that it is difficult to balance sensitivity and computational efficiency. This invention provides a gene-phenotype association analysis model. By analyzing the association data of known genes and traits, the model uses machine learning technology and grid search method to optimize parameters and finally constructs a predictive model that can effectively assess the association between rare gene mutations and traits or diseases.

[0005] On the one hand, this invention discloses a method for establishing a gene-phenotype association analysis model, comprising the following steps:

[0006] S1: Collect known related trait and gene data to form gene-trait pairs, which serve as combinations of rare gene mutation types and trait associations;

[0007] S2: Calculate the score for each rare mutation type of gene using the following rare mutation scoring formula:

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

[0009] Where S is the score for rare gene mutation types, AC is the number of alleles, AF is the allele frequency, and Pr is the probability score for deletion.

[0010] S3: Based on the scores determined in step S2, use linear regression to analyze the correlation between the scores and traits of each rare mutation type in the gene-trait pairs in S1. Calculate the intercept for each rare mutation type in each trait, using it as a weight reference value for each mutation type. Then, evaluate the weight combinations of each rare mutation type in different genes using hyperparameter optimization methods, using the correlation R between the score and the trait. 2 As an evaluation criterion, when R 2 When the value exceeds a preset threshold, it is identified as a rare mutation type weight combination and used as a parameter for 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, and calculate the rare mutation load score for each gene in each sample according to the following rare mutation load score formula:

[0012]

[0013] Where GBS is the gene mutation load score; i represents the rare gene mutation type; Beta is the weight of rare gene mutation type i; and S is the score of rare mutation type i calculated according to the rare gene mutation scoring formula.

[0014] S5: Analyze the association between rare mutation load score and phenotype using linear regression or logistic regression methods to obtain the gene-phenotype association analysis model based on rare mutation score.

[0015] In some of these schemes, in step S1, the rare gene mutation type includes at least one of the following: loss-of-function mutations (LoF), non-synonymous mutations, synonymous mutations, and other protein-coding region mutations (other mutations that may alter protein type and structure). It is understood that "other mutations that may alter protein type and structure" refers to mutations other than the aforementioned loss-of-function, non-synonymous, and synonymous mutations, which, based on common knowledge in the art, are considered to pose a risk of altering protein type and structure.

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

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

[0018] In this invention, 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, while loss-of-function mutations, synonymous mutations, and other protein-coding region mutations are treated by setting Pr=1 as a mathematical processing method to keep their original scores unchanged.

[0019] Furthermore, the random forest model is established by training it with non-synonymous mutations annotated as pathogenic, potentially pathogenic, benign, and potentially benign in the public database ClinVar20231112.

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

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

[0022] (2) Within the range of the input parameters, multiple weight combination schemes are obtained by setting the parameter combination with an interval of 0.1±0.02.

[0023] (3) Calculate the correlation R of the known gene-trait pairs determined in step S1 for each weight combination. 2 Calculate their average value, expressed as 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 schemes, the hyperparameter evaluation method is manual parameter tuning, grid search, random search, or Bayesian optimization.

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

[0026] Furthermore, the weight combination is as follows: loss-of-function mutations are 0.70-0.80, non-synonymous mutations are 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 as follows: missing-function mutation = 0.8, non-synonymous mutation = 0.7, other mutation = 0.3, and synonymous mutation = 0.

[0028] In some of these schemes, the phenotype includes traits and / or diseases.

[0029] When the phenotype is a trait, in step S5, the association between the rare mutation load fraction 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 logistic regression.

[0031] In some of these schemes, step S5 further includes: obtaining a P-value through linear regression or logistic regression, correcting the P-value for FDR to obtain an FDR value, and setting an FDR threshold to determine the degree of correlation between the analyzed gene and the target trait or disease. If the FDR obtained from gene analysis is greater than the threshold, it is determined to be a high-risk gene.

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

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

[0034] In some of these solutions, step S5 specifically includes:

[0035] (1) Based on sample ID, gene and phenotype, form an ID-gene dataset with rare mutation load scores; traverse all genes and perform linear regression or logistic regression analysis on the mutation load scores of each gene and the related traits or diseases, where the total number of people carrying the gene mutation must be greater than 100.

[0036] (2) Risk candidate genes are determined 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-described 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 used to construct the gene-phenotype association analysis model described above;

[0040] Data analysis module: used to acquire whole genome sequencing or whole exome sequencing data of the sample to be evaluated, analyze and calculate according to the above model establishment method, obtain the rare mutation load score of the sample to be evaluated, calculate the correlation between the rare mutation load score and the trait or disease, and determine the risk genes 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 the use of a risk gene in functional assessment, wherein 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:

[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 gene: MC4R;

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

[0048] (6) Genes related to the estimation of glomerular filtration rate (EGCR) based on creatinine: ARG1, CSAG1, DKK1, IL10, MOSPD2 or PEG10;

[0049] (7) Genes related to the estimation of glomerular filtration rate (EGCY) based on cystatin: CLDN10, CST3, KCNN4, KRAS or TET2;

[0050] (8) Genes related to docosahexaenoic acid (DOA) metabolism: ANGPTL3, LDLR or LIPC;

[0051] (9) Genes related to glycated hemoglobin (Hemoglobin A1C, HBA1C): 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 gene: B3GNT7;

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

[0055] (13) Genes related to ω-3 fatty acid metabolism (Omega-3 Fatty Acids, OTFA): ANGPTL3 or LIPC;

[0056] (14) Genes related to ω-6 fatty acid metabolism (Omega-6 Fatty Acids, OSFA): ABCA1, ANGPTL3, APOA5, APOB, LCAT, LDLR, LIPC, LIPG, PCSK9 or TIMD4.

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

[0058] (16) Phosphoglyceride (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) Genes related to non-HDL cholesterol and non-LDL cholesterol (RMNC): ANGPTL3, APOB, LDLR or PCSK9;

[0062] (20) Sphingomyelin (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) Genes related to total fatty acid metabolism (TFA): ANGPTL3, ANGPTL4, APOA5, APOB, G6PC, LIPC, LIPG, PCSK9 or TIMD4;

[0066] (24) Genes related to total triglyceride metabolism (TTG): 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) Genes associated with bipolar disorder (BD): 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 gene: LDLR;

[0076] (9) Crohn's Disease (CD) related gene: 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 gene: MAP2K7;

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

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

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

[0085] (18) Genes associated with venous thromboembolism (VTE): REC8, PROC, PROS1 or JAK2.

[0086] Based on common knowledge in the field, the above-mentioned preferred conditions can be combined arbitrarily to obtain various preferred embodiments of the present invention.

[0087] The reagents and raw materials used in this invention are all commercially available.

[0088] The positive and progressive effects of this invention are as follows: It is a novel method for screening candidate genes for diseases or traits by calculating rare mutation burden scores. Through effective examples, it can be seen that most of the risk genes screened have been reported, demonstrating the reliability of this method. Furthermore, compared with traditional methods, some new associated genes are obtained in disease association results, most of which have already been reported. Therefore, this method can serve as an effective supplement to traditional methods in identifying disease-related pathogenic or risk genes. Using this model, mutation types and prediction software can be used to assess mutations, and candidate risk genes can be obtained through association analysis, leading to the discovery of new pathogenic or risk genes for diseases or traits. Compared with the traditional gene-based collapsing method, it has better reproducibility and complementarity. This model can play a supplementary role in finding associated candidate risk genes for new traits or unknown diseases. Attached Figure Description

[0089] Figure 1 This is the process for determining the optimal formula.

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

[0091] Figure 3 Calculate the correlation between gene-trait pairs for hyperparameters.

[0092] Figure 4 The optimal weight combination for the four major mutations.

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

[0094] Figure 6 This report presents information on disease-related genes and compares them with traditional collapsing methods. Detailed Implementation

[0095] The present invention is further illustrated below by way of embodiments, but the invention is not limited to the scope of the embodiments described herein. Experimental methods in the following embodiments that do not specify specific conditions were performed according to conventional methods and conditions, or as selected according to the product instructions.

[0096] Example 1: Constructing a gene-trait association prediction method

[0097] (1) Known traits and genes with related relationships

[0098] To calculate the association formula between a trait and the rare mutation score of the most suitable gene, it is first necessary to identify the known gene-trait pairs that are associated. Referring to the published article (Fiziev P, McRae J, Ulirsch JC, et al. Rarepenetrant mutations confer severe risk of common diseases. Preprint. medRxiv. 2023; 2023.05.01.23289356. Published 2023 May 8. doi:10.1101 / 2023.05.01.23289356), the rare mutations of the following genes and their associated traits are used as known gene-trait association combinations, as shown in Table 1.

[0099] Table 1. Rare mutations and phenotypic combinations of known associated genes

[0100]

[0101]

[0102] (2) Determination of the optimal formula

[0103] Rare gene mutations are classified into four categories: loss of function (LoF, which mainly includes frameshift, splicing, and stopgain mutations), non-synonymous mutations, synonymous mutations, and other protein-coding region mutations (including non-frameshift insertions, deletions, and start codon loss mutations), as shown in Table 2.

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

[0105]

[0106] Based on the type of rare mutation in an individual's gene, the number of alleles (AC), the allele frequency (AF), and the probability score for deleteriousness of non-synonymous mutations (Pr, 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 root mean square error (RMSE) through linear regression.

[0107] Each mutation type is scored individually using the following seven formulas. The scores for each mutation type are then summed to obtain the rare mutation score for that individual's gene. See the detailed procedure below. Figure 1 The data were validated by screening a population of 500,000 people in the UK biobank who had the gene mutations and phenotypes shown in Table 1. The screening criteria were: for each specific trait-gene pair, the subject must carry the mutation of that gene and have a test record for that trait. Table 3 shows the amount of data used in this embodiment.

[0108] Table 3. Total number of individuals and total number of mutations for each gene trait in this invention.

[0109]

[0110]

[0111] First, there is formula (1), the specific formula is as follows, only considering the influence of allele mutation number (AC). From Table 4, and Figure 2 The results of A show that AC can initially demonstrate the correlation between gene mutations and diseases, but the effect is generally limited in comparison.

[0112] S=AC (1)

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

[0114] Table 4. Regression intercept and correlation assessment 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. They are all evaluated based on the allelic mutation frequency, but different conversion methods are used: reciprocal, square root reciprocal, ln, log2 and log10. The specific evaluation results are shown in Table 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 assessment of formula (2)

[0124]

[0125]

[0126] Table 6. Regression intercept and correlation assessment of formula (3)

[0127]

[0128]

[0129] Table 7. Regression intercept and correlation assessment of formula (4)

[0130]

[0131]

[0132] Table 8. Regression intercept and correlation assessment of formula (5)

[0133]

[0134]

[0135] Table 9. Regression intercept and correlation assessment of formula (6)

[0136]

[0137]

[0138] Through R 2, Figure 2 The results showed that the reciprocal and the square root reciprocal were not as good as the log permutation. Moreover, the results were similar for the three different bases: 2, e and 10. Since the allele frequency showed stable performance and good model fitting effect when using decimal as the statistical standard, and was verified by the inventors' experiments, log10 was adopted as the permutation function of AF.

[0139] Finally, this invention further evaluates the impact of nonsynonymous mutation probability on the overall gene-trait association. Since LoF has been proven to directly associate rare gene mutations with traits in the form of AC, and synonymous mutations have almost been proven to have weak correlations, other mutations that may change protein type and structure are difficult to quantify individually due to their relatively low frequency in the genome (see Table 2). Therefore, this invention uses AC as its scoring benchmark. For nonsynonymous mutations, AC×Pr is used as its specific value in this embodiment. Regarding the assessment of the harmful probability of nonsynonymous mutations, the machine learning-based method proposed by Shen Maoting et al. (2024) is referenced (Shen Maoting, Lin Junwei, Fan Xijie, et al. Pathogenicity prediction of nonsynonymous mutations and its feature importance analysis based on machine learning [J]. Journal of Guangzhou Medical University, 2024, 52(5):1-9.DOI:10.3969 / j.issn.2095-9664.2024.05.01.). This method integrates more than 20 different prediction methods to probabilistically assess the pathogenicity of mutations, with a probability range of 0-1. Based on the probability, this embodiment defines the weight of non-synonymous mutations as 10 times the harmful probability, 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, fatmm-MKL, fatmm-XF, Eigen, GenoCanyon, fitCons, and AlphaMissense, etc.

[0140] In this embodiment, the random forest model is established by training with 1,344,221 non-synonymous mutations annotated as pathogenic, possibly pathogenic, benign, and possibly benign from the public database ClinVar20231112 (https: / / www.ncbi.nlm.nih.gov / clinvar / ).

[0141] The random forest model used can be built and trained using conventional methods in this field, such as constructing it using the "caret" package, where method = "rf", number = 10, repeats = 5, and the model is validated using a confusion matrix. The model's internal validation accuracy is 0.9265, and its AUC is 0.974.

[0142] Therefore, according to formula (7), the regression parameters are obtained using the same method. 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 assessment of formula (7)

[0145]

[0146]

[0147] In summary, this invention, by evaluating the impact of the AC, AF, and non-synonymous mutation harmfulness frequency scores (Pr) on gene-trait pairs, ultimately determined formula (7) as the optimal formula. Figure 2 A.

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

[0149] After determining the optimal formula for rare mutation scoring, linear regression was used to obtain the intercepts for different mutation types in each trait based on the association between scores for different mutation types and different traits. These intercepts served as initial weight references for each mutation type. Since the intercepts for different mutation types varied across traits, to find universal mutation type weights, the optimal intercepts (i.e., the highest R-values) for each gene-trait pair in different mutation types were calculated using the optimal formula. 2 The intercept of the values ​​was used as one of the input parameters; then, within these ranges, a grid search was used for evaluation with an interval of 0.1, resulting in a total of 14,658 combinations. Then, weighted combinations applicable to all traits were selected, using weighted scores and the correlation R-squared of the traits. 2 As an evaluation standard.

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

[0151]

[0152] Finally, among all suitable combinations, the optimal weights are selected based on the mean and median.

[0153] By calculating 14,658 combinations, the correlation R of the 16 gene-trait pairs from Table 1 for each combination was obtained. 2 Then calculate its mean and variance, and based on R... 2 Weight combinations with an average value greater than 0.1 were selected to identify the best-performing combinations (234 in total). See details below. Figure 3 See Table 11.

[0154] Table 11R 2 Weight combinations with an average value greater than 0.1 (exemplary examples show the top 50 weight combinations with good performance).

[0155]

[0156]

[0157] Finally, by using the mean and median, the optimal weighted sums for the four mutation types were selected: LoF = 0.8, nonsynonymous mutations = 0.70, other mutations = 0.30, and synonymous mutations = 0.00. (See [link to LoF analysis]). Figure 4 See 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 formula for rare mutation load scores below is used to calculate the rare mutation load score for each person and each gene. Furthermore, through association analysis of phenotype and rare mutation load scores—namely, linear regression of traits and logistic regression of diseases—risk candidate genes for traits or diseases are screened. (See below.) Figure 2 B.

[0159]

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

[0161] (4) Screening for pathogenic or risk genes by associating harmful mutational load fraction with disease / trait through regression methods.

[0162] The rare mutation load score for each sample and each gene is calculated using the formula for rare mutation load score. Then, association analysis is performed based on specific traits (such as low-density lipoprotein) or diseases (such as type 2 diabetes) using linear regression (for traits of continuous variables) or logistic regression (for diagnosing the presence of diseases).

[0163] The specific steps are as follows: First, based on the sample ID, gene, and trait, a dataset of ID and gene is constructed using rare mutation load scores. Then, by reviewing all genes, linear regression or logistic regression analysis is performed on the mutation load scores of each gene and the related traits or diseases, requiring that the total number of people carrying the gene mutation be greater than 100. The p-values ​​obtained from the regression analysis are corrected for FDR, for example, using the Benjamini-Hochberg method. Finally, risk candidate genes are determined based on the corrected FDR values, and risk candidate genes related to the relevant traits or diseases are screened by setting FDR thresholds, where the FDR for continuous traits is <0.001 and the FDR for diseases is <0.05.

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

[0165] This embodiment provides a specific implementation scheme for 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 known associated trait and gene data as rare gene mutation types and trait association combinations, and preprocess these data, including quality control, outlier handling, and standardization, to prepare for subsequent analysis.

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

[0169] The collected mutations were labeled 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 other than the above mutation types (including non-frameshift insertion / deletion mutations, start codon loss mutations, and stop codon loss mutations).

[0170] The following formula for scoring rare gene mutations is used to calculate scores for various mutation types:

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

[0172] Where: S is the score for rare gene mutation types, AC is the number of alleles, AF is the frequency of allele mutations, and Pr is the probability score of harmfulness, 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. This model takes 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 relatively small impact, the Pr value is set to 1.

[0175] (3) Determining weight combinations through hyperparameter optimization

[0176] The weight combinations for different mutation types were evaluated using a hyperparameter method (manual parameter tuning). The optimization process is as follows:

[0177] The optimal intercepts of each gene-trait pair calculated using the rare gene mutation scoring formula for different mutation types were used as input parameters. Within these parameter ranges, parameter combinations were set with an interval of 0.1 to obtain multiple weight combination schemes. The correlation R of the known gene-trait pairs determined in step (1) corresponding to each weight combination was calculated. 2 Calculate its mean and variance, etc., according to R. 2 The weight combination applicable to all traits was selected based on the mean and median. The optimal weight combination was determined as follows: LoF = 0.8, non-synonymous mutation = 0.7, other mutations = 0.3, and synonymous mutation = 0.

[0178] (4) Calculate the rare mutation load fraction of genes.

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

[0180]

[0181] Wherein: GBS is the gene mutation load 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 fraction and traits / diseases, an association analysis model was established: for specific traits (such as low-density lipoprotein), linear regression was used with GBS as the independent variable and trait value as the dependent variable; for disease states (such as type 2 diabetes), logistic regression was used with GBS as the independent variable and disease state (0 / 1) as the dependent variable.

[0184] Based on sample ID, gene, and trait, a dataset of ID-genes is constructed using rare mutation load scores. All genes are traversed, and the mutation load score of each gene is correlated with the relevant trait or disease, wherein the total number of people carrying the gene mutation must be greater than 100.

[0185] The p-values ​​obtained from the regression analysis were corrected for FDR, and risk candidate genes were determined 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 steps (1)-(5) above, the system successfully established a gene-trait association analysis model based on rare mutation scores.

[0187] Example 1 validates the gene-trait association analysis model based on rare mutation scores.

[0188] (1) Screening for trait-associated risk genes using harmful mutation load fractions

[0189] To verify that the harmful mutation load score in the gene-trait association analysis model based on rare mutation scores in Example 2 can effectively screen for trait-related risk genes, data from a UK biobank population of 500,000 people were used to screen for gene mutations listed in Table 12 and 25 different trait phenotypes. The screening criteria were: for each specific trait-gene pair, the subject must carry the mutation of that gene and have a test record for that trait. Linear regression analysis was performed using the model in Example 2. The mutation load score, constructed using the most suitable mutation type weights and mutation factors (AF, log10, 1 / AF, and Pr), was used to obtain the reliability of the associated genes by association with the trait. Using FDR < 0.001, 149 genes associated with 24 traits were screened, of which 128 gene-trait pairs 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 for trait-related genes.

[0191]

[0192]

[0193]

[0194]

[0195]

[0196]

[0197]

[0198]

[0199] (2) Screening for disease-related genes by harmful mutation load fraction

[0200] To verify that the harmful mutation load score in the gene-trait association analysis model based on rare mutation scores in Example 2 can effectively screen for disease-related risk genes, data from a UK biobank population of 500,000 people were used to screen for gene mutations listed in Table 13 and 28 different disease phenotypes. The screening criteria were: for each specific trait-gene pair, the subject must carry the mutation of that gene and have a test record for that disease. Logistic regression analysis was performed using the model in Example 2, and 55 genes associated with 18 diseases were screened using FDR < 0.05. Among them, 45 disease-gene pairs were reported, accounting for 81.82%. Figure 6 (See Table 13). This result also demonstrates the reliability of association analysis based on gene mutation burden scores. Furthermore, to further validate its efficacy, 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. Published 2023 May 8. doi:10.1101 / 2023.05.01.23289356) was used to validate and compare the same group data screened from a population of 500,000 in the UK biobank. Collapsing screened 109 associated genes for 22 diseases, of which 83 disease-gene pairs had been reported, accounting for 76.15%. Figure 6 (A, Table 14) The reported rates are similar to those of the methods used, but the total number is greater. It was found that 33 disease-gene pairs overlapped between the two methods, all of which were reported disease-gene pairs. Figure 6 (B). Compared to the traditional collapsing method, this method should be an effective supplement to it, and can be helpful in finding more new disease-associated genes.

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

[0202]

[0203]

[0204]

[0205] Table 14 Collapsing Analysis Results

[0206]

[0207]

[0208]

[0209] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these are merely illustrative examples, and various changes or modifications can 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, Includes the following steps: S1: Collect known related trait and gene data to form gene-trait pairs, which serve as combinations of rare gene mutation types and trait associations; S2: Calculate the score for each rare mutation type of gene using the following rare mutation scoring formula: Where S is the score of rare gene mutation type, 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, use linear regression to analyze the correlation between the scores of rare mutation types and traits in gene-trait pairs in S1, calculate the intercept of each rare mutation type in each trait as a reference value for the weight of each mutation type, and evaluate the weight combination of each rare mutation type in different genes through hyperparameter optimization method. Use the correlation R² between the score and the trait as the evaluation standard. When the R² value exceeds the preset threshold, it is determined as the weight combination of rare mutation types, which is used as the parameter for subsequent calculation of rare mutation load score; S4: Apply the scoring formula determined in step S2 and the weight combination determined in step S3, and calculate the rare mutation load score for each gene in each sample according to the following rare mutation load score formula: Where GBS is the gene mutation load score; i represents the rare gene mutation type; Beta is the weight of rare gene mutation type i; and S is the score of rare mutation type i calculated according to the rare gene mutation scoring formula. S5: Analyze the association between rare mutation load score and phenotype using linear regression or logistic regression methods to obtain a gene-phenotype association analysis model based on rare mutation score.

2. The method for establishing as described in claim 1, characterized in that, In step S1, the rare gene mutation types include at least one of the following: loss-of-function mutations, non-synonymous mutations, synonymous mutations, and other protein-coding region mutations other than the above mutation types.

3. The method for establishing as described in claim 2, characterized in that, The loss-of-function mutations include at least one of frameshift mutations, splicing mutations, and truncation mutations; and / or, other protein-coding region mutations not of the above mutation types include at least one of non-frameshift insertion / deletion mutations, start codon loss mutations, and stop codon loss mutations.

4. The method for establishing as described in claim 2 or 3, characterized in that, In step S2, the Pr value is used to evaluate nonsynonymous mutations and is obtained by: using at least three of the following data analysis software programs to evaluate the pathogenicity of each rare mutation: SIFT, PolyPhen-2 (HDIVandHVAR), LRT, MutationTaster, MutationAssessor, FATHMM, PROVEAN, VEST4, MVP, MPC, PrimateAI3D, DEOGEN2, LIST.S2, CADD, DANN, fatmm-MKL, fatmm-XF, Eigen, GenoCanyon, fitCons, and AlphaMissense; using the evaluation values ​​of each data analysis system as input; and using a random forest model to predict the pathogenicity probability of the rare mutation, which is the Pr value of the nonsynonymous mutation; and / or, Pr = 1 for other mutation types.

5. The method for establishing as described in claim 4, characterized in that, The random forest model was established by training on non-synonymous mutations annotated as pathogenic, potentially pathogenic, benign, and potentially benign in the public database ClinVar20231112.

6. The method for establishing as described in claim 1, characterized in that, The hyperparameter evaluation method in step S3 includes the following steps: (1) Collect the optimal intercepts of different mutation types in each gene-trait pair calculated by the rare mutation scoring formula as parameters for input linear regression analysis; the optimal intercept refers to the correlation between the mutation type and the R-value. 2 The intercept at its highest point; (2) Within the range of the input parameters, multiple weight combination schemes are obtained by setting the parameter combination with an interval of 0.1±0.

02. (3) Calculate the correlation R² of the known gene-trait pairs determined in step S1 for each weight combination, and calculate their average value. 2 The weight combination with the highest average value is used as the weight combination applicable to all phenotypes; the R 2 The average value is greater than 0.

1.

7. The method for establishing as described in claim 6, characterized in that, In step S4, the weight combination is as follows: loss of function mutation is 0.10-1.00, non-synonymous mutation is 0.10-1.00, other mutation is 0.00-0.80, and synonymous mutation is 0.00-0.

10.

8. The method for establishing as described in claim 7, characterized in that, The weighting combination is as follows: loss-of-function mutations are weighted at 0.70-0.80, non-synonymous mutations at 0.40-1.00, other mutations at 0.00-0.60, and synonymous mutations at 0.00-0.

1.

9. The method for establishing as described in claim 8, characterized in that, The weight combination is as follows: missing-function mutation = 0.8, non-synonymous mutation = 0.7, other mutations = 0.3, and synonymous mutation = 0.

10. The method for establishing as described in claim 1, characterized in that, The phenotype includes traits and / or diseases. When the phenotype is a trait, in step S5, the association between the rare mutation load fraction 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 logistic regression.

11. The method for establishing as described in claim 10, characterized in that, Step S5 further includes: obtaining a P-value through linear regression or logistic regression, correcting the P-value with FDR to obtain an FDR value, and setting an FDR threshold to determine the correlation between the analyzed gene and the target trait or disease. If the FDR obtained from gene analysis is greater than the threshold, it is judged as a high-risk gene. FDR correction is performed using the Benjamini-Hochberg method.

12. The method for establishing as described in claim 11, characterized in that, The FDR threshold for the trait is less than 0.001, and the FDR threshold for the disease is less than 0.

05.

13. The method for establishing as described in claim 11 or 12, characterized in that, The specific steps of step S5 are as follows: (1) Based on the sample ID, gene and phenotype, form an ID-gene dataset with rare mutation load scores; traverse all genes and perform linear regression or logistic regression analysis on the mutation load scores of each gene and the related traits or diseases, where the total number of people carrying the gene mutation must be greater than 100. (2) Risk candidate genes are determined based on the corrected FDR value.

14. A gene-phenotype association analysis model, characterized in that, Obtained by the method described in any one of claims 1-13.

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

16. The use of a risk gene identified by the establishment method according to any one of claims 1-13 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) Body mass index-related gene: MC4R; (5) Calcium metabolism-related genes: ALB, ALPL, CASR, CNMN4 or FCGRT; (6) Genes related to estimating glomerular filtration rate based on creatinine: ARG1, CSAG1, DKK1, IL10, MOSPD2 or PEG10; (7) Genes related to estimating glomerular filtration rate based on cystatin: CLDN10, CST3, KCNN4, KRAS or TET2; (8) Genes related to docosahexaenoic acid metabolism: 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 gene: B3GNT7; (12) Low-density lipoprotein cholesterol-related genes: ANGPTL3, APOB, ASGR1, LDLR or PCSK9; (13) Genes related to ω-3 fatty acid metabolism: ANGPTL3 or LIPC; (14) Genes related to ω-6 fatty acid metabolism: 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 glycerol phosphate metabolism: ABCA1, ANGPTL3, APOA1, ICOSLG, LCAT, LIPC, LIPG or SCARB1; (17) Genes related to polyunsaturated fatty acid metabolism: ABCA1, ANGPTL3, APOA5, APOB, LCAT, LDLR, LIPC, LIPG, PCSK9 or TIMD4; (18) Resting heart rate related genes: KCNJ5, PTEN, RGS6 or RRAD; (19) Residual cholesterol-related genes for non-high-density lipoprotein cholesterol and non-low-density lipoprotein cholesterol: ANGPTL3, APOB, LDLR or PCSK9; (20) Genes related to sphingomyelin metabolism: ABCA1, ANGPTL3, LCAT, LDLR, LIPC, LPL, PCSK9 or SCARB1; (21) Height-related genes: GHRH, HMGA2, IHH, NPR2, NRK, PTPN11, SOCS2 or SPIN4; (22) Genes related to total cholesterol metabolism: ABCA1, ANGPTL3, APOB, LCAT, LDLR or PCSK9; (23) Genes related to total fatty acid metabolism: ANGPTL3, ANGPTL4, APOA5, APOB, G6PC, LIPC, LIPG, PCSK9 or TIMD4; (24) Genes related to total triglyceride metabolism: 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) Genes associated with bipolar disorder: 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 gene: LDLR; (9) Crohn's disease-related gene: 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 gene: MAP2K7; (15) Primary open-angle glaucoma-related gene: LTBP2; (16) Prostate cancer-related gene: ATM; (17) Type 2 diabetes-related genes: GCK, GIGYF1, C19orf57, HNF1A, MIB1, HMGXB4, IGF1R or HNF4A; (18) Genes associated with venous thromboembolic diseases: REC8, PROC, PROS1 or JAK2.

Images