Transformers-based genomic prediction method

Abstract

The application provides a kind of based on the genome prediction method of Transformer.The method comprises the following steps: step 1: obtaining original SNP site data set;Step 2: using stacking strategy to carry out feature selection to original SNP site data set to obtain feature subset step 3: all SNP sites in the feature subset are encoded to obtain an encoding matrix;Step 4: the encoding matrix is input into the trained phenotype prediction model to obtain the predicted phenotype value;Wherein, the phenotype prediction model comprises a Transformer network and a regression prediction module;Wherein, the Transformer network obtains the potential association information between SNP sites according to the input encoding matrix;The regression prediction module predicts the phenotype value according to the potential association information between SNP sites.The application improves the phenotype prediction accuracy while reducing the computational cost through feature selection.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, and in particular to a genome prediction method based on Transformer. Background Technology

[0002] Genomic prediction (GP) utilizes high-throughput genotyping technologies such as genotyping sequencing to generate genes that can significantly shorten generation intervals, thereby producing higher genetic effects. Single nucleotide polymorphisms (SNPs) are very common genetic variations in the genome and are considered to be closely related to the occurrence of traits and diseases. Genomic prediction uses SNPs to predict reproductive value (BV) or quantitative phenotype.

[0003] Common GP methods include Bayesian and Genomic Best Linear Unbiased Prediction (GBLUP), both of which are statistical learning methods. These linear models typically only consider additive inheritance, ignoring the complex nonlinear relationships that may exist between markers and phenotypes. They have limited flexibility in handling nonlinear effects in high-dimensional genomic data, resulting in high computational costs.

[0004] In addition, high-throughput data generated by genotyping is usually high-dimensional. If phenotypes are predicted directly based on high-dimensional gene sequences, not only will noise interference lead to inaccurate phenotype predictions, but there will also be a problem of large computational load. Summary of the Invention

[0005] In order to reduce computational costs and improve the accuracy of genome phenotype prediction, this invention provides a Transformer-based genome prediction method.

[0006] This invention provides a Transformer-based genome prediction method, comprising:

[0007] Step 1: Obtain the original SNP site dataset;

[0008] Step 2: Use a stacking strategy to perform feature selection on the original SNP site dataset to obtain a feature subset. ;

[0009] Step 3: For the feature subset All SNP sites are encoded to obtain the encoding matrix;

[0010] Step 4: Input the encoding matrix into the trained phenotypic prediction model to predict the phenotypic value; wherein, the phenotypic prediction model includes a Transformer network and a regression prediction module; wherein, the Transformer network obtains the potential association information between SNP sites based on the input encoding matrix; the regression prediction module predicts the phenotypic value based on the potential association information between SNP sites.

[0011] Furthermore, step 2 specifically includes:

[0012] The XGBoost algorithm, rrBLUP algorithm, and mutual information method were used to select features from the original SNP dataset. The top k SNPs from each algorithm were then merged to obtain the final feature subset. .

[0013] Furthermore, the Transformer network includes three parallel Transformer modules, each containing four attention heads; and different weights are assigned to the outputs of the three parallel Transformer modules, and the weighted sum of the outputs of the three parallel Transformer modules is used as the final output of the Transformer network.

[0014] Furthermore, the regression prediction module includes a Dense layer and a Dropout layer.

[0015] Furthermore, when training the phenotypic prediction model, the Huber loss function is used to calculate the loss between the true phenotypic value and the predicted phenotypic value.

[0016] Furthermore, step 3 specifically includes:

[0017] Step 3.1: Define all SNP sites to form an SNP sequence Using one-hot encoded sequences The different alleles at each SNP locus are encoded to generate a matrix. ;in, This indicates the total number of SNP sites. Indicates the feature dimension of SNP sites;

[0018] Furthermore, the feature subset is divided according to the following formula. All SNP sites are positionally encoded;

[0019]

[0020] in, , These are the position indexes. The 2i-th and 2i+1-th components of the encoding vector, ;

