A wheat genotype-to-phenotype prediction method based on heterogeneous deep learning

Abstract

The application provides a wheat genotype-to-phenotype prediction method based on heterogeneous deep learning, and belongs to the technical field of crop molecular breeding and artificial intelligence. The wheat genotype-to-phenotype prediction method comprises the following steps: obtaining wheat genotype data to be tested, inputting the prediction model to obtain a phenotype prediction result. A heterogeneous feature extraction network set is used to extract differential heterogeneous features from the genotype data, including a first feature representing additive genetic effects and a second feature representing gene nonlinear interaction effects, so as to accurately match the genetic expression characteristics of wheat. Then, an adaptive feature fusion network is used to dynamically weight and fuse various heterogeneous features based on learnable parameters, obtain fusion features with high representation capability, and finally output the phenotype result through a prediction network. The wheat genotype-to-phenotype prediction method solves the defects of traditional prediction models that cannot simultaneously consider the linear additive genetic effects of wheat and the local and global nonlinear interaction effects of genes, improves the prediction accuracy and generalization ability of the genotype-to-phenotype prediction, and provides efficient and accurate technical support for wheat molecular assisted breeding and trait prediction.

Description

Technical Field

[0001] This invention belongs to the field of interdisciplinary technology of crop molecular breeding and artificial intelligence, specifically involving a method for predicting wheat genotype to phenotype based on heterogeneous deep learning. Background Technology

[0002] With the rapid development of high-throughput sequencing and gene chip technologies, obtaining high-density single nucleotide polymorphism (SNP) marker information across the entire crop genome has become efficient and cost-effective, providing a massive data foundation for genome selection. Genome selection, by learning the complex mapping relationship between genotypes and phenotypes in a training population, can accurately predict the phenotype of a candidate population based solely on its genotype data, and has become one of the key technologies in modern precision agriculture and molecular design breeding.

[0003] In crop breeding, many important agronomic traits are complex quantitative traits controlled by multiple genes with minor effects. Traditional pedigree-based optimal linear unbiased prediction methods have been gradually replaced by GBLUP methods based on genomic information. However, GBLUP and its various Bayesian derivative methods are essentially still linear models, and have inherent limitations in capturing nonlinear interaction effects between genotypes.

[0004] To address the aforementioned issues, existing technologies have explored the use of deep learning methods for genome prediction, primarily including the following approaches: The first type of approach is based on fully connected networks or multilayer perceptrons, using hundreds of thousands or even millions of SNP sites as input features and learning the interaction relationships between features through multiple fully connected layers. However, this approach is prone to overfitting when faced with ultra-high-dimensional SNP data, resulting in a significant drop in the model's prediction accuracy on independent validation sets, and the large number of parameters creates a heavy computational burden.

[0005] The second type of approach is based on convolutional neural networks, which use one-dimensional convolutional kernels to slide across SNP sequences to capture information about local linkage imbalance blocks. However, this approach focuses on local feature extraction and has limited ability to capture long-distance dependencies between genes that are geographically dispersed but functionally interact.

[0006] The third type of approach employs ensemble models, such as PAGER, which combines a one-dimensional convolutional neural network, a gated recurrent unit, and a multilayer perceptron to capture local features, long-range dependencies, and nonlinear fitting, respectively. However, this type of approach uses simple vector concatenation followed by implicit fusion into a fully connected layer, failing to adequately consider the varying contributions of features from different sources to the final phenotypic prediction. This may result in the loss of crucial information or the amplification of noise. Furthermore, most existing models are directly transferred from computer vision or natural language processing, without being specifically optimized for the limited training samples and high accuracy requirements of genome prediction.

[0007] In summary, existing technologies suffer from problems such as an imbalance in the model's capture of linear and nonlinear relationships, limitations of single network structures, and inefficiency of multi-source feature fusion strategies. There is an urgent need for a prediction method that can effectively integrate the linear and nonlinear effects of genotype data, comprehensively capture multi-scale features of the genome, and achieve efficient feature fusion. Summary of the Invention

[0008] In view of the above-mentioned defects or deficiencies in the existing technology, a method for predicting wheat genotype to phenotype based on heterogeneous deep learning is provided, including the following steps: Obtain the genotype data of the crop to be tested, input the genotype data into a pre-trained prediction model, and obtain the phenotypic prediction results; The prediction model obtains the phenotypic prediction results in the following manner: The heterogeneous feature extraction network set is used to extract multiple heterogeneous features with different biological significance from the genotype data. The multiple heterogeneous features include at least a first feature reflecting additive genetic effects and a second feature reflecting nonlinear interaction effects. By using an adaptive feature fusion network, multiple heterogeneous features are dynamically weighted and fused with learnable fusion parameters to obtain fused features; Based on the fusion features, the phenotypic prediction results are output through the prediction network.

[0009] According to the technical solution provided in this application, the extraction of multiple heterogeneous features from the genotype data using a heterogeneous feature extraction network set includes the following steps: The genotype data is simultaneously input into the first, second, and third parallel branches, respectively. In the first branch, the genotype data is linearly mapped, and a linear feature vector representing the additive genetic effect is extracted as the first feature; In the second branch, local pattern capture is performed on the genotype data to extract local nonlinear feature vectors that characterize local linkage disequilibrium block patterns, which are used as a component of the second feature. In the third branch, global context capture is performed on the genotype data to extract global context feature vectors that characterize the interaction relationships between distant sites across the entire genome, as another component of the second feature; The linear feature vector output from the first branch, the local nonlinear feature vector output from the second branch, and the global context feature vector output from the third branch are input into an adaptive feature fusion network, and the three are fused into a unified fusion feature vector through learnable fusion parameters.

[0010] According to the technical solution provided in this application, the dynamic weighted fusion of the multiple heterogeneous features through an adaptive feature fusion network includes the following steps: Assign a first learnable weight coefficient to the linear feature vector output by the first branch, assign a second learnable weight coefficient to the local nonlinear feature vector output by the second branch, and assign a third learnable weight coefficient to the global context feature vector output by the third branch. Multiply the linear feature vector by the first learnable weight coefficient to obtain a weighted linear feature vector; Multiply the local nonlinear feature vector by the second learnable weight coefficient to obtain a weighted local nonlinear feature vector; Multiply the global context feature vector by the third learnable weight coefficient to obtain the weighted global context feature vector; The fused feature vector is obtained by summing the weighted linear feature vector, the weighted local nonlinear feature vector, and the weighted global context feature vector element by element.

[0011] According to the technical solution provided in this application, the step of linearly mapping the genotype data and extracting linear feature vectors representing additive genetic effects in the first branch includes the following steps: The input genotype data is standardized to obtain a standardized genotype matrix; The standardized genotype matrix is ​​input into one or more stacked fully connected layers, wherein each fully connected layer does not have a non-linear activation function, and an initial linear feature vector is obtained through layer-by-layer linear transformation; The initial linear feature vector output by each fully connected layer is subjected to batch normalization and layer scaling in sequence. The batch normalization process is used to adjust the distribution of the initial linear feature vector to zero mean and unit variance. The layer scaling process is used to multiply the normalized feature vector with a learnable scaling factor to adjust the contribution of the layer output to the final linear feature vector. The feature vectors after batch normalization and layer scaling are used as the linear feature vectors output by the first branch.

[0012] According to the technical solution provided in this application, in the second branch, performing local pattern capture on the genotype data and extracting local nonlinear feature vectors representing local linkage disequilibrium block patterns includes the following steps: The genotype data is adjusted to a one-dimensional sequence format of single nucleotide polymorphism sites to obtain a one-dimensional input sequence; The one-dimensional input sequence is input into one or more stacked one-dimensional convolutional layers. In each one-dimensional convolutional layer, a set of learnable convolutional kernels slides along the positional arrangement direction of the one-dimensional input sequence to perform convolution operations on the positional segments within each sliding window, and a non-linear activation function is applied to extract local feature maps that reflect local chain imbalance block patterns. Pooling is performed on the local feature map output by each one-dimensional convolutional layer. The feature dimension is compressed by downsampling, which preserves salient features and expands the receptive field of subsequent convolutional layers. The feature map after the last pooling layer is flattened into a one-dimensional vector, which is used as the local nonlinear feature vector output by the second branch.

