Gene-environment interaction effect characterization prediction method based on neural network

Abstract

The application discloses a gene-environment interaction effect characterization prediction method based on a neural network, and collects a whole genome SNP typing matrix of a target crop population and time-series meteorological and static soil data in a growth period; an independent gene-environment dual-channel coding network is constructed, and a gene embedding vector and an environment embedding vector are extracted respectively; an explicit learnable interaction tensor is constructed based on the two types of embeddings, compressed gene-environment interaction embedding vectors are obtained through low-rank decomposition; the gene, environment and interaction embeddings are input into a multi-branch prediction head in parallel, target trait prediction values are output, an end-to-end joint training is carried out by using a composite loss function containing a prediction loss and a variance decomposition constraint, and finally, gene-environment interaction effect prediction results are obtained. The method has high precision, is consistent with genetic priori, and supports post-quantitative gene-environment interaction attribution analysis.

Description

Technical Field

[0001] This invention belongs to the field of agricultural information technology, specifically relating to a method for characterizing and predicting gene-environment interaction effects based on neural networks. Background Technology

[0002] Complex agronomic traits such as crop yield are constrained by both genes and environment. The yield forecasting problem in traditional breeding practices can be decomposed into three mutually coupled effect components: (1) gene main effect - genetic variance caused by differences in allele frequencies; (2) environmental main effect - phenotypic differences caused by changes in environmental factors such as precipitation, temperature, and soil; (3) gene-environment interaction effect - some genotypes show unexpected advantages or disadvantages under specific environmental conditions, which cannot be predicted by the main effect superposition model. Existing forecasting methods based on linear models (such as the genome best linear unbiased prediction GBLUP) can only capture the main effect, and the interaction effect is treated as a general random residual, resulting in the prediction accuracy often hovering between 0.45 and 0.55, which cannot meet the needs of precise selection.

[0003] In recent years, deep learning methods have emerged as promising approaches in phenotypic prediction. However, current implementations generally employ implicit interaction modeling—concatenating gene and environmental features before feeding them into a multilayer perceptron. This results in complete entanglement of the three types of effects, making it impossible to quantitatively separate the contributions of each effect or perform post-hoc interpretability analysis on the prediction results. Furthermore, multivariate meteorological time-series data (daily temperature, precipitation, sunshine, etc.) and static soil geographic data (soil pH, organic matter content, etc.) are heterogeneous, and confusing encoding weakens the representational power of environmental features. Moreover, deep networks without explicit structured modeling often have tens of millions of parameters, making them prone to overfitting to limited field trial data (typically only a few thousand to tens of thousands of samples).

[0004] The core technological bottleneck of existing genotype-environment interaction prediction lies in how to model the three types of effects in an explicit, decomposable, and interpretable manner within a parameter-efficient deep learning framework, so that the prediction results are both accurate and have quantifiable genetic physical significance, thereby supporting precise selection and environmental adaptability assessment in breeding decisions. Summary of the Invention

[0005] The purpose of this invention is to provide a neural network-based method for characterizing and predicting gene-environment interaction effects. This method achieves decoupled learning and quantitative attribution of gene major effects, environmental major effects, and gene-environment interaction effects, ensuring prediction accuracy while providing clear genetic interpretability of the prediction results.

[0006] To achieve the above-mentioned objectives, the specific technical solution is as follows: A neural network-based method for characterizing and predicting gene-environment interaction effects includes the following steps: Step S1: Collect genotype data and environmental data of the target crop population. The genotype data is the whole genome SNP genotyping matrix of the candidate samples, and the environmental data includes the time-series meteorological data and static soil geographic data of each experimental site.

[0007] Step S2: Construct a dual-channel decoupled coding network, which includes a gene channel encoder and an environmental channel encoder that are independent of each other. The gene channel encoder encodes the whole genome SNP genotyping matrix to obtain a gene embedding vector. The environmental channel encoder encodes and fuses the time-series meteorological data and static soil geographic data during the growth period to obtain an environmental embedding vector.

[0008] Step S3: Construct an explicit learnable interaction tensor based on the gene embedding vector and the environment embedding vector, and perform low-rank decomposition on the explicit learnable interaction tensor to obtain a compressed gene-environment interaction embedding vector.

[0009] Step S4: Input the gene embedding vector, the environment embedding vector, and the gene-environment interaction embedding vector into the multi-branch prediction head in parallel, output the predicted value of the target trait, and perform end-to-end joint training on the dual-channel decoupled coding network, the explicit learnable interaction tensor, and the multi-branch prediction head through a composite loss function that includes a prediction loss term and a variance decomposition constraint term, to obtain the gene-environment interaction effect prediction result of the target trait.

[0010] Furthermore, the gene channel encoder described in step S2 is implemented as follows: The whole genome SNP genotyping matrix The genome is divided into several segments based on its physical location on chromosomes. Each segment is first processed through a one-dimensional convolutional layer to extract local site association patterns, then through a self-attention layer to capture long-range linkage disequilibrium relationships across segments, and finally, the gene embedding vector is obtained through global average pooling. ;in, The gene embedding vector has a dimension of . ; This is the set of learnable parameters for the gene channel encoder. The value range is [64, 512].

[0011] Furthermore, the environmental channel encoder described in step S2 adopts a timing-static dual-sub-channel structure, specifically: The time-series meteorological data of the reproductive period The data is a multivariate meteorological time series collected daily, including daily average temperature, maximum temperature, minimum temperature, precipitation, sunshine duration, and relative humidity, covering the entire growth period from the sowing date. The temporal sub-channels are encoded using a bidirectional long short-term memory network to obtain the temporal environment embedding. .

