A method and apparatus for predicting individual phenotypes based on human whole-genome genotypes
By combining deep neural networks and gradient ascent networks, the problem of low utilization efficiency of whole genome data in existing technologies is solved, and more accurate individual phenotype and disease risk prediction is achieved, which is applicable to individual phenotype prediction of various populations.
Patent Information
- Application Number
- CN202210837751.3
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2022-07-15
- Publication Date
- 2025-12-16
- Estimated Expiration
- 2042-07-15
AI Technical Summary
Existing algorithms cannot effectively utilize whole-genome genotyping data, making it difficult to accurately predict individual phenotypes or disease risks. Furthermore, the models do not perform well across different populations because linear models cannot reflect the complex nonlinear interactions and linkage disequilibrium between genes.
A deep neural network is used to convert haplotypes into functional genomic parameters. By segmenting regions and selecting representative parameters, combined with regression models and gradient ascent networks, the association between each region and individual phenotype is integrated, and prediction is performed using variable screening models and gradient ascent networks.
It improves the accuracy of individual phenotype prediction, reflects the non-linear relationship between genes at the single gene level, overcomes the limitations of linear association, and is applicable to phenotype prediction for various populations.
Smart Images
Figure CN115547408B_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The present application belongs to the technical field of phenotype prediction, and particularly relates to a method and device for predicting individual phenotypes based on human whole genome genotypes. BACKGROUND
[0002] Many human traits and diseases are controlled by a large number of genetic loci. In theory, the performance of a certain genetic trait or the incidence risk of a certain genetic disease of an individual can be accurately estimated by using the genotype data of the whole genome of the individual. At present, the rapid development of high-throughput sequencing and whole genome chips makes the cost of obtaining the whole genome data of an individual lower and lower, and the development of genome data processing software and algorithms also makes the data analysis speed faster and faster. However, the existing algorithms still cannot well complete the task of predicting phenotypes or disease risks according to genotypes, which is specifically as follows:
[0003] The existing algorithms directly use genotypes as input and traits or diseases as output, and establish the association between genetic variations and traits or diseases, but these association results still lack biological mechanism explanation and are difficult to predict disease risks. The reason for this phenomenon is that the existing algorithms all use linear models for association analysis, while there are complex nonlinear interactions between genes. In addition, there is linkage disequilibrium between genetic loci, and there are significant genetic differences between populations, so that the results or prediction models obtained in this population are difficult to be applied in another population. Therefore, although a large amount of whole genome genotyping data has been accumulated, there is still a lack of methods that can effectively use these data. SUMMARY
[0004] The purpose of the present application is to provide a method and device for predicting individual phenotypes based on human whole genome genotypes, to solve the technical problem of the prior art that there is a lack of methods that can effectively use whole genome genotyping data.
[0005] In order to achieve the above purpose, the present application adopts the following technical solutions:
[0006] In a first aspect, the present application provides a method for predicting individual phenotypes based on human whole genome genotypes, comprising:
[0007] Obtaining haplotypes in all regions of the whole genome of each individual, converting the haplotypes into functional genome parameters based on a deep neural network, and selecting at most one representative parameter from the functional genome parameters of each region;
[0008] Quantitatively analyzing the association between each region and individual phenotypes by using a regression model according to the representative parameters of the regions;
[0009] According to the risk score of each gene and the representative parameter of the regulatory region where the gene is located, the phenotype prediction value of each gene is obtained by combining the variable screening model and the first gradient ascent network prediction;
[0010] The phenotype prediction value of each gene is input into the second gradient ascent network for integration to generate a final prediction result.
[0011] Based on the above disclosure, the present application does not directly use genotype data, but converts genotype (haplotype) into functional genomic parameters through a deep neural network. Since the complex interaction relationship between genes is considered, the representative parameters are selected from the functional genomic parameters, and the association between each region and individual phenotype and the prediction value of single gene are analyzed by using the representative parameters of the region, thereby integrating the information from each functional genomic change at the single gene level, avoiding the limitation of linear association of isolated analysis of each gene site. Finally, the prediction value of all genes is integrated through the gradient ascent network to fully reflect the non-linear relationship between genes, thereby improving the prediction accuracy and having important significance for the prediction of various individual phenotypes.
[0012] In one possible design, the haplotype in all regions of the whole genome of each individual is obtained, the haplotype is converted into functional genomic parameters based on a deep neural network, and at most one representative parameter is selected from the functional genomic parameters of each region, including:
[0013] The whole genome of each individual is regionally segmented to obtain the haplotype and the corresponding base sequence in each region;
[0014] Each base sequence is converted into a plurality of corresponding functional genomic parameters based on a deep neural network;
[0015] At most one representative parameter is selected from the plurality of functional genomic parameters of each region based on a fine mapping algorithm.
[0016] Based on the above disclosure, the present application regionally segments the whole genome to make the base sequence obtained after segmentation meet the data input length requirement of the neural network. By converting genotype data into haplotype data, the calculation resources can be saved. By converting each base sequence into a plurality of functional genomic parameters through a deep neural network, the functional genomic changes caused by gene variation can be quantitatively reflected. By selecting a representative parameter, the prediction accuracy of the model can be improved.
[0017] In one possible design, the whole genome of each individual is regionally segmented to obtain the haplotype and the corresponding base sequence in each region, including:
[0018] The whole genome of each individual is divided into multiple regions according to a preset length, and each region is ensured to cover functional elements of the whole genome as completely as possible.
[0019] The genotypes are screened according to non-allele frequencies, the screened genotypes are converted into haplotypes, and the base sequences corresponding to each haplotype are obtained.
[0020] In a possible design, before selecting at most one representative parameter from the multiple functional genomic parameters of each region based on the fine mapping algorithm, the method comprises:
[0021] The multiple functional genomic parameters are integrated to obtain multiple category scores of each base sequence, and the multiple category scores are processed.
[0022] Each processed category score is converted according to an additive genetic model, a recessive genetic model and a dominant genetic model to generate a corresponding intermediate parameter, so that the fine mapping algorithm selects at most one representative parameter from the intermediate parameters of each region.
[0023] Based on the above disclosure, by integrating the functional genomic parameters to obtain multiple category scores, the amount of data information can be improved while reducing the amount of data.
[0024] In a possible design, the association between each region and the individual phenotype is quantitatively analyzed using a regression model according to the representative parameter of the region, including:
[0025] In each region of the representative parameter, the association between the region and the individual phenotype is quantitatively analyzed based on a linear regression model.
[0026] The effective region is selected from the region with linkage disequilibrium in the interior using a fine mapping algorithm.
[0027] The biological significance annotation of each region is obtained, and whether the effective region affecting the individual phenotype is more related to a certain biological significance annotation is quantitatively analyzed based on a generalized linear regression model.
[0028] Based on the above disclosure, the association between each region and the individual phenotype and whether the effective region affecting the individual phenotype is more related to a certain biological significance annotation can be determined, so that the region with a large association or the biological significance annotation can be focused on, and the data analysis efficiency can be improved.
[0029] In a possible design, the phenotype prediction value of each gene is predicted by combining a variable screening model and a first gradient ascent network according to the risk score of each gene and the representative parameter of the regulatory region where the gene is located, including:
[0030] The preprocessed gene is input into a variable screening model for calculation, and according to the output gene site weight, the exon polygenic risk score of the gene is calculated.
[0031] The exon polygenic risk score and all representative parameters in the regulatory region of the gene are input into the first gradient ascent network to obtain the phenotype prediction value of the gene.
[0032] In a possible design, the preprocessed gene is input into a variable screening model for calculation, including:
[0033] The preprocessed gene is input into a variable screening model for calculation, and according to the output gene site weight, the exon polygenic risk score of the gene is calculated.
[0034] In the regression analysis, a covariate is added, and the parameters of the covariate are not scaled by lasso, so as to output the gene site weight without linkage disequilibrium.
[0035] In a possible design, the exon polygenic risk score and all representative parameters in the regulatory region of the gene are input into the first gradient ascent network to obtain the phenotype prediction value of the gene, including:
[0036] According to the functional genomic annotation, the regulatory region of each gene is preset;
[0037] The exon polygenic risk score and all representative parameters in the regulatory region of the gene are input into the first gradient ascent network to obtain the phenotype prediction value of the gene.
[0038] In a possible design, the phenotype prediction value of each gene is input into the second gradient ascent network for integration to generate a final prediction result, including:
[0039] The phenotype prediction value of each gene is input into the second gradient ascent network for integration to generate a final prediction result, including:
[0040] In a possible design, the phenotype prediction value of each gene is input into the second gradient ascent network for integration to generate a final prediction result, including:
[0041] The parameter conversion module is configured to obtain haplotypes in all regions of the whole genome of each individual, convert the haplotypes into functional genomic parameters based on a deep neural network, and select at most one representative parameter from the functional genomic parameters of each region.
[0042] The association analysis module is configured to quantitatively analyze the association between each region and the individual phenotype according to the representative parameter of the region by using a regression model.
[0043] The first prediction module is configured to obtain a phenotype prediction value of each gene by combining a variable screening model and a first gradient ascent network according to the risk score of each gene and the representative parameter of the regulatory region where the gene is located.
[0044] The second prediction module is configured to input the phenotype prediction value of each gene into a second gradient ascent network for integration to generate a final prediction result.
[0045] In a third aspect, the present application provides a computer device, comprising a memory, a processor and a transceiver connected in sequence and in communication, wherein the memory is configured to store a computer program, the transceiver is configured to receive and send messages, and the processor is configured to read the computer program and execute the method in any possible design of the first aspect.
[0046] In a fourth aspect, the present application provides a computer readable storage medium, wherein the computer readable storage medium stores instructions, and when the instructions are executed on a computer, the method in any possible design of the first aspect is executed.
[0047] In a fifth aspect, the present application provides a computer program product comprising instructions, and when the instructions are executed on a computer, the computer is caused to execute the method in any possible design of the first aspect. BRIEF DESCRIPTION OF DRAWINGS
[0048] Figure 1 A flowchart of the method for predicting individual phenotypes based on human whole-genome genotypes provided by the present application is shown in the figure.
[0049] Figure 2 A flowchart of step S1 provided by the present application is shown in the figure.
[0050] Figure 3 A flowchart of step S2 provided by the present application is shown in the figure.
[0051] Figure 4 A flowchart of step S3 provided by the present application is shown in the figure. DETAILED DESCRIPTION
[0052] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the present application will be briefly introduced below in combination with the drawings and the description of the embodiments or the prior art. Obviously, the following description of the drawings is only some embodiments of the present application, and other drawings can also be obtained by those skilled in the art without creative labor. It should be noted that the description of these embodiments is used to help understand the present application, but does not constitute a limitation on the present application.
[0053] Embodiments
[0054] In order to overcome the technical problems of the prior art that there is no method that can effectively utilize whole genome genotyping data, the present embodiment provides a method for predicting individual phenotypes based on human whole genome genotype, which integrates information from each functional genome change at the single gene level, avoids the limitation of linear correlation of isolated analysis of each gene site, and finally integrates all gene prediction values through gradient ascent network to fully reflect the non-linear relationship between genes, thereby improving the prediction accuracy and having important significance for predicting various individual phenotypes.
[0055] It should be noted that the method provided by the present application can be applied to various prediction scenarios of human phenotypes, including but not limited to height, weight, etc., and the sampling population of the genotype data used can also be diverse. In order to facilitate the description, the present embodiment will be described in detail by taking the use of European population whole genome chip typing data to construct a model for predicting individual height level as an example. Of course, it can be understood that based on the description below, those skilled in the art have the motivation to apply the present embodiment to data of other populations to construct a model for predicting other traits, so other embodiments are also within the protection scope of the present application. For the convenience of understanding, some software tools, algorithm models and algorithm parameters that can be used are specifically described in the present embodiment. Of course, it can be understood that replacing these software tools, algorithm models and algorithm parameters with other similar tools does not represent a creative change to the present embodiment, and is within the protection scope of the present application.
[0056] The method for predicting individual phenotypes based on human whole genome genotype provided by the present application will be described in detail below.
[0057] As shown in Figures 1-4 the flowchart of the method provided by the present application, the method includes but is not limited to being implemented by steps S1-S4:
[0058] Step S1. Obtain haplotypes in all regions of the whole genome of each individual, convert the haplotypes into functional genome parameters based on a deep neural network, and select at most one representative parameter from the functional genome parameters of each region;
[0059] In the embodiment, preferably, the deep neural network adopts a sei network model.
[0060] In step S1, haplotypes in all regions of the whole genome of each individual are obtained, the haplotypes are converted into functional genomic parameters based on a deep neural network, and at most one representative parameter is selected from the functional genomic parameters of each region, including:
[0061] Step S11. The whole genome of each individual is segmented into regions, and haplotypes and corresponding base sequences in each region are obtained, specifically including:
[0062] Step S111. The whole genome of each individual is segmented into regions according to a preset length, and each region is ensured to cover functional elements of the whole genome as completely as possible.
[0063] It should be noted that when the sei network model is adopted, since the input length specified by the sei model is 4096 bp, the whole genome is segmented into regions of 4096 bp in length in the embodiment, and the regions are ensured to cover all functional elements of the genome as completely as possible. Specifically, first, based on the chromatin state annotation of 222 human tissues in epimap, for each tissue, all chromosomal regions annotated as “transcription start site (TSS) transcription region (TX), enhancer (enhancer), promoter (promoter)” in at least half of the samples are marked as active regions, and the active regions of all tissues are taken and collected, and the regions annotated as genomic gap in Hg38 genome are removed. Then, for the series of active regions, if the length is less than 4096 bp, the active region is defined as a 4096 bp region centered on the active region, and if the length is greater than 4096 bp, a 4096 bp region is gradually defined from the midpoint of the region. Finally, the remaining genomic regions are covered by non-repeating 4096 bp blocks.
[0064] Step S112. The genotypes are screened according to the allele frequency, the screened genotypes are converted into haplotypes, and the base sequences corresponding to each haplotype are obtained.
[0065] It should be noted that the embodiment assumes that the.bgen format genotype data of all individuals has been obtained in advance, and haplotyping has been performed. Specifically, the genotypes are screened according to the allele frequency, including: for each region, first use the bgenix tool to obtain all variations located in the region from the.bgen file, remove the variations with quality control parameters INFO<0.8, hardy-wenberg test p value less than 10 -6, the total number of variants is less than 10, and the individuals with more than 50% of the variants lost in the region are removed. It should be understood that the above values are based on the sample population and the type of trait to be predicted in this embodiment. If other samples or other traits are used, the above values can be adjusted according to the needs, which are not limited here.
[0066] Preferably, after screening the genotype, the screened.bgen data can be converted into.haps data using plink, and the.vcf data of all haplotypes and the haplotype genotype of all individuals can be obtained by awk and datamash. The base sequence corresponding to all haplotypes can be obtained by awk and bcftools, and stored in a fasta file.
[0067] Step S12. Converting each base sequence into a plurality of functional genomic parameters based on a deep neural network;
[0068] In step S12, the fasta file storage is input into the sei algorithm network, and the data reading step uses the fasta format processing function provided by the selene toolkit. The sei algorithm network includes three steps when converting data: 1) a convolutional network with linear and nonlinear dual paths to encode sequence information, 2) a residual dilation convolution layer, and 3) a spatial basis function transformation and an output layer. Through these three steps, the sei network can calculate 21906 functional genomic parameters corresponding to each base sequence.
[0069] Preferably, after obtaining the functional genomic parameters, the method comprises:
[0070] 1) integrating the plurality of functional genomic parameters to obtain a plurality of category scores of each base sequence, and processing the plurality of category scores;
[0071] Specifically, the sei calculates 21906 functional genomic parameters corresponding to each base sequence to obtain a functional genomic parameter matrix . Then, based on the matrix, 40 sequence category score matrices are calculated according to the following formula:
[0072] ; (1)
[0073] wherein, is a projection matrix obtained based on UMAP dimension reduction and clustering analysis, which calculates the possibility of each base sequence belonging to a certain sequence category (40 sequence categories in total) according to the functional genomic parameters of the base sequence, and is represented by . The final output matrix For n rows and 40 columns, that is, each base sequence corresponds to 40 sequence category scores. The greater the score, the greater the likelihood that the base sequence belongs to this sequence category.
[0074] Preferably, after the sei prediction of all base sequences in a region and the hg38 reference sequence of the sequence are completed, the method further comprises completing the standardization of the scores using the following formula:
[0075] ; (2)
[0076] Wherein, score histone represents the histone mark parameter, ref represents the reference sequence of each region, represents the sum of all histone mark parameters corresponding to the reference sequence in each region.
[0077] Preferably, for the standardized category scores, the method further comprises:
[0078] Eliminate scores in the sequence category score matrix that do not meet the preset conditions, such as the category scores of low signal, heterochromatin, and polycomb, and calculate the coefficient of variation of each remaining sequence category score , the calculation formula is as follows:
[0079] ; (3)
[0080] Wherein, represents the standard deviation of each sequence category score, represents the average of the standard deviation of each sequence category score.
[0081] 2) Transform each category score after processing according to the additive genetic model, the recessive genetic model, and the dominant genetic model to generate corresponding intermediate parameters, so that the fine mapping algorithm selects at most one representative parameter from the intermediate parameters of each region.
[0082] Specifically, step 2) comprises:
[0083] Select the first 10 sequence category scores with the largest coefficient of variation , respectively transform them according to the additive genetic model, the recessive genetic model, and the dominant genetic model to obtain 30 intermediate parameters. Of course, it can be understood that the above values can be adjusted according to specific needs, which are not limited here.
[0084] Specifically, transforming according to the additive genetic model includes taking the average of the category scores of the two haplotypes; transforming according to the recessive genetic model includes taking the minimum of the category scores of the two haplotypes; and transforming according to the dominant genetic model includes taking the maximum of the category scores of the two haplotypes.
[0085] Step S13. Selecting at most one representative parameter from the plurality of functional genomic parameters of each region based on a fine mapping algorithm.
[0086] Preferably, the fine mapping algorithm in the embodiment adopts the SUSIE algorithm, which selects at most one representative parameter from the plurality of functional genomic parameters of each region, including:
[0087] Variable screening is performed on the 30 intermediate parameters using the SUSIE algorithm, and the parameter settings are all selected as default values. Among them, the result variable is selected as the residual of the following regression model:
[0088] ; (4)
[0089] Among them, , , , , , respectively represent height, age, gender, sample source, genetic principal component and other user-specified covariates. SUSIE divides the independent variables in a regression model containing a large number of linearly correlated independent variables into several sets and estimates the posterior probability that each set contains at least one truly effective independent variable. Within each set, SUSIE further estimates the posterior probability that each independent variable is a truly effective independent variable. The mathematical expression of this process is as follows:
[0090] ; (5)
[0091] Among them, represents the coefficient of each independent variable estimated from the data, represents the coefficient of the truly effective independent variable, represents a 0 / 1 variable, and the task of SUSIE is to estimate the posterior probability that each item of is 1.
[0092] From the SUSIE results, a group of intermediate parameters with the largest probability is selected, and one intermediate parameter with the largest posterior probability is selected. If the inclusion probability of all groups is less than 0.5, skip this region and do not perform analysis. The selected one intermediate parameter is used as the representative parameter of the region for subsequent analysis.
[0093] Step S2. Quantitative analysis of the association between each region and individual phenotype according to the representative parameters of the region using a regression model;
[0094] In step S2, the association between each region and individual phenotype is quantitatively analyzed according to the representative parameters of the region using a regression model, including:
[0095] Step S21. In each region of the representative parameters, the association between the region and individual phenotype is quantitatively analyzed based on a linear regression model, and the linear regression model is expressed as follows:
[0096] ; (6)
[0097] wherein, , , , , , The regression coefficient and p value of the representative parameter in the regression model are used as the association result between the region and the phenotype.
[0098] Step S22. Use a fine mapping algorithm to select the onset region from the region with linkage disequilibrium inside;
[0099] Specifically, the European hg38 version of the human genome data is used for LDetect analysis to obtain the whole genome linkage disequilibrium block. The variations inside these blocks have linkage disequilibrium, and the blocks are basically free of linkage disequilibrium. In each block, all the representative parameters of the region are input into SUSIE, and the analysis method is the same as that in the first step of SUSIE. The posterior probability > 0.8 output by SUSIE is used as the threshold to determine the onset region.
[0100] Step S23. Obtain the biological significance annotation of each region, and quantitatively analyze whether the onset region affecting individual phenotype is more related to a certain biological significance annotation based on a generalized linear regression model.
[0101] Specifically, first, three types of biological significance annotations are collected, and all regions of each individual are annotated as follows:
[0102] 1) Tissue-specific activity: Obtain chromatin state annotations of 222 human tissues from epimap, and obtain the active positions of the genome in each tissue according to the foregoing method.
[0103] 2) Cell-specific activity: Obtain human single-cell chromatin accessibility sequencing peak data, and obtain the active positions of the genome in each cell.
[0104] 3) Biological function: download all biological pathway annotation data from msigdb database, exclude pathways containing less than 10 genes or more than 500 genes, construct gene x pathway annotation matrix (each element is 0 or 1), perform hierarchical clustering based on this matrix and sort all paths. From one end of the sorting, gradually include the path into the candidate set to the other end. Whenever the adjacent two paths have more than 80% overlap, delete the smaller path and include the larger path into the candidate set until all paths are deleted or included. For each region, if it exists within the 100000bp range of the genes in each path, it is annotated as affecting the path.
[0105] Based on the above disclosure, for each annotation, all regions have two states of 1 and 0. For each annotation, construct the following generalized regression analysis model to analyze the following problem: whether the regions affecting the phenotype are more related to this annotation; wherein the expression of the generalized regression analysis model is as follows:
[0106] ; (7)
[0107] wherein, indicates whether the region meets the annotation (1 / 0), indicates the posterior probability that the region is an effective region calculated by SUSIE, indicates the total number of haplotypes in the region, indicates the total number of samples included in the analysis of the region, indicates the sequence class score corresponding to the reference sequence of the region, indicates the chromosome recombination frequency in the region, indicates the background selection intensity in human history of the region, indicates the conservation of the region, indicates the block number as a random effect term. Among them, the regression coefficient and p value of the term are used to measure the correlation between the trait heritability and the annotation.
[0108] Step S3. According to the risk score of each gene and the representative parameter of the regulatory region where the gene is located, the phenotype prediction value of each gene is predicted by combining the variable screening model and the first gradient ascent network.
[0109] In step S3, according to the risk score of each gene and the representative parameter of the regulatory region where the gene is located, the phenotype prediction value of each gene is predicted by combining the variable screening model and the first gradient ascent network, including:
[0110] Preferably, the variable screening model uses the snpnet algorithm model, and the first gradient ascent network uses the XGboost gradient ascent network.
[0111] Step S31. Calculate the input of each gene after preprocessing into the variable screening model, according to the output of the gene site weight, calculate the exon polygenic risk score of the gene in the, including:
[0112] Preprocess each gene, input the preprocessed multiple gene sites into the snpnet algorithm model, and perform lasso batch screening iterative regression analysis;
[0113] Add covariants in the regression analysis, and do not lasso scale the parameters of the covariants, to output the gene site weight removed from linkage disequilibrium.
[0114] Specifically, a bed file of the exon range of all coding protein genes in GENECODE v29 is obtained, for each gene, the bgenix and plink are used to screen the variations in the exon region, and the variations with INFO parameter <0.8, hardy-wenberg test p value less than 10 -6 , the total number of variation genotypes is less than 10, and the individual losing more than 50% of the variations in the region is removed. The above screened gene sites are input into the snpnet algorithm for batch screening iterative lasso regression analysis. Among them, batch screening iterative lasso is a multiple iteration algorithm, and the initial value is an empty set as the candidate variable set, and the following steps are repeated in multiple iterations:
[0115] 1) Calculate the inner product of each variable and the residual, and gradually find the variable that explains the most residual (i.e. the variable with the largest inner product with the residual) and add it to the candidate variable set.
[0116] 2) Solve LASSO on it, that is, find the parameter λ and the regression coefficient :
[0117] ; (8)
[0118] Wherein, y and x represent the phenotype to be predicted and the parameters used for prediction, β represents the regression coefficient, λ represents the scaling coefficient,
[0119] 3) Find the smallest parameter λ and its corresponding set of regression coefficients that meet the following boundary conditions among all solutions :
[0120] ; (9)
[0121] Wherein, represents in the kth iteration, The corresponding residual is x, and n represents the total number of parameters used for prediction.
[0122] In the LASSO regression model, gender, age, sample source, genotype principal component and other user-provided covariates are added, and the regression parameters of these covariates are not scaled by LASSO. According to the output gene site weight, the polygenic risk score of all individuals on the gene exon can be calculated.
[0123] Step S32. Input the exon polygenic risk score and all representative parameters in the regulatory region of the gene into the first gradient ascent network to obtain the phenotype prediction value of the gene, including:
[0124] 1) According to the functional genomics annotation, the regulatory region of each gene is preset;
[0125] Specifically, the enhancer-gene relationship annotation of 222 human tissues provided by epimap is obtained. For each tissue, at least the enhancers marked as associated in half of the samples are taken and the union is taken. For each gene, the enhancers associated with it in all tissues are taken and the union is taken, and then the upstream and downstream 100,000bp range of the gene is added as the entire regulatory region of the gene. Then, for each gene, all representative parameters in its regulatory region are taken into the next step.
[0126] 2) The exon polygenic risk score and all representative parameters of the gene in the regulatory region are input as a model, and the residual of the regression model of the trait to be predicted is used as the model prediction label. The first gradient ascent network is used to predict the phenotype prediction value of the gene.
[0127] Specifically, for each gene, the above exon polygenic risk score and all representative parameters in the regulatory region are input into XGboost to establish a gradient boosting model (R language xgboost function). XGboost uses independent variables to build a tree predictor, and combines existing trees and their prediction residuals in iterations to achieve the effect of gradient ascent. In the iteration, the tree generated in the mth iteration follows the following rules:
[0128] ; (10)
[0129] wherein, is the regularization coefficient of the mth iteration, is the residual of the (m-1)th iteration, is the function obtained by training to predict the (m-1)th iteration residual. In this embodiment, the result variable is selected as the residual of the following regression model:
[0130] ; (11)
[0131] wherein, , , , , , respectively represent height, age, gender, sample source, genetic principal component and other user-specified covariates.
[0132] Preferably, the XGboost selects the regression: squarederror mode, the eta value is set to 1, the gamma value is set to 1, and the maximum depth is set to 3. The XGboost calculates the height prediction value of each gene on each individual. After deleting all genes with empty output models, the prediction value of each gene is respectively correlated with the residual of the to-be-predicted trait as the result variable. Preferably, the prediction values of the top 5000 genes with the smallest p values are selected to enter the next step of analysis.
[0133] Step S4. The phenotype prediction value of each gene is simultaneously input into a second gradient ascent network for integration to generate a final prediction result.
[0134] In step S4, the phenotype prediction value of each gene is simultaneously input into a second gradient ascent network for integration to generate a final prediction result, including:
[0135] The phenotype prediction value of each gene is input as a model, and the true value of the to-be-predicted trait is input as a prediction label, and the final prediction result of the individual phenotype is obtained by using the second gradient ascent network.
[0136] Specifically, the height prediction value of the 5000 genes selected in the above step is input into a new XGboost model, and the label is the original value of the to-be-predicted phenotype. If the user provides the covariates that need to be corrected in the previous step , and these covariates are determined by expert knowledge to be used for phenotype prediction, they are also included in the XGboost. Since the height is a continuous variable, the XGboost selects the regression: squarederror mode, the eta value is set to 1, the gamma value is set to 1, and the maximum depth is set to 5. Since the number of variables in this step is large, the calculation demand is large, and in the case of insufficient user computing resources, the maximum depth can be adjusted to 3.
[0137] Preferably, after using the model for prediction, the method further includes model verification and application, specifically as follows:
[0138] In the independent validation population, repeat the above steps S1-S3 to obtain the gene level prediction measure, and input into the XGboost model obtained in the last step. If the standardization of the phenotype to be predicted is not performed during the training in step S4, the value output by the XGboost is the same as the prediction value of the real phenotype unit, which may be affected by the systematic bias between different populations. If the phenotype to be predicted is standardized during the training, the output value is the Z value, which needs to be converted into the real value according to the general distribution of the height of the target population, that is,
[0139] ; (12)
[0140] Based on the above disclosure, the embodiment does not directly use genotype data, but converts genotype (haplotype) into functional genomic parameters through a deep neural network, solving the defects that the existing technology cannot locate the effective site and accurately predict the disease according to the effective site when directly analyzing the correlation between genotype and disease. Since the complex interaction relationship between genes is considered, the representative parameters are selected from the functional genomic parameters, and the correlation between each region and the individual phenotype and the prediction value of the single gene are analyzed by using the representative parameters of the region, thereby integrating the information from each functional genomic change at the single gene level, avoiding the limitation of linear correlation of isolated analysis of each gene site. Finally, the prediction value of all genes is integrated through the gradient ascent network to fully reflect the nonlinear relationship between genes, thereby improving the prediction accuracy and having important significance for the prediction of various individual phenotypes.
[0141] In a second aspect, the present application provides a device for predicting individual phenotypes based on human whole genome genotypes, comprising:
[0142] A parameter conversion module is configured to obtain haplotypes in all regions of the whole genome of each individual, convert the haplotypes into functional genomic parameters based on a deep neural network, and select at most one representative parameter from the functional genomic parameters of each region.
[0143] An association analysis module is configured to quantitatively analyze the correlation between each region and the individual phenotype by using a regression model according to the representative parameters of the region.
[0144] A first prediction module is configured to predict the phenotype prediction value of each gene by combining a variable screening model and a first gradient ascent network according to the risk score of each gene and the representative parameter of the regulatory region where the gene is located.
[0145] A second prediction module is configured to integrate each gene phenotype prediction value into a second gradient ascent network to generate a final prediction result.
[0146] The working process, working details and technical effects of the foregoing device provided by the second aspect of the embodiment can be referred to the method as described in the first aspect or any possible design in the first aspect, and will not be described here again.
[0147] In a third aspect, the present application provides a computer device, comprising a memory, a processor and a transceiver connected in sequence and in communication, wherein the memory is configured to store a computer program, the transceiver is configured to transceive messages, and the processor is configured to read the computer program and execute the method as described in any possible design of the first aspect.
[0148] For example, the memory can include, but is not limited to, a random access memory (RAM), a read-only memory (ROM), a flash memory, a first input first output (FIFO) memory and / or a first input last output (FILO) memory, etc.; the processor can be, but is not limited to, a microprocessor of STM32F105 series; the transceiver can be, but is not limited to, a WiFi wireless transceiver, a Bluetooth wireless transceiver, a GPRS wireless transceiver and / or a ZigBee wireless transceiver, etc. In addition, the computer device can further include, but is not limited to, a power module, a display screen and other necessary components.
[0149] The working process, working details and technical effects of the foregoing computer device provided by the third aspect of the embodiment can be referred to the method as described in the first aspect or any possible design in the first aspect, and will not be described here again.
[0150] In a fourth aspect, the present application provides a computer readable storage medium, wherein instructions are stored on the computer readable storage medium, and when the instructions are run on a computer, the method as described in any possible design of the first aspect is executed.
[0151] The computer readable storage medium refers to a carrier for storing data, which can include, but is not limited to, a floppy disk, an optical disk, a hard disk, a flash memory, a USB flash disk and / or a memory stick, etc., and the computer can be a general-purpose computer, a special-purpose computer, a computer network or other programmable devices.
[0152] The working process, working details and technical effects of the aforementioned computer readable storage medium provided by the fourth aspect of the embodiment can be referred to the method as described in the first aspect or any possible design in the first aspect, and will not be repeated here.
[0153] In a fifth aspect, the present application provides a computer program product containing instructions which, when executed on a computer, cause the computer to carry out the method as described in any possible design of the first aspect.
[0154] The working process, working details and technical effects of the aforementioned computer program product containing instructions provided by the fifth aspect of the embodiment can be referred to the method as described in the first aspect or any possible design in the first aspect, and will not be repeated here.
[0155] Finally, it should be noted that: the above only for the preferred embodiments of the present application, and not for limiting the scope of protection of the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the scope of protection of the present application.
Claims
1. A method for predicting individual phenotypes based on human whole-genome genotype, characterized in that, include: Obtain haplotypes from all regions of the whole genome of each individual, convert haplotypes into functional genomic parameters based on deep neural networks, and select at most one representative parameter from the functional genomic parameters of each region; Based on representative parameters of each region, regression models are used to quantitatively analyze the association between each region and individual phenotypes, including: Within each region representing a parameter, the association between that region and the individual phenotype is quantitatively analyzed based on a linear regression model. The effective region is selected from areas with internal chain imbalances using a fine-grained localization algorithm; Obtain the biological significance annotation for each region. Quantitatively analyze, based on a generalized linear regression model, whether the regions influencing individual phenotypes are more closely related to a particular biological significance annotation. For each annotation, construct the following generalized regression analysis model to analyze the following question: whether regions influencing phenotypes are more closely related to this annotation; the expression of the generalized regression analysis model is as follows: ；（7） in, Indicate whether the region conforms to the annotation (1 / 0). This indicates the posterior probability calculated by SUSIE that the region is the effective region. This indicates the total number of haplotypes in the region. This indicates the total number of samples included in the analysis for that region. This represents the sequence category score corresponding to the reference sequence in this region. This indicates the frequency of chromosome recombination within that region. This indicates the intensity of background selection in human history for that region. This indicates the degree of conservatism in the region. This represents the block number as a random effect term, where, in the regression results, The regression coefficients and p-values of the terms are used to measure the association between the heritability of a trait and a given annotation; Based on the risk score of each gene and the representative parameters of the regulatory region in which the gene is located, the phenotypic predicted value of each gene is obtained by combining the variable screening model and the first gradient ascent network prediction, including: After preprocessing each gene, it is input into the variable selection model for calculation. Based on the output gene locus weights, the exon polygenic risk score of that gene is calculated. The calculation process for each gene after preprocessing includes: Each gene is preprocessed, and the preprocessed gene loci are input into the SNPNet algorithm model for Lasso batch screening and iterative regression analysis. Incorporate covariates into the regression analysis without lasso scaling of the covariate parameters to output gene locus weights that remove linkage disequilibrium. The exon polygenic risk score and all representative parameters of the regulatory region containing the gene are input into the first gradient ascent network to obtain the phenotypic prediction value of the gene, including: Based on functional genome annotation, the regulatory regions of each gene are predefined; The exon polygenic risk score and all representative parameters of the gene in the regulatory region are used as model inputs, and the residuals of the regression model of the trait to be predicted are used as model prediction labels. The phenotypic prediction value of the gene is obtained by using the first gradient ascent network. The phenotypic predictions for each gene are simultaneously input into a second gradient ascent network for integration, generating the final prediction results, including: The phenotypic prediction value of each gene is used as the model input, and the actual value of the trait to be predicted is used as the prediction label. The final prediction result of the individual phenotypic is obtained by using the second gradient ascent network. Before selecting at most one representative parameter from multiple functional genomic parameters for each region based on a fine-mapping algorithm, the process also includes: Multiple functional genomic parameters are integrated to obtain multiple category scores for each base sequence, and these multiple category scores are then processed. Each category score after processing is transformed according to the additive, recessive, and dominant genetic models to generate corresponding mediating parameters, so that the fine-mapping algorithm can select at most one representative parameter from the mediating parameters of each region.
2. The method according to claim 1, characterized in that, Haplotypes are obtained from all regions of the whole genome of each individual. These haplotypes are then converted into functional genomic parameters using a deep neural network. At most one representative parameter is selected from the functional genomic parameters of each region, including: The entire genome of each individual is divided into regions to obtain the haplotype and corresponding base sequence of each region; Each base sequence is converted into multiple corresponding functional genomic parameters based on deep neural networks; Based on the fine-mapping algorithm, at most one representative parameter is selected from multiple functional genomic parameters in each region.
3. The method according to claim 2, characterized in that, The entire genome of each individual is segmented into regions to obtain the haplotype and corresponding base sequence for each region, including: The entire genome of each individual is divided into multiple regions according to a preset length, and each region is guaranteed to cover the functional elements of the entire genome as completely as possible. Genotypes are screened based on non-allelic frequencies, the screened genotypes are converted into haplotypes, and the base sequence corresponding to each haplotype is obtained.
4. A computer device, characterized in that, The device includes a memory, a processor, and a transceiver that are sequentially and communicatively connected, wherein the memory is used to store a computer program, the transceiver is used to send and receive messages, and the processor is used to read the computer program and execute the method as described in any one of claims 1-3.