[0013] According to the technical solution provided in this application, the third branch involves performing global context capture on the genotype data and extracting global context feature vectors that characterize interactions between distant loci across the entire genome, including the following steps: The genotype data is converted into an embedding representation vector sequence, and the position encoding is superimposed on the embedding representation vector of each site to obtain an input embedding sequence with position information; The input embedding sequence with location information is input into one or more stacked self-attention layers. In each self-attention layer, the attention weight between any two sites in the input embedding sequence is calculated, and the embedding representation vectors of all sites are weighted and aggregated to capture the interaction relationship between sites at any distance across the whole genome, thus obtaining the global context representation sequence. The global context representation sequence output by the last self-attention layer is subjected to global average pooling or the representation vector of the first position of the sequence is taken to obtain the global context feature vector output by the third branch.

[0014] According to the technical solution provided in this application, the step of outputting the phenotypic prediction result through a prediction network based on the fused features includes the following steps: The fused feature vector is input into one or more stacked fully connected layers, each with a non-linear activation function. Through layer-by-layer non-linear mapping, the fused feature vector is transformed into a high-level abstract representation. The high-level abstract representation output from the last fully connected layer is input into an output fully connected layer without an activation function, and the phenotypic prediction result is obtained through linear mapping.

[0015] According to the technical solution provided in this application, the pre-training process of the prediction model includes the following steps: Obtain the genotype data and corresponding true phenotypic values ​​of the training population; The genotype data of the training population is input into the prediction model to be trained to obtain the predicted phenotypic values; A composite loss function is constructed, which is composed of a weighted sum of a mean squared error loss term and a correlation loss term based on the Pearson correlation coefficient. The mean squared error loss term is used to measure the magnitude of the error between the predicted phenotypic value and the true phenotypic value, and the correlation loss term is used to measure the degree of linear correlation between the predicted phenotypic value and the true phenotypic value. Based on the composite loss function, the loss value is calculated, and the parameters of the prediction model are updated through the backpropagation algorithm until the preset training termination condition is met, thus obtaining the pre-trained prediction model.

[0016] According to the technical solution provided in this application, the expression of the composite loss function is: L=α×L mse +β×L corr Where L is the total loss value, L mse L is the mean square error between the predicted phenotypic value and the true phenotypic value. corr The correlation loss term based on the Pearson correlation coefficient is denoted by α and β, which are preset first and second weighting coefficients, respectively. L corr The calculation method is as follows: L corr =1-|r|, where r is the Pearson correlation coefficient between the predicted phenotypic value and the true phenotypic value.

[0017] According to the technical solution provided in this application, during the pre-training process, L2 regularization constraints are applied to the model parameters of the heterogeneous feature extraction network set and / or the prediction network to suppress overfitting of the model when the number of training samples is less than the number of features.

[0018] Compared with the prior art, the beneficial effects of this application are as follows: I. Achieving explicit fusion of linear additive effects and nonlinear interaction effects. By using a heterogeneous feature extraction network set, the first feature reflecting additive genetic effects and the second feature reflecting nonlinear interaction effects are extracted from genotype data, respectively. The prior knowledge of linear additive effects in genotype data is explicitly preserved and utilized in a parallel manner. At the same time, deep learning is used to mine complex nonlinear epistatic effects. This overcomes the technical shortcomings of existing technologies that either emphasize linear or nonlinear effects, leading to an imbalance between the two, and effectively improves the model's generalization ability to new individuals.

[0019] II. Improving the fusion efficiency of multi-source features through adaptive weighted fusion. An adaptive feature fusion network is adopted to dynamically weight and fuse multiple heterogeneous features with different biological meanings using learnable fusion parameters. This enables the model to automatically learn the contribution of different features to phenotypic prediction, rather than using implicit fusion methods such as simple vector concatenation. This effectively avoids the problem of key information being buried or noise being amplified, thus improving the efficiency of feature fusion and prediction accuracy.

[0020] Third, significantly improve prediction robustness under the "large p, small n" condition. Addressing the typical problem of the number of features far exceeding the number of samples in genome prediction, the design of a heterogeneous feature extraction network set effectively controls model complexity. Combined with an adaptive fusion mechanism, it reduces the risk of overfitting, enabling the model to maintain high prediction accuracy and stability even under small sample training conditions, better meeting the needs of practical breeding applications. Attached Figure Description

[0021] Other features, objects, and advantages of this application will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings: Figure 1 The flowchart illustrates the steps of the wheat genotype-to-phenotype prediction method based on heterogeneous deep learning provided in this application. Detailed Implementation

[0022] The present application will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, only the parts relevant to the invention are shown in the accompanying drawings.

[0023] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0024] As mentioned in the background section, this application proposes a wheat genotype-to-phenotype prediction method based on heterogeneous deep learning, such as... Figure 1 As shown, it includes the following steps: S1. Obtain the genotype data of the crop to be tested, input the genotype data into the pre-trained prediction model, and obtain the phenotypic prediction results; The prediction model obtains the phenotypic prediction results in the following manner: S11. Using a heterogeneous feature extraction network set, extract multiple heterogeneous features with different biological meanings from the genotype data. The multiple heterogeneous features include at least a first feature reflecting additive genetic effects and a second feature reflecting nonlinear interaction effects. S12. Through an adaptive feature fusion network, multiple heterogeneous features are dynamically weighted and fused using learnable fusion parameters to obtain fused features; S13. Based on the fusion features, the phenotypic prediction results are output through the prediction network.

[0025] Specifically, the core idea of ​​this method is to abandon the existing approach of using only a single type of neural network to extract homogeneous features, and instead construct a predictive model that can explicitly extract and fuse features with different biological significance. The key terms and steps are explained in detail below.

[0026] In heterogeneous deep learning, "heterogeneity" refers to the fundamental differences in structure, function, and captured biological effects between the feature extraction networks within the prediction model, rather than a simple stacking of the same type of networks. Specifically, the model does not use only fully connected networks or only convolutional networks, but integrates at least two structurally different networks in parallel, each responsible for mining features of different properties from the genotype data. Genotype data refers to single nucleotide polymorphism (SNP) marker information covering the entire crop genome, obtained through high-throughput sequencing or gene chip technology. Before being input into the prediction model, this data is typically processed into a numerical matrix, where rows represent different crop individuals, columns represent different SNP loci, and the element values ​​in the matrix are the genotype codes for that locus. For example, using a 0, 1, 2 encoding format, these represent homozygous reference, heterozygous, and homozygous substitute alleles, respectively. Phenotypic prediction results refer to the quantitative prediction values ​​of target agronomic traits (such as yield, plant height, thousand-grain weight), which are continuous numerical values.

[0027] One of the key innovations of this application is the use of a heterogeneous feature extraction network set to extract multiple heterogeneous features with different biological meanings from the genotype data. Here, the heterogeneous feature extraction network set refers to a collection of at least two structurally different neural network branches, which are connected in parallel and independently process the input genotype data. Multiple heterogeneous features with different biological meanings mean that the feature vectors extracted by these branches respectively characterize biological effects at different levels in the genotype data. These multiple heterogeneous features include at least a first feature reflecting additive genetic effects and a second feature reflecting nonlinear interaction effects. Additive genetic effects refer to the independent contribution of alleles at a single locus to the phenotype; it is the classic basis of numerical trait inheritance, manifested as a linear relationship between genotype and phenotype. Therefore, the first feature aims to capture this linear additive effect. Nonlinear interaction effects encompass epistatic effects between genes (i.e., nonlinear interactions between different loci) and complex interactions between genes and the environment. The second feature is precisely designed to capture these complex relationships that cannot be effectively modeled by traditional linear models. Through this parallel extraction method, this application is the first to explicitly combine linear additive prior knowledge with nonlinear deep learning capabilities at the model architecture level, overcoming the information loss or overfitting problems caused by existing technologies that either emphasize linearity or nonlinearity.

[0028] Next, the process utilizes an adaptive feature fusion network to dynamically weight and fuse multiple heterogeneous features using learnable fusion parameters, thus solving the challenge of efficiently integrating multi-source features. Existing techniques typically employ simple vector concatenation, an implicit fusion method that fails to distinguish the differences in contribution of different features to the final prediction task. The introduced adaptive feature fusion network is a trainable module containing learnable fusion parameters that are automatically updated through backpropagation during model training. Its working principle is as follows: each heterogeneous feature is assigned a weight coefficient, the magnitude of which dynamically reflects the importance of the feature in predicting a specific phenotype. Then, all features are weighted and summed according to these weights, rather than being treated as equal weights. The dynamic aspect here is reflected in the fact that the final weight values ​​are adaptively determined based on the input data and the prediction target, rather than being fixed values ​​set manually. Through this mechanism, if a linear effect is more important for a certain trait, the weight of the first feature will be automatically increased; conversely, if the nonlinear interaction has a greater impact, the weight of the second feature will be strengthened. This effectively avoids information redundancy and conflict, enabling the fused features to retain the most valuable information for the prediction task to the greatest extent.