[0012] The static soil geographic data This includes the latitude and longitude of the test site, altitude, soil pH, soil organic matter content, and soil texture grade; the static sub-channels are encoded using a multilayer perceptron to obtain the static environment embedding: .

[0013] The temporal environment embedding and the static environment embedding are concatenated as follows to obtain the environment embedding vector: ;in, The environment embedding vector has a dimension of . , The value range is [64, 512]; and These are the learnable parameters for the timing sub-channel and the static sub-channel, respectively.

[0014] Furthermore, the specific method for constructing the explicit learnable interaction tensor and its low-rank decomposition in step S3 is as follows: Performing the outer product operation between the gene embedding vector and the environment embedding vector yields a second-order explicit learnable interaction tensor: ;in, For explicit learnable interaction tensors This represents the outer product of vectors.

[0015] Furthermore, to reduce the number of parameters and avoid overfitting, the explicit learnable interaction tensor is parameterized using a CP low-rank decomposition to obtain a compressed gene-environment interaction embedding vector: ;in, Gene-environment interaction embedding vector; Let be the rank of the CP decomposition, with a value range of [8, 64]. and The first and second sides respectively, the gene side and the environment side There are learnable projection vectors, where ⊙ represents element-wise multiplication.

[0016] Through the aforementioned low-rank decomposition form of CP, the original size of Full interaction tensor compression to 3D vectors, the parameter complexity of interaction modeling is from Down to .

[0017] Furthermore, in step S4, the multi-branch prediction head adopts a principal-effect-interaction parallel structure. The gene embedding vector, the environment embedding vector, and the gene-environment interaction embedding vector are linearly summed after passing through independent fully connected branches to obtain the predicted value of the target trait. ;in, The predicted value for the target trait; , , These are the learnable weight vectors for the gene major effect branch, the environment major effect branch, and the gene-environment interaction branch, respectively. The bias scalar; the parallel structure ensures that the gene main effect, the environment main effect, and the gene-environment interaction effect are decoupled in the forward propagation path, which facilitates the application of subsequent variance decomposition constraints and the ex-post attribution analysis of each effect component.

[0018] Furthermore, the composite loss function containing the prediction loss term and the variance decomposition constraint term described in step S4 is constructed as follows: The first term is the prediction loss term, which calculates the mean squared error between the predicted value and the true phenotypic value for each sample in each training batch: ;in, This represents the number of samples in the training batch. For the first The true phenotypic value of each sample.

[0019] The second term is the variance decomposition constraint term. For the predicted values ​​of the current training batch, based on the parallel decoupling characteristics of the main effect-interaction branches, the contributions of the gene main effect branch, the environment main effect branch, and the gene-environment interaction branch to the variance of the final predicted value are calculated separately: ;in, , , These represent the sample variances of the gene major effect branch, the environment major effect branch, and the gene-environment interaction branch outputs within the batch.

[0020] Based on the above variance decomposition relationship, a variance consistency constraint is constructed as a regularization term to minimize the relative deviation between the sum of the three variances and the total prediction variance: .

[0021] Furthermore, the composite loss function is composed of a weighted average of the prediction loss term, the variance decomposition constraint term, and the parameter regularization term: ;in, The variance decomposition constraint weights have a value range of [0.01, 1.0]. The parameter is the regularization weight, and its value ranges from [1e-5, 1e-3]. This refers to the set of all learnable parameters in the dual-channel decoupled coding network, the explicit learnable interaction tensor, and the multi-branch prediction head.

[0022] By minimizing the composite loss function, the three parallel branches are forced to maintain low correlation, so that the variance of each branch output approximately corresponds to the true variance contribution of gene main effect, environment main effect and gene-environment interaction effect, respectively, thereby giving the gene-environment interaction effect prediction results a physical meaning that can be quantitatively attributed.

[0023] Furthermore, the explicit learnable interaction tensor and its low-rank decomposition described in step S3 can also be implemented using an environment-gated form. Specifically, the environment-gated form involves element-wise multiplying the gene-environment interaction embedding vector by the gating vector obtained after the environment embedding vector undergoes a gating transformation. ;in, and These are the learnable weight matrix and bias vector for environment gating, respectively. This is the Sigmoid activation function, and ⊙ represents element-wise multiplication.

[0024] The described environmental gating method realizes the physical meaning of "environmental modulation gene effect" in a way with lower parameter complexity, and is suitable for implementation scenarios with high environmental data dimensionality and large explicit tensor computation overhead.

[0025] Furthermore, the explicit learnable interaction tensor and its low-rank decomposition described in step S3 can also be replaced by a multi-head cross-attention form, which specifically involves: dividing the gene embedding vector into gene-side sequences according to genomic fragments and dividing the environment embedding vector into environment-side sequences according to reproductive stage, using the gene-side sequences as queries and the environment-side sequences as keys and values.

[0026] Calculate multi-head cross-attention: ;in, , , The first The query, key, and value matrix of each attention head; The key vector dimension for each attention head is defined, and the number of attention heads ranges from [2, 16]. The gene-environment interaction embedding vector is obtained by concatenating the outputs of each attention head and then linearly projecting the result.

[0027] The multi-head cross-attention approach is suitable for fine-grained interaction scenarios that require modeling "specific genomic segments being stimulated by the environment at specific reproductive stages".