[0021] Step 3.2: Add positional encoding to the one-hot encoded genomic data to obtain the feature representation of the encoded genome. :

[0022] .

[0023] The beneficial effects of this invention are:

[0024] The method of this invention first uses a stacking strategy to select important SNP sites for feature selection, and then captures the location information of SNP sites through a combination of positional encoding and one-hot encoding. Next, a Transformer architecture is used to further learn the nonlinear interactions between SNP sites, thereby predicting phenotypes through the genome.

[0025] Furthermore, the Transformer architecture contains multiple Transformer modules, and the outputs of these modules have different weights, which enhances the robustness of the model.

[0026] Experiments show that the present invention is simple and effective, and compared with other existing methods, the present invention has better performance in predicting genome phenotypes. Attached Figure Description

[0027] Figure 1 A flowchart of the Transformer-based genome prediction method provided in this embodiment of the invention;

[0028] Figure 2 A schematic diagram of feature selection using a stacking strategy provided in an embodiment of the present invention;

[0029] Figure 3 This is a schematic diagram of a Transformer network learning potential associations provided in an embodiment of the present invention;

[0030] Figure 4 This is a schematic diagram of the regression prediction module provided in an embodiment of the present invention;

[0031] Figure 5 This diagram illustrates the impact of different feature selection numbers on prediction results, as provided in an embodiment of the present invention. Detailed Implementation

[0032] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of the embodiments of this invention will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0033] Example 1

[0034] like Figure 1 As shown, this embodiment of the invention provides a Transformer-based genome prediction method, including the following steps:

[0035] S101: Obtain the original SNP site dataset;

[0036] S102: A feature subset is obtained by using a stacking strategy to select features from the original SNP site dataset. ;

[0037] Specifically, such as Figure 2 As shown, in this embodiment, the XGBoost algorithm, the rrBLUP algorithm, and the mutual information method are used to select features from the original SNP site dataset. The top k SNP sites from the three feature selection algorithms are merged to obtain the final feature subset. .

[0038] The core idea of ​​the stacking strategy is to combine the predictions of multiple base learners to obtain the final prediction. It requires training multiple base learners on the same training dataset, and then feeding the predictions of the base learners into a higher-level meta-model to make the final prediction.

[0039] XGBoost (eXtreme Gradient Boosting) is a machine learning algorithm based on gradient boosting trees. In this embodiment of the invention, this algorithm is used to evaluate the importance of each SNP site by analyzing the splitting process of decision tree nodes and the usage of each SNP site. The basic form of the algorithm's loss function consists of an empirical loss term and a regularization term. Loss Function The expression is:

[0040]

[0041]

[0042] in, This represents the experience loss term, which represents the true phenotypic value. and predicted phenotypic values The error between; Let represent the regularization term, and let represent the sum of the complexities of all t trees. This is also the method used in XGBoost to control overfitting. Represents the i-th sample The phenotypic prediction after step t. Represents the i-th sample Phenotypic predictions after step t-1. Represents the t-th tree For the i-th sample Phenotypic predicted values.

[0043] rrBLUP (ridge regression best linear unbiased prediction) is a statistical method for feature selection. This algorithm combines the ideas of ridge regression and linear best unbiased prediction. In this embodiment of the invention, the algorithm is used to estimate the relationship between random effects and fixed effects for feature selection.

[0044]

[0045] in, Indicates phenotypic value, yes The genotype matrix, where n is the number of samples and m is the total number of SNP loci. It is the mean of all phenotypic values ​​y. yes The label association effect and Follows a normal distribution , This indicates a random effect.

[0046] Mutual information (MIS) algorithms are used to measure the correlation and dependency between two variables. In this embodiment of the invention, the MIS method is used to calculate the correlation between genotype and phenotype variables. In this embodiment, SNP sites with higher MIS values ​​are considered more correlated with phenotype values, thereby selecting the most relevant subset of features. The MIS formula is as follows:

[0047]

[0048] in, This is the joint probability distribution of the genotype variable X and the phenotypic variable Y. It should be noted that when X and Y are statistically independent, the MI is zero, i.e.: .