[0029] Finally, based on the fused features, the phenotypic prediction result is output through a prediction network. The prediction network typically consists of several fully connected layers, responsible for mapping the high-dimensional fused features to the final scalar prediction value.

[0030] The technical effect of this implementation method is that by combining heterogeneous feature extraction and adaptive fusion, the model can simultaneously capture the linear additive effect and nonlinear interaction effect of genotype data, and automatically learn the best combination of different effects. Thus, in the typical scenario of genome prediction with "large p and small n" (the number of features is much larger than the number of samples), the prediction accuracy and model generalization ability of new individuals are significantly improved, and the defects of overfitting of a single type of model are avoided.

[0031] In a preferred embodiment, the extraction of multiple heterogeneous features from the genotype data using a heterogeneous feature extraction network set includes the following steps: The genotype data is simultaneously input into the first, second, and third parallel branches, respectively. In the first branch, the genotype data is linearly mapped, and a linear feature vector representing the additive genetic effect is extracted as the first feature; In the second branch, local pattern capture is performed on the genotype data to extract local nonlinear feature vectors that characterize local linkage disequilibrium block patterns, which are used as a component of the second feature. In the third branch, global context capture is performed on the genotype data to extract global context feature vectors that characterize the interaction relationships between distant sites across the entire genome, as another component of the second feature; The linear feature vector output from the first branch, the local nonlinear feature vector output from the second branch, and the global context feature vector output from the third branch are input into an adaptive feature fusion network, and the three are fused into a unified fusion feature vector through learnable fusion parameters.

[0032] Specifically, firstly, the step of simultaneously inputting the genotype data into the first, second, and third parallel branches clarifies the direction of the data flow. This means that the genotype matrix input to the model is copied three times, and the same batch of data is simultaneously fed into three structurally distinct branches for parallel processing, ensuring the efficiency and independence of feature extraction. The three branches are connected in parallel and each independently extracts features from the genotype data, meaning that the forward propagation process of any branch does not depend on the output of other branches; their functions are decoupled.

[0033] The first branch is defined as performing a linear mapping on the genotype data to extract linear feature vectors representing additive genetic effects. Here, linear mapping refers to a purely linear transformation of the input data through a network structure composed of stacked fully connected layers without nonlinear activation functions. Its biological significance lies in the fact that additive effect models in quantitative genetics are inherently linear; therefore, a linear branch that does not introduce nonlinear distortion can preserve and amplify this effect intact, avoiding the problem of deep learning models forgetting linear priors due to forced nonlinear transformations.

[0034] The second branch is defined as capturing local patterns in the genotype data, extracting local nonlinear feature vectors characterizing local linkage disequilibrium block patterns. In genomics, due to varying recombination rates, adjacent SNP loci on chromosomes often form linkage disequilibrium blocks, meaning that the alleles at these loci are inherited together non-randomly. The SNP combination patterns within these blocks have a significant impact on phenotype. The second branch is specifically responsible for capturing these short-range local sequence patterns. Its local pattern capture is typically implemented using a one-dimensional convolutional neural network, where the convolution kernel slides along the SNP sequence direction, performing convolution operations on loci within each window to detect biologically significant local motifs. Because a nonlinear activation function is introduced, the features extracted by this branch are called local nonlinear feature vectors.

[0035] The third branch is defined as capturing the global context of the genotype data, extracting global contextual feature vectors that characterize long-distance interactions between loci across the entire genome. Besides local block effects, phenotypes are also influenced by long-range regulatory relationships between genes scattered across different locations in the genome, and even on different chromosomes. These long-distance dependencies are difficult to capture due to the limited receptive field of convolutional neural networks. The third branch addresses this problem by introducing mechanisms with global receptive fields (such as self-attention networks). It calculates the association strength between any two SNP sites in the sequence, thereby constructing a genome-wide interaction map and extracting global contextual feature vectors.

[0036] Finally, the linear feature vector output from the first branch, the local nonlinear feature vector output from the second branch, and the global context feature vector output from the third branch are input into the adaptive feature fusion network. Through learnable fusion parameters, the three are fused into a unified fusion feature vector. This explicitly indicates that the specific feature objects to be fused are three types of heterogeneous feature vectors, representing a concretization of multiple heterogeneous features.

[0037] The technical advantage of this implementation lies in constructing a multi-scale, comprehensive feature extraction system that ranges from local to global and from linear to nonlinear. The three branches correspond to three different dimensions of biological mechanisms: additive effects, local LD ​​block effects, and long-distance epistatic effects, enabling the model to more comprehensively and deeply mine genotype information. This explicit decomposition based on biological priors, compared to blind end-to-end black-box models, not only improves predictive performance but also enhances the model's interpretability.

[0038] In a preferred embodiment, the dynamic weighted fusion of the multiple heterogeneous features using an adaptive feature fusion network includes the following steps: Assign a first learnable weight coefficient to the linear feature vector output by the first branch, assign a second learnable weight coefficient to the local nonlinear feature vector output by the second branch, and assign a third learnable weight coefficient to the global context feature vector output by the third branch. Multiply the linear feature vector by the first learnable weight coefficient to obtain a weighted linear feature vector; Multiply the local nonlinear feature vector by the second learnable weight coefficient to obtain a weighted local nonlinear feature vector; Multiply the global context feature vector by the third learnable weight coefficient to obtain the weighted global context feature vector; The fused feature vector is obtained by summing the weighted linear feature vector, the weighted local nonlinear feature vector, and the weighted global context feature vector element by element.

[0039] Specifically, the learnable weight coefficients are determined by constructing a lightweight weight generation network within the adaptive feature fusion network, rather than simply setting three independent trainable scalars. The specific weight determination and dynamic weighted fusion process is as follows.

[0040] First, the linear feature vector output from the first branch, the local nonlinear feature vector output from the second branch, and the global context feature vector output from the third branch are denoted as F1, F2, and F3, respectively, all with dimension d. The first step in the weight generation network is to concatenate these three feature vectors along their respective feature dimensions to obtain a joint feature vector Fconcat with a dimension of 3d. This concatenation operation preserves complete information from all branches, providing a basis for subsequently determining the relative importance of each branch.

[0041] Next, the joint feature vector Fconcat is input into a weight generation subnetwork. This subnetwork is structured as follows: First, Fconcat is linearly transformed through a fully connected layer, compressing its dimension from 3d to a smaller intermediate dimension h, i.e., dimensionality reduction mapping, to extract the compact information most relevant to the weight decisions. Then, a non-linear activation function (such as ReLU) is applied to this dimensionality-reduced vector, introducing non-linearity to enhance the expressive power of the weight generation network. Immediately afterwards, this vector is passed through a second fully connected layer, mapping its dimension from h to 3, resulting in an original score vector containing three values. These three values ​​correspond to the original importance scores of the first, second, and third branches, respectively.

[0042] To ensure these three scores become reasonable weight coefficients, normalization is required. Specifically, the original score vector is input into the Softmax activation function, the exponent value corresponding to each score is calculated, and each exponent value is divided by the sum of the three exponent values ​​to obtain three weight coefficients with a sum of 1. These three weight coefficients are the first learnable weight coefficient w1, the second learnable weight coefficient w2, and the third learnable weight coefficient w3. Mathematically, this is expressed as: [w1, w2, w3] = Softmax(W2 · ReLU(W1 · Fconcat + b1) + b2), where W1 and W2 are the weight matrices of the weight generation sub-network, and b1 and b2 are bias terms. These parameters are all learnable fusion parameters, which are optimized and updated along with other parameters of the prediction model during model training via backpropagation. The function of the entire weight generation sub-network is to dynamically determine which branch's features are more critical to the prediction task of the current sample based on the specific content of the three heterogeneous features of the current input sample, and assign it a higher weight coefficient.