[0028] Furthermore, for target crops or target traits with limited training sample size, the method also includes a cross-scenario transfer training step: The dual-channel decoupled coding network and the explicit learnable interaction tensor are pre-trained on a large-scale public breeding dataset to obtain source domain pre-training parameters. When fine-tuning on a small-scale training set of the target crop or target trait, all parameters of the gene channel encoder and the environment channel encoder except the last layer are frozen, and only the last layer, the explicit learnable interaction tensor, and the multi-branch prediction head are updated.

[0029] Furthermore, a KL divergence constraint term for the source-target domain parameter distribution is added to the composite loss function: ;in, The loss is calculated using the composite loss function over the target domain; The KL divergence between the parameter distributions of the source domain and the parameter distributions of the target domain; The migration constraint weights have a value range of [0.01, 0.5].

[0030] Through the cross-scenario transfer training step, even when the number of training samples in the target domain is less than 10% of the number of training samples in the source domain, the prediction performance can still be maintained without significant degradation.

[0031] Compared with the prior art, the beneficial effects of this invention are: The prediction accuracy of this invention is significantly improved. In wheat grain yield prediction, the Pearson correlation coefficient reaches 0.74, which is about 0.22 higher than the GBLUP method and about 0.13 higher than the standard MLP. By using explicit interaction tensors and variance decomposition constraints, the interaction effects are decoupled from the main effects of genes and environment. The measured variance contributions of the three types of effects are consistent with the genetic priors of many years of field trials. For small sample target traits, cross-scenario transfer training can maintain performance without degradation even when the sample size is only 10% of the source domain. Attached Figure Description

[0032] Figure 1 This is a flowchart of the gene-environment interaction effect characterization and prediction method based on neural networks of the present invention. Detailed Implementation

[0033] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention are described clearly and completely below. Obviously, the described embodiments are only a part of the embodiments of this invention, not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0034] Example 1 like Figure 1The diagram shows the flowchart of the gene-environment interaction effect characterization and prediction method based on neural networks in this embodiment. Using wheat grain yield (kg / mu) as the target trait, gene-environment interaction effects were modeled and predicted for 600 wheat inbred lines from 8 ecological experimental sites in northern winter wheat and southern spring wheat regions across nearly 5 growth cycles. The training sample size in this embodiment is approximately 4800 records (600 inbred lines × 8 experimental sites, including some missing records), which are divided into training, validation, and test sets in a 7:1:2 ratio. The training set is further used for a complete end-to-end training process.

[0035] The method includes the following steps S1 to S4.

[0036] Step S1: Collect genotype data and environmental data of the target crop population. The genotype data is the whole genome SNP genotyping matrix of the candidate samples, and the environmental data includes the time-series meteorological data and static soil geographic data of each experimental site.

[0037] Specifically, for the wheat grain yield prediction scenario, the data collection method is as follows: (1) Genotypic data: 600 wheat inbred lines were genotyped using the 90KiSelect SNP chip. Quality control was performed (minor allele frequency MAF>0.05, deletion rate<10%, Hardy-Weinberg equilibrium test). Approximately 82,000 high-quality SNP loci were retained, encoded additively (0 for reference homozygosity, 1 for heterozygosity, and 2 for alternative homozygosity), forming a genotype matrix. .

[0038] (2) Meteorological data during the growing season were collected daily by automatic weather stations at each experimental site, covering the entire growing season of 270 days from the sowing date (around early October) to the harvest date (around the end of June). Each time step included six variables: daily average temperature, daily maximum temperature, daily minimum temperature, daily precipitation, sunshine duration, and relative humidity, forming a meteorological time series tensor XE^seq∈R^(8×270×6); (3) Static soil geographic data included eight dimensions: latitude and longitude of the experimental site (2-dimensional), altitude (1-dimensional), soil pH value (1-dimensional), soil organic matter content (1-dimensional), and soil texture grade (3-dimensional according to one-hot encoding of sandy soil / loam / clay soil), forming a static environment matrix. (4) The target trait y is the measured grain yield of each sample at the corresponding test point, in kg / mu.

[0039] Step S2: Construct a dual-channel decoupled coding network, which includes a gene channel encoder and an environmental channel encoder that are independent of each other. The gene channel encoder encodes the whole genome SNP genotyping matrix to obtain a gene embedding vector. The environmental channel encoder encodes and fuses the time-series meteorological data and static soil geographic data during the growth period to obtain an environmental embedding vector.

[0040] The gene channel encoder is implemented by: converting the whole genome SNP genotyping matrix... The genome is divided into several segments based on its physical location on chromosomes. Each segment is first processed through a one-dimensional convolutional layer to extract local site association patterns, then through a self-attention layer to capture long-range linkage disequilibrium relationships across segments, and finally, the gene embedding vector is obtained through global average pooling. ;in, The gene embedding vector has a dimension of . ; This is the set of learnable parameters for the gene channel encoder. The value range is [64, 512].

[0041] In this embodiment, 82,000 SNP loci are divided into 21 genomic segments based on their physical locations on 21 wheat chromosomes (21 pairs of chromosomes, 1A-7D). Each segment contains an average of approximately 3,900 SNP loci. Each segment is first processed through three layers of one-dimensional convolution (kernel sizes of 7, 5, and 3, output channel numbers of 32, 64, and 128, and the activation function is ReLU) to extract local locus association patterns. The convolutional outputs of the 21 segments are then concatenated along the segment dimension and fed into two layers of self-attention layers (8 attention heads and 128 model dimensions) to capture long-range linkage disequilibrium relationships across chromosomes. Finally, global average pooling is used to obtain the gene embedding vector. Dimension =256. The total number of parameters in the gene channel encoder is approximately 3.2 million.