[0049] In this embodiment, the importance scores of each SNP site are calculated using the three feature selection methods described above, and then sorted in descending order to obtain the top k SNP sites for each feature selection method. Finally, a union operation is performed on all the top k SNP sites to obtain the feature subset. The above process can be represented by the following formula:

[0050]

[0051] in, This represents the original SNP locus dataset. For example, when k is 200, it means that the first 200 SNP loci are selected for each method.

[0052] S103: For the feature subset All SNP sites are encoded to obtain an encoding matrix, which is then used as input for subsequent models.

[0053] As one possible implementation method, this embodiment of the invention uses a combination of one-hot encoding and positional encoding to encode all SNP sites. The specific process is as follows:

[0054] S1031: To represent the locus information of genomic data, one-hot encoding is used to encode the different alleles of each SNP. A length of [missing information] is defined. SNP sequences of the sample After one-hot encoding, a matrix is ​​generated. Representing SNP sequences .in This indicates the size of the sample, which is the total number of SNP loci. The feature dimension of the SNP site is represented, that is, the one-hot encoding dimension. In this embodiment of the invention, the SNP is defined to have four values: AA, Aa, aa and missing value, so d=4.

[0055] Furthermore, the feature subset is divided according to the following formula. All SNP sites are positionally encoded;

[0056]

[0057] in, , These are the position indexes. The 2i-th and 2i+1-th components of the encoding vector, .

[0058] Positional encoding effectively captures positional information of a sequence for the model by assigning a unique vector representation to each position in the sequence, thereby improving modeling performance.

[0059] S1032: The feature representation of the encoded genome is obtained by adding positional coding to the genomic data after one-heat coding. :

[0060]

[0061] S104: Convert the encoding matrix (i.e., feature representation) The input is fed into a trained phenotypic prediction model to predict phenotypic values. The phenotypic prediction model includes a Transformer network (also known as Transformer Blocks) and a regression prediction module. The Transformer network obtains potential association information between SNP sites based on the input encoding matrix. The regression prediction module predicts phenotypic values ​​based on the potential association information between SNP sites.

[0062] Specifically, the self-attention network Transformer completely abandons recursion and convolution, and can compute the relationship between any two positions in the input sequence in a single step. Therefore, it can focus on the most influential genomic regions highly correlated with traits, more effectively capturing non-local, long-term dependencies in the genome, which helps extract more relevant breeding signals. Furthermore, this architecture allows for parameter sharing, enabling the detection of significant gene variations with different effects, thus allowing for better modeling of complex nonlinear interactions in genotype-phenotype relationships.

[0063] As one possible implementation method, such as Figure 3 As shown, the Transformer network in this embodiment includes three parallel Transformer modules, each containing four attention heads; and different weights are assigned to the outputs of the three parallel Transformer modules; finally, the fused feature matrix is ​​obtained through the Gelu activation function.

[0064] The feature extraction module uses a Transformer network with an embedding layer dimension of 128, the activation function being Gelu, 3 Transformer blocks, 4 attention heads, and a dropout ratio of 0.3. During the training and prediction phases of the regression prediction module, [the following is used:] The activation function is used in conjunction with the Huber loss function.

[0065] Specifically, the implementation of the Transformer network is as follows:

[0066] Multi-head attention layer: This layer is the main mechanism of the Transformer module to capture long-range dependencies. It uses multiple parallel self-attention heads to extract features and finally fuses the outputs of multiple parallel self-attention heads.

[0067] The process of extracting features for each self-attention head is as follows: The features of the encoded genome are represented... With the learnable weight matrix respectively , and Multiplication yields the query matrix Key matrix Sum matrix Then, the self-attention score is calculated through the Softmax operation, and finally compared with the value matrix. Multiplication yields the representation of the self-attention head. :

[0068]

[0069]

[0070] in This represents the j-th self-attention head.