[0043] After determining the weight coefficients, a weighted fusion operation is performed. The linear feature vector F1 is multiplied by the first learnable weight coefficient w1 to obtain a weighted linear feature vector; the local nonlinear feature vector F2 is multiplied by the second learnable weight coefficient w2 to obtain a weighted local nonlinear feature vector; and the global context feature vector F3 is multiplied by the third learnable weight coefficient w3 to obtain a weighted global context feature vector. These three weighting operations scale the magnitude of each feature vector as a whole. Finally, the three weighted feature vectors are summed element-wise, that is, the elements at corresponding positions are added together to obtain a unified fused feature vector with the same dimension d.

[0044] The technical advantage of this implementation lies in the following: by introducing a dedicated weight generation subnetwork, the process of determining weight coefficients is elevated from a simple parameter lookup to a dynamic decision-making process based on the input data content. This subnetwork can automatically evaluate the contribution of linear, locally nonlinear, and global context features to the phenotypic prediction of the current input genotype data according to its specific characteristics, and dynamically allocate weights accordingly. This data-dependent dynamic fusion strategy is more adaptable to the differences in genetic mechanisms among different samples than fixed weights or simple vector concatenation, thereby significantly improving the information density of the fused features and the overall accuracy of the prediction model. Simultaneously, through Softmax normalization, the sum of the weight coefficients is kept constant at 1, which provides a good mathematical constraint for model training and ensures training stability.

[0045] In a preferred embodiment, the step of linearly mapping the genotype data in the first branch to extract linear feature vectors representing additive genetic effects includes the following steps: The input genotype data is standardized to obtain a standardized genotype matrix; The standardized genotype matrix is ​​input into one or more stacked fully connected layers, wherein each fully connected layer does not have a non-linear activation function, and an initial linear feature vector is obtained through layer-by-layer linear transformation; The initial linear feature vector output by each fully connected layer is subjected to batch normalization and layer scaling in sequence. The batch normalization process is used to adjust the distribution of the initial linear feature vector to zero mean and unit variance. The layer scaling process is used to multiply the normalized feature vector with a learnable scaling factor to adjust the contribution of the layer output to the final linear feature vector. The feature vectors after batch normalization and layer scaling are used as the linear feature vectors output by the first branch.

[0046] First, the input genotype data is standardized to obtain a standardized genotype matrix. The original input genotype data is an SNP encoding matrix with values ​​of 0, 1, and 2. The standardization process involves calculating the mean and standard deviation of each SNP locus (i.e., each column of the matrix) across the entire training population. Then, the mean is subtracted from each value in that column, and the result is divided by the standard deviation, yielding a standardized value with a mean of 0 and a standard deviation of 1. This step eliminates the numerical distribution differences between different SNP loci caused by varying allele frequencies, ensuring that the input features of all loci are on the same numerical scale, providing a numerically stable starting point for subsequent linear transformations.

[0047] The standardized genotype matrix is ​​then input into one or more stacked fully connected layers. For each fully connected layer, the mathematical operation is Y = XW + b, where X is the input matrix, W is the weight matrix, and b is the bias vector. The key constraint is that no non-linear activation function is applied after each fully connected layer; the output Y is directly used as the input to the next layer or as the final output of the current layer. When multiple such fully connected layers without activation functions are stacked sequentially, their composite transformation is mathematically equivalent to the input data passing through a single fully connected layer, because the combination of multiple linear transformations is still a linear transformation. The aim of this design is to always keep the feature space within the original linear transformation range of the genotype data, ensuring that the extracted feature vectors are linear combinations of the original SNP information, thereby explicitly capturing additive genetic effects—the independent and additive contribution of each allele to the phenotype.

[0048] To ensure training stability and prevent overfitting in purely linear networks, batch normalization and layer scaling are performed sequentially after the linear transformation of each fully connected layer. Batch normalization involves calculating the mean and variance of each feature dimension of the layer's output in the current mini-batch, normalizing it to a distribution with a mean of 0 and a variance of 1, and then scaling and shifting. This effectively addresses the internal covariate shift problem caused by increasing network depth. Layer scaling introduces a learnable scaling factor, typically set to a small initial value (e.g., 0.1). This scaling factor is multiplied by the batch-normalized feature vector, adjusting the overall magnitude of the layer's output. In the early stages of training, a smaller scaling factor limits the contribution of linear branches to the fusion layer, causing the model to preferentially rely on other more easily trained nonlinear branches. Once training stabilizes, this factor can adaptively increase, thereby adjusting the contribution of the layer's output to the final linear feature vector.

[0049] Finally, the final feature vector after the above multi-layer linear transformation, batch normalization and layer scaling is used as the linear feature vector output by the first branch.

[0050] The technical advantage of this implementation lies in the following: by stacking fully connected layers without activation functions, and supplementing this with batch normalization and layer scaling, a stable and controllable purely linear feature extraction channel is successfully constructed within a deep learning framework. This not only solidifies the prior knowledge that additive effects are linear relationships, avoiding spurious correlations that may be introduced by nonlinear transformations, but also provides a progressive training strategy from weak to strong through layer scaling, effectively alleviating the training conflict between linear and nonlinear branches. This allows the entire heterogeneous model to more fully utilize linear additive information, improving generalization prediction capabilities.

[0051] In a preferred embodiment, the step of performing local pattern capture on the genotype data and extracting local nonlinear feature vectors characterizing local linkage disequilibrium block patterns in the second branch includes the following steps: The genotype data is adjusted to a one-dimensional sequence format of single nucleotide polymorphism sites to obtain a one-dimensional input sequence; The one-dimensional input sequence is input into one or more stacked one-dimensional convolutional layers. In each one-dimensional convolutional layer, a set of learnable convolutional kernels slides along the positional arrangement direction of the one-dimensional input sequence to perform convolution operations on the positional segments within each sliding window, and a non-linear activation function is applied to extract local feature maps that reflect local chain imbalance block patterns. Pooling is performed on the local feature map output by each one-dimensional convolutional layer. The feature dimension is compressed by downsampling, which preserves salient features and expands the receptive field of subsequent convolutional layers. The feature map after the last pooling layer is flattened into a one-dimensional vector, which is used as the local nonlinear feature vector output by the second branch.

[0052] First, the genotype data is converted into a one-dimensional sequence format for single nucleotide polymorphism (SNP) sites, resulting in a one-dimensional input sequence. The original input genotype data is in matrix format, where each row represents a crop individual and each column represents a SNP site, arranged in order of their physical location on the chromosome or their position on the genetic map. Extracting each row of this matrix forms a one-dimensional sequence with a length equal to the total number of SNP sites. The value at each position in the sequence is the genotype code for that site: 0, 1, or 2. This conversion provides a standard one-dimensional sequence input for subsequent convolution operations, preserving the order information between sites.

[0053] Next, the one-dimensional input sequence is fed into one or more stacked one-dimensional convolutional layers. Each convolutional layer contains a set of learnable convolutional kernels, each a weight vector of length k, where k is the kernel size. During convolution, each kernel slides sequentially backward along the positional alignment of the one-dimensional input sequence, starting from the first position, with a stride of one or more positions. At each sliding window position, the kernel weights are element-wise multiplied with the corresponding input value within the window, and then summed to obtain an output value. After the convolutional kernel has traversed the entire input sequence, a one-dimensional feature map is generated. By applying a non-linear activation function (such as ReLU) after the convolution operation, this layer can capture non-linear patterns of genotype combinations within the locus fragment. Multiple convolutional kernels can learn a variety of different local motifs, and stacking multiple convolutional layers can extract hierarchical local patterns from simple to complex.

[0054] In biology, SNPs that are physically adjacent on chromosomes often exhibit a high degree of linkage disequilibrium, meaning that certain allele combinations tend to be inherited together, forming linkage disequilibrium blocks. The sliding of the convolution kernel along the sequence can precisely capture these local window patterns of site combinations, thus the extracted feature maps can reflect the patterns of local linkage disequilibrium blocks.

[0055] Then, pooling is performed on the local feature maps output by each one-dimensional convolutional layer. Common pooling operations include max pooling and average pooling. Max pooling takes the maximum value within the pooling window as the output, preserving the most salient feature signal; average pooling takes the average value within the window. Pooling compresses the feature dimension and reduces computation by downsampling the output feature map to a lower resolution than the input feature map, while simultaneously expanding the effective receptive field of subsequent convolutional layers, enabling them to indirectly perceive a wider range of sequence context.