[0042] In other implementation schemes, Adjustable within the range of [64,512]: A smaller value (e.g., 64) is suitable for scenarios with a small sample size and a high risk of overfitting. Larger values ​​(e.g., 512) are suitable for scenarios with sufficient sample size and complex trait genetic structure.

[0043] The environmental channel encoder adopts a time-series-static dual-channel structure, specifically: the time-series meteorological data during the reproductive period The data is a multivariate meteorological time series collected daily, including daily average temperature, maximum temperature, minimum temperature, precipitation, sunshine duration, and relative humidity, covering the entire growth period from the sowing date. The temporal sub-channels are encoded using a bidirectional long short-term memory network to obtain the temporal environment embedding. .

[0044] The static soil geographic data This includes the latitude and longitude of the test site, altitude, soil pH, soil organic matter content, and soil texture grade; the static sub-channels are encoded using a multilayer perceptron to obtain the static environment embedding: .

[0045] The temporal environment embedding and the static environment embedding are concatenated as follows to obtain the environment embedding vector: ;in, The environment embedding vector has a dimension of . , The value range is [64, 512]; and These are the learnable parameters for the timing sub-channel and the static sub-channel, respectively.

[0046] In this embodiment, the timing sub-channels employ a two-layer bidirectional long short-term memory network (128 hidden units per direction, therefore...). The input is a 256-dimensional meteorological time-series tensor representing the wheat growth period with 270 days × 6 variables. To enhance the impact of key growth stages (wheat overwintering, jointing, heading, and grain-filling stages), a time-dimensional attention pooling layer is superimposed on the BiLSTM output to obtain a fixed-length temporal context embedding. (Dimension 256).

[0047] The static subchannel employs a 3-layer multilayer perceptron (hidden layer dimensions of 64, 128, and 64 respectively, with ReLU activation function) to encode the 8-dimensional static soil geographic vector into a static environment embedding. (Dimension 64).

[0048] The two are concatenated to obtain the environment embedding vector. Dimension =256+64=320. In other implementations, dE can be adjusted in the range of [64,512].

[0049] Step S3: Construct an explicit learnable interaction tensor based on the gene embedding vector and the environment embedding vector, and perform low-rank decomposition on the explicit learnable interaction tensor to obtain a compressed gene-environment interaction embedding vector.

[0050] The specific method for constructing and decomposing the explicit learnable interaction tensor is as follows: The gene embedding vector and the environment embedding vector are multiplied by an outer product to obtain a second-order explicit learnable interaction tensor. ;in, For explicit learnable interaction tensors This represents the outer product of vectors.

[0051] To reduce the number of parameters and avoid overfitting, the explicit learnable interaction tensor is parameterized using a CP low-rank decomposition to obtain a compressed gene-environment interaction embedding vector: ;in, Gene-environment interaction embedding vector; Let be the rank of the CP decomposition, with a value range of [8, 64]. and The first and second sides respectively, the gene side and the environment side A learnable projection vector, where ⊙ represents element-wise multiplication; through the aforementioned CP low-rank decomposition form, the original size is... Full interaction tensor compression to 3D vectors, the parameter complexity of interaction modeling is from Down to .

[0052] In this embodiment, the outer product formula constructs an explicit learnable interaction tensor. The elements of this tensor Characterizing gene embedding Embedded in the environment Second-order interactions between dimensions.

[0053] If all elements of TGE (256×320=81920 dimensions) are directly used as subsequent input, it will lead to an explosion of prediction head parameters and easy overfitting; therefore, CP low-rank decomposition is used to parameterize and compress TGE.

[0054] Specifically, in this embodiment, the rank R of the CP decomposition is set to 32, corresponding to 32 sets of learnable projection vectors. Each of them , That is, gene-side projection matrix Environmental side projection matrix Gene-environment interaction embedding vectors were calculated using the CP decomposition formula. .

[0055] The CP low-rank decomposition reduces the parameter complexity of interaction modeling from... Compress to The number of parameters is reduced by approximately 77%; at the same time, the interpretability of explicit interaction modeling is preserved. Each dimension can be traced back to a specific gene-side projection direction. Projection direction with the environment The coupling. In other implementations, the rank It can be adjusted within the range of [8, 64].

[0056] To further illustrate the specific calculation process of the CP low-rank decomposition described in this embodiment, a complete forward propagation example is given: Suppose a sample in the current batch has a gene embedding vector The values ​​of the first 5 dimensions are , Environment Embedding Vector The values ​​of the first 5 dimensions are ; the first 5 dimensions of the first group projection vector u1 on the gene side and environmental side Group projection vectors The first 5 dimensions For example, then (Approximated by the first 5 dimensions) ≈ (The actual calculation is a full 256-dimensional dot product; only the first 5 dimensions are shown here.) v1 zE (approximated by the first 5 dimensions) ≈ 0.387; according to the CP decomposition formula, The first dimension of ×E is ( )·( zE)≈-0.011×0.387≈-0.0043; Traverse in the same way A complete 32-dimensional interaction embedding vector can be obtained. ×E.