[0071] After passing through the multi-head attention layer, the output of the multi-head attention layer is added to the original input via a residual connection to obtain the feature representation. :

[0072]

[0073] in This represents a normalization operation. Layer normalization is performed after the multi-head attention layer to stabilize learning.

[0074] Adding residual connections between the inputs and outputs of multi-head attention layers and feedforward layers helps alleviate the vanishing gradient problem during training. Therefore, the feature representation of the final output of a Transformer block... as follows:

[0075]

[0076] The feedforward neural network can be expressed by the following formula:

[0077]

[0078] Unlike a single Transformer block, merging three Transformer blocks allows for simultaneous attention to multiple aspects of the sequence, thus enriching the feature representation. The feature representation is obtained by merging the outputs of the three Transformer blocks through a higher-level connection. :

[0079]

[0080] in This represents the feature representation after passing through the Transformer network. , and All of these are hyperparameters.

[0081] As one possible implementation method, such as Figure 4 As shown, the regression prediction module in this embodiment includes a Dense layer and a Dropout layer.

[0082] Specifically, the Dense layer is used to reduce the dimensionality of the Transformer network's output. Specifically, the feature vector output by the Transformer network serves as the input to the first linear layer in the fully connected layers, and the output of the last linear layer is the phenotypic prediction, which is used as the genome prediction result. The specific formula is:

[0083]

[0084] in, Is with Binary masks with the same shape This represents element-wise multiplication. This represents the predicted phenotypic value.

[0085] Furthermore, to address the overfitting problem associated with stacked multi-layer networks, the Gelu activation function is used after the first linear layer, and Dropout layers are added between subsequent linear layers. The Dropout layers remove neurons from the network with a certain probability. For stochastic gradient descent, the introduction of random neuron removal ensures that each iteration trains a different network. The Dropout layers effectively solve the overfitting problem and improve the model's generalization ability.

[0086] It is understandable that the genome prediction model is pre-trained. In this embodiment of the invention, the Adamoptimizer stochastic optimization algorithm is used for learning. Meanwhile, considering that phenotypic data often exhibits imbalance, i.e., the distribution of phenotypic data mostly follows a normal Gaussian distribution, with fewer phenotypic data at the extremes and more in the middle, such a data distribution leads to poor model performance on phenotypic data with low data volume. To address this, this embodiment of the invention uses the Huber loss function to calculate the loss between the true phenotypic value and the predicted phenotypic value during model training. This is because, compared to squared error loss, Huber loss is not sensitive to outliers in the data and is differentiable at 0. Essentially, Huber loss is absolute error, but it becomes squared error when the error is very small. How small the error becomes quadratic error is controlled by the hyperparameter δ (delta). When the Huber loss is between [0-δ, 0+δ], it is equivalent to MSE, while in [-∞, δ] and [δ, +∞] it is MAE. The Huber loss is formalized as follows:

[0087]

[0088] in These are the actual phenotypic values. These are predicted phenotypic values. These are hyperparameters. When |y−f(x)| > δ, the gradient of the Huber loss function is always approximately δ, which ensures that the model updates parameters at a relatively fast speed; when |y−f(x)| ≤ δ, the gradient gradually decreases, which ensures that the model obtains the global optimum more accurately.

[0089] This invention employs a stacking strategy to integrate different feature selection methods to calculate the importance of each SNP site, and utilizes a Transformer network to obtain potential associations between SNP sites, thereby making full use of known biological information to accurately predict phenotypes.

[0090] To verify the effectiveness of the method of the present invention, the present invention also provides the following experimental data.

[0091] (a) Dataset

[0092] Soybean Dataset: The soybean dataset used is from the SoyNAM project, which contains over 5000 recombinant inbred lines and 4236 common SNPs between the estimated data and the original quality assurance data. Genotypic and phenotypic data are provided in the “SoyNAM” R package. Phenotypic data were selected from locations in Illinois in 2013 and 2012, focusing on five characteristics: plant height, oil content, water content, protein content, and yield.