[0056] Finally, the feature map obtained after the last convolution and pooling layer is flattened into a one-dimensional vector. The flattening operation refers to concatenating all elements in the feature map sequentially into a long vector. This vector is the local nonlinear feature vector output by the second branch, which integrates multi-level, multi-scale local LD ​​block information extracted from the genotype sequence.

[0057] The technical advantage of this implementation lies in the following: By employing multi-level one-dimensional convolution and pooling operations, a hierarchical feature extraction channel is constructed, ranging from local motifs to larger-scale block patterns. This design precisely matches the biological characteristic of linked imbalance blocks being continuously distributed and hierarchically structured on chromosomes, enabling the model to automatically discover SNP combination patterns significantly related to the target phenotype without the need for manually pre-defining haplotype blocks. Simultaneously, the translation invariance introduced by pooling operations makes the model robust to small shifts in block boundaries, further enhancing the stability of feature extraction and the generalization ability of prediction.

[0058] In a preferred embodiment, the step of performing global context capture on the genotype data and extracting global context feature vectors characterizing interactions between distant loci across the entire genome in the third branch includes the following steps: The genotype data is converted into an embedding representation vector sequence, and the position encoding is superimposed on the embedding representation vector of each site to obtain an input embedding sequence with position information; The input embedding sequence with location information is input into one or more stacked self-attention layers. In each self-attention layer, the attention weight between any two sites in the input embedding sequence is calculated, and the embedding representation vectors of all sites are weighted and aggregated to capture the interaction relationship between sites at any distance across the whole genome, thus obtaining the global context representation sequence. The global context representation sequence output by the last self-attention layer is subjected to global average pooling or the representation vector of the first position of the sequence is taken to obtain the global context feature vector output by the third branch.

[0059] First, the genotype data is converted into an embedding representation vector sequence, and positional encoding is superimposed on the embedding representation vector of each locus to obtain an input embedding sequence with positional information. The input genotype data is an SNP encoding matrix, where the original value of each locus is a discrete integer encoding such as 0, 1, or 2. This low-dimensional discrete representation is difficult to directly express the complex semantic relationships between loci, so embedding transformation is required. The specific operation of embedding transformation is to assign a learnable d-dimensional embedding vector to each possible genotype value (0, 1, 2). This embedding vector is randomly assigned during model initialization and continuously optimized through backpropagation during training. After replacing the discrete encoding of each locus in the input sequence with its corresponding embedding vector, the entire sequence is converted into an embedding representation vector sequence. In addition, the self-attention mechanism itself is not aware of the positional order of elements in the sequence during computation; that is, it cannot distinguish the computational difference between the embedding vectors of the first locus and the tenth locus. In order for the model to utilize the physical order information of SNP loci on the chromosome, positional information needs to be explicitly injected into the embedding vector. The specific implementation of positional encoding involves generating a unique d-dimensional positional encoding vector for each location. This vector can be calculated by combining sine and cosine functions, or it can be set as a learnable parameter vector. Then, the embedding representation vector of each location is directly added to its corresponding element-wise. The result of superimposing the positional information is the input embedding sequence with positional information.

[0060] Next, the input embedding sequence with location information is fed into one or more stacked self-attention layers. The core computational process of a self-attention layer consists of three steps. The first step is to generate query vectors, key vectors, and value vectors. The d-dimensional vector of each site in the input embedding sequence is multiplied by three learnable weight matrices Wq, Wk, and Wv, respectively, to obtain the corresponding query vector q, key vector k, and value vector v, which are typically d in dimension. The parameters of these three matrices are central to the self-attention layer's learning, determining how the model evaluates the relationships between sites and how it aggregates information. The second step is to calculate attention weights, i.e., to evaluate the interaction strength between any two sites. For each site i, its query vector qi is used to perform a dot product with the key vectors kj of all sites j in the sequence, resulting in a raw attention score that reflects the importance of site j to site i. To stabilize training, the dot product result is usually divided by a scaling factor (often the square root of dk, where dk is the dimension of the key vectors) to prevent the dot product value from becoming too large and causing the softmax gradient to vanish. Then, a Softmax function is applied to the scaled scores of all loci j, normalizing them to a probability distribution such that the sum of the attention weights of all loci to loci i is 1. The third step is weighted aggregation. The calculated attention weights are used to sum the value vectors vj of all loci in the sequence to obtain the output vector of loci i. This operation makes the output vector of each locus no longer its isolated representation, but a context-aware representation that integrates information from all other loci across the entire genome. Applying the above process to each locus in the sequence yields the global context representation sequence. In practical applications, self-attention layers typically employ a multi-head attention mechanism, that is, executing multiple independent sets of the above computational processes in parallel. Each set is called an attention head, and each set has its own independent Wq, Wk, and Wv matrices. Different attention heads can focus on different types of inter-locus relationships. Finally, the outputs of each head are concatenated and mapped back to the original dimension through a fully connected layer, thereby comprehensively capturing the complex interaction relationships between distant loci from multiple perspectives. The sequential stacking of multiple self-attention layers enables the model to learn multi-level global interaction patterns from low to high order.

[0061] Finally, after obtaining the global context representation sequence output by the last self-attention layer, this sequence needs to be compressed into a fixed-dimensional global feature vector. There are two ways to achieve this. One is global average pooling, which averages the output vectors of all positions in the sequence across each feature dimension to obtain a d-dimensional vector. The other is to use the representation vector of the first position in the sequence, i.e., using only the output vector corresponding to the first position as the global feature vector. Both methods can effectively summarize the information of the entire sequence to obtain the global context feature vector output by the third branch.

[0062] The technical advantages of this implementation are as follows: The third branch, through a self-attention mechanism, completely breaks the physical limitations of the receptive field of convolutional neural networks, endowing the model with the ability to directly model the interaction relationship between any two SNP sites within a single layer. This allows the interaction information of long-distance regulatory elements scattered across different locations in the genome, and even different chromosomes, to be effectively captured, compensating for the second branch's limitation of only perceiving local LD ​​block patterns. Simultaneously, the multi-head attention mechanism enables the model to learn different types of long-distance relationships in parallel from multiple representation subspaces, enhancing the richness and robustness of feature extraction. By combining the local information of the second branch with the global contextual information of the third branch, the model can construct a complete genome feature map from local to global perspectives, more comprehensively depicting the complex mapping from genotype to phenotype.

[0063] In a preferred embodiment, the step of outputting the phenotypic prediction result through a prediction network based on the fused features includes the following steps: The fused feature vector is input into one or more stacked fully connected layers, each with a non-linear activation function. Through layer-by-layer non-linear mapping, the fused feature vector is transformed into a high-level abstract representation. The high-level abstract representation output from the last fully connected layer is input into an output fully connected layer without an activation function, and the phenotypic prediction result is obtained through linear mapping.

[0064] Specifically, firstly, the fused feature vector is input into one or more stacked fully connected layers. Each fully connected layer is followed by a non-linear activation function. Through layer-by-layer non-linear mapping, the fused feature vector is transformed into a high-level abstract representation. The fused feature vector here is a unified feature vector output by the adaptive feature fusion network, typically with dimension d. It integrates linear additive information from the first branch, local LD ​​block pattern information from the second branch, and genome-wide long-range interaction information from the third branch. This fused feature vector is first fed into the first fully connected layer. The mathematical operation of a fully connected layer is Y = XW + b, where X is the input fused feature vector, W is the weight matrix, and b is the bias vector. The output dimension of this layer is typically smaller than the input dimension d, i.e., dimensionality reduction mapping is performed to progressively compress information and extract more abstract features. After the linear transformation of the fully connected layer, a non-linear activation function follows. Commonly used non-linear activation functions include ReLU, with the functional form f(x) = max(0, x), which sets the input values ​​less than zero to zero and leaves the values ​​greater than zero unchanged. The introduction of nonlinear activation functions endows predictive networks with nonlinear fitting capabilities, enabling them to learn complex nonlinear interactions between feature dimensions. Without nonlinear activation functions, a predictive network, regardless of the number of layers, is mathematically equivalent to a single linear regression model, unable to fit the complex nonlinear relationships that may exist in the genotype-phenotype mapping. When a predictive network contains multiple fully connected layers, these layers are sequentially connected, with the output of the previous layer serving as the input to the next, and the output dimension decreasing layer by layer, forming a funnel-shaped structure from high to low dimension. This gradual dimensionality reduction design allows the network to extract the most relevant abstract features to the target phenotype layer by layer, filtering out redundant information.