[0057] Step S4: Input the gene embedding vector, the environment embedding vector, and the gene-environment interaction embedding vector into the multi-branch prediction head in parallel, output the predicted value of the target trait, and perform end-to-end joint training on the dual-channel decoupled coding network, the explicit learnable interaction tensor, and the multi-branch prediction head through a composite loss function that includes a prediction loss term and a variance decomposition constraint term, to obtain the gene-environment interaction effect prediction result of the target trait.

[0058] The multi-branch prediction head adopts a principal-effect-interaction parallel structure, specifically: ;in, The predicted value for the target trait; , , These are the learnable weight vectors for the gene major effect branch, the environment major effect branch, and the gene-environment interaction branch, respectively. The bias scalar; the parallel structure ensures that the gene main effect, the environment main effect, and the gene-environment interaction effect are decoupled in the forward propagation path, which facilitates the application of subsequent variance decomposition constraints and the ex-post attribution analysis of each effect component.

[0059] In this embodiment, the major branch of the gene is a single-layer fully connected ( The environment's main branch is a single-layer fully connected branch. The gene-environment interaction branch is a single-layer fully connected layer. The output dimensions are all scalars; the bias scalar b is randomly initialized to zero. During forward propagation, the three branches compute their respective scalar outputs in parallel to obtain the predicted wheat grain yield. (kg / mu).

[0060] In a typical sample of this embodiment, the outputs of the three branches are as follows: ≈245.3、 ≈108.6 ≈21.8、 ≈10.2, then ≈245.3+108.6+21.8+10.2=385.9 (kg / mu), of which the main contribution of gene effect is about 63.6%, the main contribution of environment effect is about 28.1%, the contribution of gene-environment interaction is about 5.7%, and the bias contribution is about 2.6%, which is consistent with the genetic prior of wheat yield prediction that "the main effect of gene effect is dominant, the main effect of environment effect is secondary, and the interaction effect is supplementary".

[0061] The composite loss function, which includes a prediction loss term and a variance decomposition constraint term, is constructed as follows: The first term is the prediction loss term, which calculates the mean squared error between the predicted value and the true phenotypic value for each sample in each training batch: ;in, This represents the number of samples in the training batch. For the first The true phenotypic value of each sample.

[0062] The second term is the variance decomposition constraint term. For the predicted values ​​of the current training batch, based on the parallel decoupling characteristics of the main effect-interaction branches, the contributions of the gene main effect branch, the environment main effect branch, and the gene-environment interaction branch to the variance of the final predicted value are calculated separately: ;in, , , These represent the sample variances of the gene major effect branch, the environment major effect branch, and the gene-environment interaction branch outputs within the batch.

[0063] Based on the above variance decomposition relationship, a variance consistency constraint is constructed as a regularization term to minimize the relative deviation between the sum of the three variances and the total prediction variance: The composite loss function is composed of a weighted average of the predicted loss term, the variance decomposition constraint term, and the parameter regularization term. ;in, The variance decomposition constraint weights have a value range of [0.01, 1.0]. The parameter is the regularization weight, and its value ranges from [1e-5, 1e-3]. This refers to the set of all learnable parameters in the dual-channel decoupled coding network, the explicit learnable interaction tensor, and the multi-branch prediction head.

[0064] By minimizing the composite loss function, the three parallel branches are forced to maintain low correlation, so that the variance of each branch output approximately corresponds to the true variance contribution of gene main effect, environment main effect and gene-environment interaction effect, respectively, thereby giving the gene-environment interaction effect prediction results a physical meaning that can be quantitatively attributed.

[0065] In this embodiment, the variance decomposition constraint weights are taken as follows: Parameter regularization weights Learnable parameter set This includes gene channel encoder parameters (approximately 3.2 million), environmental channel encoder parameters (approximately 800,000), and CP low-rank decomposition parameters. and Matrix (approximately 18,000), multi-branch prediction head weights (Approximately 600); The optimizer used was AdamW with an initial learning rate of 1×10⁻⁶. - ³, decay according to cosine annealing strategy; batch size The total number of training rounds is 200. The early stopping strategy is based on the condition that the predicted loss Lpred on the validation set does not decrease for 10 consecutive rounds.

[0066] In the actual training process of wheat in this embodiment, the variance decomposition constraint term Lvar gradually converged from an initial value of approximately 0.42 (i.e., the ratio of the sum of the variances of the three branch outputs to the total prediction variance deviated by approximately 42%) to below 0.05 at the end of training, indicating that the covariance among the outputs of the three parallel branches has been significantly suppressed, and the variance decomposition results of the main effect-interaction have physical consistency. In other embodiments, Available Adjustments within the scope: A smaller variance can weaken the variance decoupling constraint and improve prediction accuracy; A larger value indicates a greater emphasis on interpretability and facilitates post-hoc attribution analysis; Available Adjustments can be made within the specified range.

[0067] To verify the effectiveness of the method in this embodiment, the inventors conducted comparative experiments using the following three schemes on the same wheat dataset: Scheme A is the classic genome best linear unbiased prediction (GBLUP) method, which only captures gene major effects and does not explicitly model gene-environment interactions.

[0068] Scheme B is the standard multilayer perceptron (MLP) method, which concatenates the genotype matrix and environmental feature vectors and directly inputs them into a 5-layer fully connected network (implicit interaction modeling).

[0069] Scheme C is the method described in this embodiment. The three schemes are trained independently on the same training set / validation set / test set partition, predict wheat grain yield, and use 5-fold cross-validation to statistically analyze various indicators.

[0070] Regarding prediction accuracy, the Pearson correlation coefficients of schemes A, B, and C on the test set... The values ​​are 0.52, 0.61, and 0.74, respectively, with relative root mean square errors of 0.31, 0.27, and 0.21. The method in this embodiment improves upon GBLUP. With an absolute value of approximately 0.22 and an improvement of approximately 0.13 compared to MLP, the prediction accuracy has a significant advantage.

[0071] Regarding the consistency of variance decomposition, scheme A can only provide the variance contribution of the main gene effect but cannot separate gene-environment interaction; although scheme B has higher overall accuracy, the three types of effects are entangled due to implicit modeling, and the variance decomposition constraint term Lvar is measured to be 0.38 (that is, the sum of the three variances deviates from the total predicted variance by about 38%), resulting in poor physical consistency; after training, scheme C's Lvar converges to below 0.04, and the variance contributions of the main gene effect, the main environmental effect, and the gene-environment interaction are about 62%, 25%, and 13%, respectively, which are highly consistent with the empirical values ​​observed in wheat field trials over many years: "yield heritability is about 0.55-0.65, environmental effect is about 20%-30%, and interaction effect is about 10%-15%".

[0072] Regarding cross-environment generalization ability, the "leave one environment for cross-validation" method (leaving one ecological point out of training each time and predicting the new environment during testing) was used for evaluation. The r values ​​of schemes A, B, and C decreased to 0.38, 0.42, and 0.61, respectively. The method in this embodiment still maintains strong predictive ability in unseen new environments, verifying the generalization advantage brought by explicit modeling of gene-environment interactions.

[0073] The method in this embodiment has approximately 4.2 million learnable parameters, and the training time for a single round is approximately 18 seconds (for a single NVIDIA A100 GPU). The total training time for 200 rounds is approximately 1 hour, which is roughly the same as the standard MLP scheme and significantly lower than the brute-force implementation that directly constructs a complete outer product tensor (the brute-force implementation has approximately 82 million parameters, takes more than 8 hours to train, and suffers from severe overfitting).

[0074] Example 2 This embodiment replaces the explicit learnable interaction tensor and its CP low-rank decomposition step with an environment-gated implementation.

[0075] Specifically, in the wheat scenario: the gene-environment interaction embedding vector is obtained by element-wise multiplying the gating vector (obtained after gating the environment embedding vector) with the gene embedding vector. ;in, and These are the learnable weight matrix and bias vector for environment gating, respectively. This is the Sigmoid activation function, and ⊙ represents element-wise multiplication.

[0076] The described environmental gating method realizes the physical meaning of environmental modulation gene effects in a way with lower parameter complexity, and is suitable for implementation scenarios with high environmental data dimensionality and large explicit tensor computation overhead.

[0077] In the alternative implementation of the environment gating, the learnable weight matrix of the environment gating Bias vector During forward propagation, the environment is embedded into the vector. (320-dimensional) after linear transformation We obtain a 256-dimensional vector, which is then compressed to the [0,1] interval using the Sigmoid activation function to form a gated vector. (256-dimensional), then combined with the gene embedding vector Element-wise multiplication yields the gene-environment interaction embedding vector. (256 dimensions).

[0078] To illustrate the specific calculation process of the optional implementation of the environmental gating, a complete forward propagation example is given: Assume a sample in the current batch has a gene embedding vector... The values ​​of the first 5 dimensions are , Environment Embedding Vector The values ​​of the first 5 dimensions are .

[0079] by The first 5 rows and first 5 columns of the submatrix: Bias The first 5 dimensions For illustration purposes.

[0080] Approximate calculation of the first 5 dimensions of gated linear transformation (only taking...) (First 5 dimensions involved): Components ; quantity ; and so on, the components can be obtained. , quantity , quantity (The actual calculation is a full 320-dimensional dot product).

[0081] The result of the gated linear transformation is then activated by the Sigmoid function: , , , , Finally, sort by element and... Multiplying the first 5 dimensions: , , , , After traversing all 256 dimensions, the gene-environment interaction embedding vector under this optional implementation method is obtained.

[0082] The number of learnable parameters for the environmental gating method is: Compared to CP low-rank decomposition It has slightly more computational overhead, but does not require intermediate explicit tensor construction and has lower forward computational overhead; its physical meaning is "dimensional gating modulation of the main effect of the gene by environmental signal", which is suitable for implementation scenarios with high environmental data dimensionality.

[0083] Example 3 This embodiment also replaces the explicit learnable interaction tensor and its CP low-rank decomposition step with a multi-head cross-attention form.

[0084] Specifically, in the wheat scenario: the gene embedding vector is segmented into gene-side sequences according to genomic fragments, and the environment embedding vector is segmented into environment-side sequences according to the growth stage. Using the gene-side sequences as queries and the environment-side sequences as keys and values, multi-head cross-attention is calculated. ;in, , , The first The query, key, and value matrix of each attention head; The key vector dimension for each attention head is defined, and the number of attention heads ranges from [2, 16]. The gene-environment interaction embedding vector is obtained by concatenating the outputs of each attention head and then linearly projecting the result.

[0085] The multi-head cross-attention approach is suitable for fine-grained interaction scenarios that require modeling "specific genomic segments being stimulated by the environment at specific reproductive stages".

[0086] In the alternative implementation of the multi-head cross-attention method, this embodiment divides the gene embedding vector zG into gene-side sequences according to 21 wheat chromosome segments. (Each segment has a dimension of 256 / 21, rounded down and padded with zeros for alignment). The environmental embedding vector zE is segmented into environmental-side sequences according to the four key growth stages of wheat (overwintering, jointing, heading, and grain-filling). .

[0087] The multi-head cross-attention algorithm has 8 attention heads (h=8) and 64 dimensions for each key vector. The gene-side sequence, after linear projection, is used as the query Q, and the environment-side sequence, after linear projection, is used as the key K and value V. The multi-head cross-attention output is calculated, then pooled after a fully connected layer to obtain the gene-environment interaction embedding vector. .

[0088] The advantage of the multi-headed cross-attention approach lies in its ability to explicitly model the fine-grained interactions between specific genomic segments (such as segments on chromosomes carrying yield QTLs) and specific growth stages (such as high-temperature stress during the heading stage), thereby supporting post-hoc visual attribution through the attention weight matrix. The training process, composite loss function, and variance decomposition constraints of the optional implementation of the multi-headed cross-attention approach are completely consistent with the main approach, and the specific calculation process will not be elaborated here.

[0089] Example 4 Beyond wheat grain yield prediction, for target traits with small training sample sizes (for example, the measured sample sizes of the two traits "Fusarium head blight resistance score" and "grain protein content" in the wheat dataset described in this embodiment are only about 600 and 800 respectively), cross-scenario transfer training steps are implemented.

[0090] The dual-channel decoupled coding network and the explicit learnable interaction tensor are pre-trained on a large-scale public breeding dataset to obtain source domain pre-training parameters. When fine-tuning on a small-scale training set for the target crop or target trait, all parameters in the gene channel encoder and the environment channel encoder except for the last layer are frozen. Only the last layer, the explicit learnable interaction tensor, and the multi-branch prediction head are updated, and a KL divergence constraint term for the source-target domain parameter distribution is added to the composite loss function. ;in, The loss is calculated using the composite loss function over the target domain. The KL divergence between the parameter distributions of the source domain and the parameter distributions of the target domain; The transfer constraint weights have a value range of [0.01, 0.5]. Through the cross-scene transfer training step, even when the number of training samples in the target domain is less than 10% of the number of training samples in the source domain, the prediction performance can still be maintained without significant degradation.

[0091] Specifically, in this embodiment, a public maize breeding dataset containing approximately 48,000 samples (covering 10 ecopoints, 3 growth cycles, and approximately 1,600 maize inbred lines) is used as the source domain. The dual-channel decoupled coding network and the explicit learnable interaction tensor are pre-trained. Wheat data is not involved in the pre-training stage. When fine-tuning the wheat target domain, all parameters in the gene channel encoder except for the last self-attention layer are frozen, and all parameters in the BiLSTM and static MLP in the environment channel encoder are frozen. Only the last self-attention layer, the U and V matrices of the CP low-rank decomposition, and the multi-branch prediction head are updated.

[0092] migration constraint weights In the fine-tuning experiment of wheat scab resistance scoring traits, compared to training from scratch, cross-scenario transfer training increased the Pearson correlation coefficient on the test set from 0.41 to 0.58 and decreased the relative root mean square error from 0.32 to 0.24, significantly improving prediction performance. Similarly, for the wheat grain protein content trait, cross-scenario transfer training increased the Pearson correlation coefficient from 0.46 to 0.61. In other implementation schemes, the transfer constraint weights... Available Adjustments can be made within the specified range.

[0093] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for characterizing and predicting gene-environment interaction effects based on neural networks, characterized in that, The method includes the following steps: Step S1: Collect genotype data and environmental data of the target crop population. The genotype data is the whole genome SNP genotyping matrix of the candidate samples, and the environmental data includes the time-series meteorological data and static soil geographic data of each experimental site. Step S2: Construct a dual-channel decoupled coding network, which includes a gene channel encoder and an environmental channel encoder that are independent of each other. The gene channel encoder encodes the whole genome SNP genotyping matrix to obtain a gene embedding vector. The environmental channel encoder encodes and fuses the time-series meteorological data and static soil geographic data during the growth period to obtain an environmental embedding vector. Step S3: Construct an explicit learnable interaction tensor based on the gene embedding vector and the environment embedding vector, and perform low-rank decomposition on the explicit learnable interaction tensor to obtain a compressed gene-environment interaction embedding vector. Step S4: Input the gene embedding vector, the environment embedding vector, and the gene-environment interaction embedding vector into the multi-branch prediction head in parallel, output the predicted value of the target trait, and perform end-to-end joint training on the dual-channel decoupled coding network, the explicit learnable interaction tensor, and the multi-branch prediction head through a composite loss function that includes a prediction loss term and a variance decomposition constraint term, to obtain the gene-environment interaction effect prediction result of the target trait.

2. The method according to claim 1, characterized in that, The gene channel encoder described in step S2 is implemented by: converting the whole genome SNP genotyping matrix... The genome is divided into several segments based on its physical location on chromosomes. Each segment is first processed through a one-dimensional convolutional layer to extract local site association patterns, then through a self-attention layer to capture long-range linkage disequilibrium relationships across segments, and finally, the gene embedding vector is obtained through global average pooling. ;in, The gene embedding vector has a dimension of . ; This is the set of learnable parameters for the gene channel encoder. The value range is [64, 512].

3. The method according to claim 2, characterized in that, The environmental channel encoder described in step S2 adopts a timing-static dual-sub-channel structure, specifically: The time-series meteorological data of the reproductive period The data is a multivariate meteorological time series collected daily, including daily average temperature, maximum temperature, minimum temperature, precipitation, sunshine duration, and relative humidity, covering the entire growth period from the sowing date. The temporal sub-channels are encoded using a bidirectional long short-term memory network to obtain the temporal environment embedding. ; The static soil geographic data This includes the latitude and longitude of the test site, altitude, soil pH value, soil organic matter content, and soil texture grade; The static sub-channels are encoded using a multilayer perceptron to obtain the static environment embedding: ; The temporal environment embedding and the static environment embedding are concatenated as follows to obtain the environment embedding vector: ;in, The environment embedding vector has a dimension of . , The value range is [64, 512]; and These are the learnable parameters for the timing sub-channel and the static sub-channel, respectively.

4. The method according to claim 3, characterized in that, The specific method for constructing the explicit learnable interaction tensor and its low-rank decomposition in step S3 is as follows: Perform an outer product operation between the gene embedding vector and the environment embedding vector to obtain a second-order explicit learnable interaction tensor. ;in, For explicit learnable interaction tensors Represents the cross product of vectors; To reduce the number of parameters and avoid overfitting, the explicit learnable interaction tensor is parameterized using a CP low-rank decomposition to obtain a compressed gene-environment interaction embedding vector: ;in, Gene-environment interaction embedding vector; Let be the rank of the CP decomposition, with a value range of [8, 64]. and The first and second sides respectively, the gene side and the environment side A learnable projection vector, where ⊙ represents element-wise multiplication; through the aforementioned CP low-rank decomposition form, the original size is... Full interaction tensor compression to 3D vectors, the parameter complexity of interaction modeling is from Down to .

5. The method according to claim 4, characterized in that, In step S4, the multi-branch prediction head adopts a principal-interaction parallel structure. The gene embedding vector, the environment embedding vector, and the gene-environment interaction embedding vector are linearly summed after passing through independent fully connected branches to obtain the predicted value of the target trait. ;in, The predicted value for the target trait; , , These are the learnable weight vectors for the gene major effect branch, the environment major effect branch, and the gene-environment interaction branch, respectively. The bias scalar; the parallel structure ensures that the gene main effect, the environment main effect, and the gene-environment interaction effect are decoupled in the forward propagation path, which facilitates the application of subsequent variance decomposition constraints and the ex-post attribution analysis of each effect component.

6. The method according to claim 5, characterized in that, The composite loss function in step S4, which includes the prediction loss term and the variance decomposition constraint term, includes the prediction loss term. Variance decomposition constraint terms and parameter regularization terms Weighted composition: ;in, The variance decomposition constraint weights have a value range of [0.01, 1.0]. The parameter is the regularization weight, and its value ranges from [1e-5, 1e-3]. This refers to the set of all learnable parameters in the dual-channel decoupled coding network, the explicit learnable interaction tensor, and the multi-branch prediction head.

7. The method according to claim 6, characterized in that, The prediction loss term calculates the mean squared error between the predicted value and the true phenotypic value for each sample in the training batch: ;in, This represents the number of samples in the training batch. For the first The true phenotypic values ​​of each sample; The variance decomposition constraint term, for the predicted values ​​of the current training batch, calculates the contributions of the gene main effect branch, the environment main effect branch, and the gene-environment interaction branch to the variance of the final predicted value, based on the parallel decoupling characteristic of the main effect-interaction branches: ;in, , , The variances of the gene major effect branch, the environment major effect branch, and the gene-environment interaction branch within the batch are respectively represented. Based on the above variance decomposition relationship, a variance consistency constraint is constructed as a regularization term to minimize the relative deviation between the sum of the variances of the gene major effect branch, the environment major effect branch, and the gene-environment interaction branch and the total prediction variance. .

8. The method for characterizing and predicting gene-environment interaction effects based on neural networks according to claim 4, characterized in that, The specific method for constructing the explicit learnable interaction tensor and its low-rank decomposition in step S3 can also be replaced by an environment-gated approach. Specifically, the environment-gated approach involves element-wise multiplying the gene-environment interaction embedding vector by the gating vector obtained after the environment embedding vector undergoes a gating transformation. ;in, and These are the learnable weight matrix and bias vector for environment gating, respectively. This is the Sigmoid activation function, and ⊙ represents element-wise multiplication.

9. The method according to claim 4, characterized in that, The specific method for constructing the explicit learnable interaction tensor and its low-rank decomposition in step S3 can also be replaced by a multi-head cross-attention approach. Specifically, the multi-head cross-attention approach involves: segmenting the gene embedding vector into gene-side sequences according to genomic fragments, segmenting the environment embedding vector into environment-side sequences according to reproductive stage, and using the gene-side sequences as queries and the environment-side sequences as keys and values ​​to calculate the multi-head cross-attention. ;in, , , The first The query, key, and value matrix of each attention head; The key vector dimension for each attention head is defined, and the number of attention heads ranges from [2, 16]. The gene-environment interaction embedding vector is obtained by concatenating the outputs of each attention head and then linearly projecting the result.

10. The method according to claim 9, characterized in that, For target crops or target traits with limited training sample size, the method further includes a cross-scenario transfer training step: pre-training the dual-channel decoupled coding network and the explicit learnable interaction tensor on a large-scale public breeding dataset to obtain source domain pre-training parameters; when fine-tuning on a small-scale training set of the target crop or target trait, freezing all parameters in the gene channel encoder and the environment channel encoder except for the last layer, and updating only the last layer, the explicit learnable interaction tensor, and the multi-branch prediction head.

Images