[0093] Winter wheat dataset: Ten winter wheat breeding populations adapted to the Pacific Northwest of the United States were used, including one association map (AMP), two F5 biparental populations, two double haploid (DH) biparental populations, and five winter wheat line populations used as yield validation populations in preliminary yield trials (Prel). The AMP comprised 456 lines, totaling 1768 samples and 11089 SNPs, evaluated between 2015 and 2018 in Lind (LND) and Pullman (PUL), Washington, USA.

[0094] Table 1. Basic information about the soybean and wheat datasets.

[0095]

[0096] Crop phenotypic information is normalized by performing min-max normalization on the resulting matrix using the following formula, ultimately yielding the normalized matrix. The formula is as follows:

[0097]

[0098] in, This indicates the sample size for that crop. A numerical value representing a certain phenotypic trait. This represents the minimum value of this phenotypic trait. This represents the maximum value of this phenotypic trait.

[0099] (II) Evaluation Indicators

[0100] The predictive performance of the method of the present invention was tested by conducting ten-fold cross-validation experiments.

[0101] (1) Ten-fold cross-validation

[0102] In each fold, the phenotypic samples in the dataset are divided into (training set: validation set: test set) = (9:1):1. The average of the 10-fold cross-validation results is taken after each 10-fold cross-validation. Finally, ten 10-fold cross-validations are performed, and the average is taken to evaluate the model's performance, ensuring that the error in the experimental results is minimized.

[0103] (2) Evaluation indicators

[0104] In regression prediction problems, in regression analysis, the mean squared error (MSE), mean absolute error (MAE), and R-squared (R²) are used. The Pearson correlation coefficient (Pcc) and other metrics are used to evaluate the model's performance.

[0105]

[0106]

[0107]

[0108]

[0109] in This represents the number of samples in the test set. Indicates the predicted phenotypic value. Represents the true phenotypic value.

[0110] (III) Experimental Results

[0111] The parameter settings for the method of the present invention are shown in Table 2.

[0112] Table 2 Parameter settings for the method of the present invention

[0113]

[0114] (1) Tenfold cross-validation analysis

[0115] Comparative experiments were conducted on the dataset. All experiments underwent 10-fold cross-validation ten times, and the mean was taken to ensure fairness. Specific experimental results are listed in the table below:

[0116] We conducted experiments using rrBLUP, Bayes Lasso, XGBoost, RandomForest, MLP, dualCNN, and GPFTSS, respectively. The experimental results are as follows:

[0117] In the winter wheat dataset, yield traits showed the highest Pearson correlation coefficient in the GPFTSS model. Deep learning models (MLP, dualCNN) performed worse than statistical learning models on the winter wheat dataset. This may be due to the architecture of deep learning models; for example, the dualCNN model still uses the convolutional kernel parameters of the original model, which may differ from the intrinsic characteristics of the winter wheat dataset, thus resulting in poorer performance.

[0118] Table 3 Prediction Indicators of Different Models for Winter Wheat

[0119]

[0120] On the winter wheat dataset, the MLP and dualCNN methods performed the worst, while the rrBLUP and Bayesian methods yielded comparable results. The Pcc of the method in this invention was 0.87. The value is 0.74. In the ensemble model, it can be seen that the method of this invention, whether in terms of Pcc or the explanatory coefficients (…),… All of these are optimal.

[0121] For plant height, oil content, water content, protein, and yield traits, the rrBLUP, BRR, and Bayes A methods performed poorly, while MLP, dualCNN, and XGBoost yielded the next best results, and GPFTSS yielded the highest results.

[0122] Table 4. Predictions (Pcc) for various traits using different models on the soybean dataset.

[0123]

[0124] (2) The effect of the number of features selected on the experimental results

[0125] Table 5. Impact of different SNP counts on prediction accuracy in the winter wheat dataset.

[0126]

[0127] Table 6. Impact of different SNP numbers on prediction accuracy in soybean dataset.

[0128]

[0129] The results showed that as the number of SNP sites increased, the phenotypic prediction performance of both datasets first increased and then decreased, with the best performance observed when the number of SNP sites was 500. Figure 5 As shown.