[0065] Then, the high-level abstract representation output from the last fully connected layer is input into an output fully connected layer without an activation function, and the phenotypic prediction result is obtained through linear mapping. After feature abstraction and dimensionality reduction in multiple nonlinear fully connected layers, the last fully connected layer is designed as a purely linear layer without an activation function. The input to this layer is the high-level abstract representation vector output from the penultimate layer, which has undergone multiple nonlinear transformations and its dimension has been reduced to a small value, such as 64 or 32. The weight matrix of the output fully connected layer linearly maps this vector to a scalar value, which is the final phenotypic prediction result. The reason for not using an activation function in the output layer is that the phenotypic prediction task is a regression task, and its prediction target is a continuous numerical value (such as yield, plant height, thousand-grain weight, etc.), and the output value range should be the entire real number range. If an activation function such as ReLU is applied to the output layer, it will force the output value to be non-negative, which conflicts with the requirement that the actual phenotypic value may be negative (after standardization) or that it needs to cover the positive and negative ranges. Applying activation functions such as Sigmoid or Tanh will restrict the output value to a specific interval, which is also unsuitable for regression prediction of unbounded continuous values. Therefore, a linear output layer is the correct choice for regression prediction tasks.

[0066] The technical advantages of this implementation are as follows: The prediction network, through a combination of multiple nonlinear fully connected layers and linear output layers, constructs a regression prediction module that possesses both strong nonlinear fitting capabilities and the ability to output continuous, unbounded prediction values. The funnel-shaped structure with layer-by-layer dimensionality reduction effectively compresses the feature dimension, reduces computational overhead, and simultaneously acts as an implicit feature selection mechanism, enabling the network to automatically focus on the most critical abstract features for the prediction task. Maintaining linearity in the output layer ensures that the predicted values ​​are mathematically within the same numerical space as the true phenotypic values, avoiding systematic biases introduced by mismatched activation functions in the output layer. This design enables the entire heterogeneous deep learning prediction model to complete the end-to-end mapping from raw SNP genotype data to continuous phenotypic prediction values.

[0067] In a preferred embodiment, the pre-training process of the prediction model includes the following steps: Obtain the genotype data and corresponding true phenotypic values ​​of the training population; The genotype data of the training population is input into the prediction model to be trained to obtain the predicted phenotypic values; A composite loss function is constructed, which is composed of a weighted sum of a mean squared error loss term and a correlation loss term based on the Pearson correlation coefficient. The mean squared error loss term is used to measure the magnitude of the error between the predicted phenotypic value and the true phenotypic value, and the correlation loss term is used to measure the degree of linear correlation between the predicted phenotypic value and the true phenotypic value. Based on the composite loss function, the loss value is calculated, and the parameters of the prediction model are updated through the backpropagation algorithm until the preset training termination condition is met, thus obtaining the pre-trained prediction model.

[0068] First, obtain the genotype data and corresponding true phenotypic values ​​of the training population. The training population is the foundation of genome selection; it consists of a group of crop individuals whose genotyping and precise phenotypic measurements have been completed. Genotype data is the whole-genome SNP coding matrix for each individual in the population, with values ​​of 0, 1, or 2, representing homozygous reference, heterozygous, and homozygous substitute alleles, respectively. True phenotypic values ​​refer to the measured values ​​of the target agronomic traits to be predicted, obtained from field trials or laboratory measurements of each individual in the population, such as yield data corrected through multi-year, multi-location trials, precise plant height, and thousand-grain weight. These true phenotypic values ​​typically require preprocessing before being input into the training process, including removing outliers and standardizing the phenotypic values ​​of all individuals to ensure a mean of zero and a variance of one, thus eliminating dimensional differences between different traits and making loss calculations more numerically stable. The sample size of the training population directly affects the model's generalization ability. In typical genome prediction scenarios, the number of training samples is usually much smaller than the number of SNP loci, exhibiting a "large p, small n" characteristic.

[0069] Next, the genotype data of the training population is input into the prediction model to be trained, yielding predicted phenotypic values. At the start of training, all learnable parameters of the prediction model, including the weight matrices and bias vectors of each branch in the heterogeneous feature extraction network set, the fusion parameters in the adaptive feature fusion network, and the weight matrix and bias vector of the prediction network, are randomly initialized. The SNP encoding matrices of the training population are input into the model batch by batch. After linear mapping in the first branch, local pattern capture in the second branch, global context capture in the third branch, dynamic weighted fusion in the adaptive feature fusion network, and regression mapping in the prediction network, the predicted phenotypic value for each individual is finally output. This complete forward propagation process transforms the raw genotype data into the model's estimate of the target phenotypic value.

[0070] Then, a composite loss function is constructed, which consists of a weighted sum of a mean squared error loss term and a correlation loss term based on the Pearson correlation coefficient. The mean squared error loss term measures the magnitude of the error between the predicted phenotypic value and the true phenotypic value. It is calculated as follows: for all individuals in each training batch, the square of the difference between the predicted value and the true value is calculated, and then the average of the squared differences of all individuals is taken. The smaller the mean squared error, the smaller the absolute error of the model's prediction, which is the most basic optimization objective in regression tasks. However, in breeding practice, the focus is often not on the prediction error of a single individual, but on the consistency between the ranking of predicted values ​​and the ranking of true values, i.e., the accuracy of the prediction. To directly optimize this indicator, a correlation loss term based on the Pearson correlation coefficient is additionally introduced into the loss function. The Pearson correlation coefficient r ranges from -1 to 1; the closer its absolute value is to 1, the higher the linear correlation between the predicted and true values. The correlation loss term is calculated as Lcorr = 1 - |r|, meaning the closer the absolute value of the correlation coefficient is to 1, the closer this loss term is to zero. The mean squared error loss term and the correlation loss term are weighted and summed using a first weighting coefficient α and a second weighting coefficient β to form a composite loss function L = α × Lmse + β × Lcorr. The two weighting coefficients are pre-defined hyperparameters used to adjust the model's emphasis between minimizing absolute error and maximizing ranking correlation.

[0071] Finally, based on the composite loss function, the loss value is calculated, and the parameters of the prediction model are updated using the backpropagation algorithm until the preset training termination condition is met, resulting in a pre-trained prediction model. After calculating the loss value L in each batch of forward propagation, the backpropagation algorithm calculates the partial derivative of the loss value with respect to each learnable parameter, i.e., the gradient, layer by layer, starting from the output layer and moving in the reverse direction of the network. Then, the gradient descent optimizer is used to update each parameter according to the calculated gradient and the preset learning rate, so that the loss value is reduced in the next forward propagation. Optimizers such as Adam and SGD can be used. One complete iteration of all batches of training data is called a training epoch. The model needs to undergo multiple training epochs of iteration, monitoring the loss value and prediction accuracy on the validation set in each epoch. When the loss value on the validation set no longer decreases, or when the preset maximum number of training epochs is reached, training terminates, the model parameters at this point are saved, and the pre-trained prediction model is obtained.

[0072] The technical advantages of the above solution are as follows: By employing a composite loss function consisting of a mean squared error loss term and a Pearson correlation coefficient loss term, the pre-training process can simultaneously optimize two key indicators: the magnitude of the model's prediction error and the accuracy of its prediction ranking. The mean squared error loss term ensures that the model's predicted values ​​are numerically close to the true values, avoiding serious biases; while the correlation coefficient loss term guides the model to focus on the consistency of predicted rankings among individuals, directly optimizing the accuracy requirements in breeding selection. This dual-objective optimization strategy, compared to using only a single mean squared error loss function, more comprehensively meets the practical needs of genomic selection applications, significantly enhancing the model's practical value in breeding material screening. Simultaneously, by jointly optimizing all parameters of the heterogeneous feature extraction network set, the adaptive feature fusion network, and the prediction network through the backpropagation algorithm, it ensures that all components of the entire model can work collaboratively, enabling the linear features of the first branch, the local nonlinear features of the second branch, the global context features of the third branch, and the adaptive fusion weights to achieve optimal configuration under a unified optimization objective.

[0073] In a preferred embodiment, the expression for the composite loss function is: L=α×L mse +β×L corr Where L is the total loss value, L mse L is the mean square error between the predicted phenotypic value and the true phenotypic value. corr The correlation loss term based on the Pearson correlation coefficient is denoted by α and β, which are preset first and second weighting coefficients, respectively. L corr The calculation method is as follows: L corr =1-|r|, where r is the Pearson correlation coefficient between the predicted phenotypic value and the actual phenotypic value.

[0074] Specifically, the composite loss function consists of the product of the first weighting coefficient and the mean squared error loss term, plus the product of the second weighting coefficient and the correlation loss term. The first and second weighting coefficients are pre-set hyperparameters used to adjust the model's emphasis between the two optimization objectives. When the breeding task focuses more on the accurate estimation of individual phenotypic values, a larger first weighting coefficient can be set; when the task focuses more on the accurate ranking and selection of superior individuals in the candidate population, a larger second weighting coefficient can be set. The two coefficients typically satisfy normalization constraints, but can also be set independently according to actual needs.

[0075] The mean squared error (MSE) loss term measures the magnitude of the absolute error between the predicted phenotypic value and the true phenotypic value. Its calculation steps are as follows: for each crop individual in the current training batch, calculate the difference between its predicted and true phenotypic values, square this difference, then sum the squared differences of all individuals in the batch, and divide by the total number of individuals in the batch. The result is the MSE loss value. The MSE loss value is always non-negative, reaching its theoretical minimum of zero only when all predicted values ​​are exactly equal to the true values. Because the error is squared, predicted values ​​that deviate further from the true value are penalized more, causing the model to prioritize correcting predictions with particularly large deviations during training. The MSE loss term ensures that the predicted values ​​output by the model remain within a similar scale to the true phenotypic values ​​in absolute terms.

[0076] The correlation loss term based on the Pearson correlation coefficient measures the degree of linear correlation between predicted and actual phenotypic values. The calculation steps are as follows: First, calculate the mean of predicted phenotypic values ​​for all individuals in the current training batch, and the mean of actual phenotypic values ​​for all individuals. Then, for each individual, calculate the difference between their predicted value and the predicted mean, and the difference between their actual value and the actual mean. Multiply these two differences for each individual, and sum the products of all individuals in the batch to obtain the numerator. Simultaneously, calculate the sum of squares of the differences between predicted values ​​and the predicted mean for all individuals, and the sum of squares of the differences between actual values ​​and the actual mean for all individuals. Take the square root of each of these sums and multiply them to obtain the denominator. Divide the numerator by the denominator; the quotient is the Pearson correlation coefficient, which ranges from -1 to +1. The closer the Pearson correlation coefficient is to positive one, the stronger the positive linear correlation between the predicted and actual values. That is, individuals with high predicted values ​​also have high actual values, and individuals with low predicted values ​​also have low actual values, and the better the consistency of the ranking among individuals. The closer it is to zero, the weaker the linear relationship between the two. The closer it is to negative one, the more negative the linear correlation between the two.

[0077] After obtaining the Pearson correlation coefficient, its absolute value is calculated, and then subtracted from one to obtain the correlation loss value. The logic behind this design is that the closer the absolute value of the correlation coefficient is to one, the closer the correlation loss value is to zero, indicating a greater consistency between the model's predicted ranking and the actual ranking; conversely, the closer the absolute value of the correlation coefficient is to zero, the closer the correlation loss value is to one, indicating a weaker predictive ranking ability of the model. By adding this term to the total loss function and minimizing it, the model is guided to maximize the consistency between the predicted ranking and the actual ranking. It should be noted that the reason for taking the absolute value of the correlation coefficient is that in breeding practice, the focus is on the accuracy of prediction, i.e., the strength of the correlation itself, without distinguishing between positive and negative correlation directions. As long as the correlation is strong, it indicates that the model has effectively captured the inherent mapping law from genotype to phenotype.

[0078] Finally, the product of the first weighting coefficient and the mean squared error loss value is added to the product of the second weighting coefficient and the correlation loss value to obtain the total loss value for this training batch. During the backpropagation phase of model training, the optimizer aims to minimize this total loss value, calculates the gradient of the total loss value with respect to all learnable parameters of the model, and updates the parameters.

[0079] In a preferred embodiment, during the pre-training process, L2 regularization constraints are applied to the model parameters of the heterogeneous feature extraction network set and / or the prediction network to suppress overfitting of the model when the number of training samples is less than the number of features.

[0080] Specifically, in crop genome prediction scenarios, the number of training samples typically ranges from hundreds to thousands of crop individuals, while the number of input features, i.e., SNP loci, can reach tens of thousands or even hundreds of thousands. This situation, where the number of features far exceeds the number of samples, is known in statistics and machine learning as the "large feature number, small sample number problem." Under these conditions, deep neural networks containing a large number of learnable parameters are prone to overlearning noise or individual-specific information in the training data. This leads to the model performing well on the training set but experiencing a significant drop in prediction accuracy on independent validation sets or candidate populations that were not used in the training, a phenomenon known as overfitting. Overfitting is one of the core challenges faced by deep learning models for genome prediction, directly limiting the model's cross-population generalization ability from the training population to the breeding population.

[0081] To address the aforementioned issues, during pre-training, L2 regularization constraints are applied to the model parameters of either or both of the heterogeneous feature extraction network set and the prediction network. L2 regularization, also known as weight decay, works by adding an additional penalty term to the original loss function. The magnitude of this penalty term is proportional to the sum of the squares of all constrained parameters.

[0082] In practice, the first step is to determine the range of parameters for which regularization constraints need to be applied. The learnable parameters in the heterogeneous feature extraction network set include the weight matrices of each fully connected layer in the first branch, the kernel weights of each convolutional layer in the second branch, and the query weight matrix, key weight matrix, and value weight matrix of the self-attention layer in the third branch. The learnable parameters of the prediction network include the weight matrices of each fully connected layer. During implementation, regularization can be applied to all of the above parameters, or it can be differentiated based on the number of parameters in each branch. For example, regularization can be applied only to the convolutional and self-attention branches with larger parameter numbers, while no regularization or a weaker regularization can be applied to the linear branches with relatively smaller parameter numbers.

[0083] The L2 regularization penalty is calculated as follows: The value of each of the selected learnable parameters subject to the constraint is squared, then the sum of all squared values ​​is multiplied by a preset regularization strength coefficient to obtain the L2 regularization penalty value. This penalty value is always non-negative, reaching its theoretical minimum of zero when all constrained parameters are zero. The regularization strength coefficient is a crucial hyperparameter, directly affecting the strength of regularization. If the coefficient is set too high, the model parameters will be over-compressed and approach zero, leading to underfitting and an inability to fully learn the mapping from genotype to phenotype. If the coefficient is set too low, the regularization effect is weak and cannot effectively suppress overfitting. In practice, this coefficient is usually determined through grid search or cross-validation on the validation set, with a common value range between 0.00001 and 0.01.

[0084] The final total training loss consists of two parts: one is the calculated composite loss function value, which is the product of the first weight coefficient and the mean squared error loss value plus the product of the second weight coefficient and the correlation loss value; the other part is the L2 regularization penalty value mentioned above. During model training, the backpropagation algorithm aims to minimize this total loss value, while simultaneously calculating the gradients of both losses with respect to all learnable parameters of the model and updating the parameters accordingly.

[0085] The principle behind L2 regularization in suppressing overfitting is that by forcing model parameters as close to zero as possible, it limits the search range of the parameter space, causing the model to tend to learn a smoother mapping function with a balanced weight distribution. Under regularization constraints, the model will not assign an extremely large weight coefficient to a specific SNP site unless that site shows a stable association with the phenotype in a large number of training samples. This effectively suppresses the model's learning of individual-specific or noisy data in the training samples, instead extracting consistent genetic patterns across the entire training population, thereby improving the generalization prediction ability for new individuals.

[0086] In addition to L2 regularization, other regularization strategies can be combined during pre-training to further enhance the anti-overfitting effect. For example, a random neuron dropout layer can be inserted after batch normalization of the fully connected layers. In each training batch, the output of some neurons is randomly set to zero with a certain probability, forcing the network to learn redundant and robust feature representations and reducing complex cooperative adaptations between neurons. An early stopping strategy can also be implemented, continuously monitoring the prediction accuracy on the validation set during training. When the validation set accuracy no longer improves in several consecutive training epochs, training is terminated early, retaining the model parameters at the moment of highest validation set accuracy, thus preventing the model from continuing to over-optimize on the training set. These strategies can be used in conjunction with L2 regularization to form a multi-layered regularization system, jointly suppressing overfitting from different perspectives.