[0130] (3) Ablation test

[0131] Table 7. Impact of Three Feature Selection Methods on Model Predictive Ability

[0132]

[0133] By using different feature selection methods, the results show that using a stacking strategy for feature selection performs better in subsequent Transformer models than using XGBoost, Mutual, and rrBLUP individually.

[0134] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A Transformer-based genome prediction method, characterized in that, include: Step 1: Obtain the original SNP site dataset; Step 2: Use a stacking strategy to perform feature selection on the original SNP site dataset to obtain a feature subset. Specifically, this includes: using the XGBoost algorithm, the rrBLUP algorithm, and the mutual information method to perform feature selection on the original SNP site dataset, and merging the top k SNP sites from the three feature selection algorithms to obtain the final feature subset. ; The XGBoost algorithm evaluates the importance of each SNP site by analyzing the splitting process of decision tree nodes and the usage of each SNP site; the loss function of the XGBoost algorithm... The expression is: in, This represents the experience loss term, which represents the true phenotypic value. and predicted phenotypic values The error between; Let represent the regularization term, and let represent the sum of the time complexity of all t trees. Represents the i-th sample The phenotypic prediction after step t. Represents the i-th sample Phenotypic predictions after step t-1. Represents the t-th tree For the i-th sample Phenotypic predicted values; Feature selection is performed by estimating the relationship between random and fixed effects using the rrBLUP algorithm. The calculation formula is as follows: in, Indicates phenotypic value, yes The genotype matrix, where n is the number of samples and m is the total number of SNP loci. It is the mean of all phenotypic values ​​y. yes The label association effect and Follows a normal distribution , Indicates random effects; Mutual information is used to calculate the correlation between genotype and phenotypic variables. SNPs with higher mutual information values ​​are considered more correlated with phenotypic values, thus selecting the most relevant subset of features. The mutual information formula is: in, It is the joint probability distribution of genotype variable X and phenotype variable Y; when X and Y are statistically independent, ; After sorting the SNP sites in descending order of importance scores to obtain the top k SNP sites for each feature selection method, a union operation is performed on all the top k SNP sites to obtain the feature subset. The above process can be expressed by the following formula: in, This represents the original SNP site dataset; Step 3: For the feature subset All SNP sites are encoded to obtain the encoding matrix; Step 4: Input the encoding matrix into the trained phenotypic prediction model to predict phenotypic values; wherein, the phenotypic prediction model includes a Transformer network and a regression prediction module; wherein, the Transformer network obtains the potential association information between SNP sites based on the input encoding matrix; and the regression prediction module predicts phenotypic values ​​based on the potential association information between SNP sites.

2. The genome prediction method based on Transformer according to claim 1, characterized in that, The Transformer network consists of three parallel Transformer modules, each containing four attention heads. Different weights are assigned to the outputs of the three parallel Transformer modules, and the weighted sum of the outputs of the three parallel Transformer modules is used as the final output of the Transformer network.

3. The genome prediction method based on Transformer according to claim 1, characterized in that, The regression prediction module includes a Dense layer and a Dropout layer.

4. The genome prediction method based on Transformer according to claim 1, characterized in that, When training the phenotypic prediction model, the Huber loss function is used to calculate the loss between the true phenotypic value and the predicted phenotypic value.

5. The genome prediction method based on Transformer according to claim 1, characterized in that, Step 3 specifically includes: Step 3.1: Define all SNP sites to form an SNP sequence Using one-hot encoded sequences The different alleles at each SNP locus are encoded to generate a matrix. ;in, This indicates the total number of SNP sites. Indicates the feature dimension of SNP sites; Furthermore, the feature subset is divided according to the following formula. All SNP sites are positionally encoded; in, , These are the position indexes. The 2i-th and 2i+1-th components of the encoding vector, ; Step 3.2: Add positional encoding to the one-hot encoded genomic data to obtain the feature representation of the encoded genome. : 。

Images