[0087] The above description is merely a preferred embodiment of this application and an explanation of the technical principles employed. Those skilled in the art should understand that the scope of the invention involved in this application is not limited to technical solutions formed by specific combinations of the above-described technical features, but should also cover other technical solutions formed by arbitrary combinations of the above-described technical features or their equivalents without departing from the inventive concept. For example, technical solutions formed by substituting the above features with (but not limited to) technical features with similar functions disclosed in this application.

Claims

1. A method for predicting wheat genotype to phenotype based on heterogeneous deep learning, characterized in that, Includes the following steps: Obtain the genotype data of the crop to be tested, input the genotype data into a pre-trained prediction model, and obtain the phenotypic prediction results; The prediction model obtains the phenotypic prediction results in the following manner: The heterogeneous feature extraction network set is used to extract multiple heterogeneous features with different biological significance from the genotype data. The multiple heterogeneous features include at least a first feature reflecting additive genetic effects and a second feature reflecting nonlinear interaction effects. By using an adaptive feature fusion network, multiple heterogeneous features are dynamically weighted and fused with learnable fusion parameters to obtain fused features; The adaptive feature fusion network obtains the learnable fusion parameters in the following way: the multiple heterogeneous features are concatenated and input into a weight generation subnetwork. The weight generation subnetwork is then passed through a dimensionality reduction fully connected layer, a non-linear activation function, a fully connected layer that maps to a score equal to the number of features, and Softmax normalization to obtain the learnable weight coefficients corresponding to each feature. Each feature is multiplied by its corresponding learnable weight coefficient and summed to obtain the fused feature. Based on the fusion features, the phenotypic prediction results are output through the prediction network; The method of extracting multiple heterogeneous features from the genotype data using a heterogeneous feature extraction network set includes the following steps: The genotype data is simultaneously input into the first, second, and third parallel branches, respectively. In the first branch, the genotype data is linearly mapped, and a linear feature vector representing the additive genetic effect is extracted as the first feature; In the second branch, local pattern capture is performed on the genotype data to extract local nonlinear feature vectors that characterize local linkage disequilibrium block patterns, which are used as a component of the second feature. In the third branch, global context capture is performed on the genotype data to extract global context feature vectors that characterize the interaction relationships between distant sites across the entire genome, as another component of the second feature; The linear feature vector output from the first branch, the local nonlinear feature vector output from the second branch, and the global context feature vector output from the third branch are input into an adaptive feature fusion network, and the three are fused into a unified fusion feature vector through learnable fusion parameters.

2. The wheat genotype-to-phenotype prediction method based on heterogeneous deep learning according to claim 1, characterized in that, The dynamic weighted fusion of the multiple heterogeneous features using an adaptive feature fusion network includes the following steps: Assign a first learnable weight coefficient to the linear feature vector output by the first branch, assign a second learnable weight coefficient to the local nonlinear feature vector output by the second branch, and assign a third learnable weight coefficient to the global context feature vector output by the third branch. Multiply the linear feature vector by the first learnable weight coefficient to obtain a weighted linear feature vector; Multiply the local nonlinear feature vector by the second learnable weight coefficient to obtain a weighted local nonlinear feature vector; Multiply the global context feature vector by the third learnable weight coefficient to obtain the weighted global context feature vector; The fused feature vector is obtained by summing the weighted linear feature vector, the weighted local nonlinear feature vector, and the weighted global context feature vector element by element.

3. The wheat genotype-to-phenotype prediction method based on heterogeneous deep learning according to claim 1, characterized in that, In the first branch, performing a linear mapping on the genotype data and extracting a linear feature vector representing the additive genetic effect includes the following steps: The input genotype data is standardized to obtain a standardized genotype matrix; The standardized genotype matrix is ​​input into one or more stacked fully connected layers, wherein each fully connected layer does not have a non-linear activation function, and an initial linear feature vector is obtained through layer-by-layer linear transformation; The initial linear feature vector output by each fully connected layer is subjected to batch normalization and layer scaling in sequence. The batch normalization process is used to adjust the distribution of the initial linear feature vector to zero mean and unit variance. The layer scaling process is used to multiply the normalized feature vector with a learnable scaling factor to adjust the contribution of the layer output to the final linear feature vector. The feature vectors after batch normalization and layer scaling are used as the linear feature vectors output by the first branch.

4. The wheat genotype-to-phenotype prediction method based on heterogeneous deep learning according to claim 1, characterized in that, In the second branch, local pattern capture is performed on the genotype data to extract local nonlinear feature vectors representing local linkage disequilibrium block patterns, including the following steps: The genotype data is adjusted to a one-dimensional sequence format of single nucleotide polymorphism sites to obtain a one-dimensional input sequence; The one-dimensional input sequence is input into one or more stacked one-dimensional convolutional layers. In each one-dimensional convolutional layer, a set of learnable convolutional kernels slides along the positional arrangement direction of the one-dimensional input sequence to perform convolution operations on the positional segments within each sliding window, and a non-linear activation function is applied to extract local feature maps that reflect local chain imbalance block patterns. Pooling is performed on the local feature map output by each one-dimensional convolutional layer. The feature dimension is compressed by downsampling, which preserves salient features and expands the receptive field of subsequent convolutional layers. The feature map after the last pooling layer is flattened into a one-dimensional vector, which is used as the local nonlinear feature vector output by the second branch.

5. The wheat genotype-to-phenotype prediction method based on heterogeneous deep learning according to claim 1, characterized in that, In the third branch, global context capture is performed on the genotype data to extract global context feature vectors representing interactions between long-distance sites across the entire genome, including the following steps: The genotype data is converted into an embedding representation vector sequence, and the position encoding is superimposed on the embedding representation vector of each site to obtain an input embedding sequence with position information; The input embedding sequence with location information is input into one or more stacked self-attention layers. In each self-attention layer, the attention weight between any two sites in the input embedding sequence is calculated, and the embedding representation vectors of all sites are weighted and aggregated to capture the interaction relationship between sites at any distance across the whole genome, thus obtaining the global context representation sequence. The global context representation sequence output by the last self-attention layer is subjected to global average pooling or the representation vector of the first position of the sequence is taken to obtain the global context feature vector output by the third branch.

6. The wheat genotype-to-phenotype prediction method based on heterogeneous deep learning according to claim 1, characterized in that, The step of outputting the phenotypic prediction result through a prediction network based on the fused features includes the following steps: The fused feature vector is input into one or more stacked fully connected layers, each with a non-linear activation function. Through layer-by-layer non-linear mapping, the fused feature vector is transformed into a high-level abstract representation. The high-level abstract representation output from the last fully connected layer is input into an output fully connected layer without an activation function, and the phenotypic prediction result is obtained through linear mapping.

7. The wheat genotype-to-phenotype prediction method based on heterogeneous deep learning according to claim 1, characterized in that, The pre-training process of the prediction model includes the following steps: Obtain the genotype data and corresponding true phenotypic values ​​of the training population; The genotype data of the training population is input into the prediction model to be trained to obtain the predicted phenotypic values; A composite loss function is constructed, which is composed of a weighted sum of a mean squared error loss term and a correlation loss term based on the Pearson correlation coefficient. The mean squared error loss term is used to measure the magnitude of the error between the predicted phenotypic value and the true phenotypic value, and the correlation loss term is used to measure the degree of linear correlation between the predicted phenotypic value and the true phenotypic value. Based on the composite loss function, the loss value is calculated, and the parameters of the prediction model are updated through the backpropagation algorithm until the preset training termination condition is met, thus obtaining the pre-trained prediction model.

8. The wheat genotype-to-phenotype prediction method based on heterogeneous deep learning according to claim 7, characterized in that, The expression for the composite loss function is: L=α×L mse +β×L corr Where L is the total loss value, L mse L is the mean square error between the predicted phenotypic value and the true phenotypic value. corr The correlation loss term based on the Pearson correlation coefficient is denoted by α and β, which are preset first and second weighting coefficients, respectively. L corr The calculation method is as follows: L corr =1-|r|, where r is the Pearson correlation coefficient between the predicted phenotypic value and the actual phenotypic value.

9. The wheat genotype-to-phenotype prediction method based on heterogeneous deep learning according to claim 7, characterized in that, During the pre-training process, L2 regularization constraints are applied to the model parameters of the heterogeneous feature extraction network set and / or the prediction network to suppress overfitting of the model when the number of training samples is less than the number of features.

Images