Key pathogenic factor screening method based on deep embedding and gene network constraint

Abstract

The application belongs to the technical field of biometric recognition, and relates to a key pathogenic factor screening method based on deep embedding and gene network constraint. The method first acquires a gene expression matrix; a gene interaction network corresponding to the gene expression matrix is constructed to generate a sparse adjacency matrix and a sparse mask matrix; a graph variational autoencoder model with gene network constraint is constructed, the input of the model is the gene expression matrix, the sparse mask matrix M is used in the form of element-level multiplication to limit the weight connection of the encoder, the decoder reconstructs the input gene expression matrix according to the latent structure characteristics, and a multi-task prediction module at the output end predicts the infection stage probability distribution and pathogen load through parallel classification branches and regression branches; the constructed model is trained; feature attribution analysis is performed on the input genes according to the task output of the model, the importance scores of the genes are calculated, and a key gene candidate set is generated. The method significantly improves the biological rationality and stability of the screening result.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of biological feature recognition, and particularly relates to a key pathogenic factor screening method based on deep embedding and gene network constraint, which is used for identifying key pathogenic factors of corn in the infection process of Fusarium graminearum. BACKGROUND

[0002] Corn is an important food and feed crop in the world, and its yield and quality have been seriously threatened by pathogenic fungi such as Fusarium graminearum for a long time. The pathogenic fungi can infect corn ears and grains, causing ear rot and scab, causing significant economic losses, and can also produce toxin pollution, threatening food safety. Revealing the defense response mechanism and key regulatory genes of corn in the process of pathogen infection is of great significance for disease resistance breeding and disease control.

[0003] In recent years, with the development of high-throughput sequencing technology, multi-level omics data such as transcriptome, methylation, and proteome have been widely used in crop and pathogen interaction research. Among them, transcriptome sequencing (RNA-seq) provides rich information for studying gene expression regulation, but its data dimension is extremely high, and the sample size is limited, which leads to the problems of overfitting and ignoring nonlinear relationships of traditional statistical methods and machine learning models.

[0004] Existing key gene screening methods mainly include differential expression analysis, weighted co-expression network analysis (WGCNA), and machine learning models based on feature selection (such as LASSO regression, support vector machine recursive feature elimination SVM-RFE, etc.). However, these methods still have the following shortcomings: (1) strong linear assumption of features: difficult to capture complex nonlinear regulation patterns; (2) ignoring the structure information of gene interaction: not considering the topological constraints of biological networks, the stability of screening results is poor; (3) lack of interpretability: although traditional deep learning models have strong feature learning ability, it is difficult to trace the importance of specific genes.

[0005] In recent years, deep embedding learning methods such as graph variational autoencoder (Variational Autoencoder, VAE) have shown advantages in multi-omics feature dimension reduction and phenotype prediction, which can effectively alleviate the "dimension disaster" problem. However, ordinary VAE models cannot incorporate biological network priors, and it is also difficult to achieve balanced optimization between tasks in multi-task prediction. SUMMARY

[0006] In view of the defects and deficiencies existing in the prior art, the purpose of the present application is to provide a key pathogenic factor screening method based on deep embedding and gene network constraint, which not only integrates gene network prior knowledge in the model structure, significantly improves the biological rationality and stability of the screening result, but also realizes nonlinear dimension reduction of high-dimensional transcriptome data and key gene explainability screening, and provides a new computing framework for crop disease molecular mechanism research.

[0007] To achieve the above-mentioned purpose, the technical scheme adopted by the present application is as follows:

[0008] The key pathogenic factor screening method based on deep embedding and gene network constraint comprises the following steps:

[0009] Step 1. Obtain the gene expression data of crops under the condition of pathogen infection, and obtain the gene expression matrix X after preprocessing;

[0010] Step 2. Construct a gene interaction network corresponding to the gene expression matrix of step 1 based on a public protein interaction database, generate a sparse adjacency matrix A, and form a sparse mask matrix M by adding a unit matrix I;

[0011] Step 3. Construct a graph variational autoencoder model with gene network constraint, the input of the graph variational autoencoder model is the preprocessed gene expression matrix X, the sparse mask matrix M is used in the form of element-level multiplication to limit the weight connection of the encoder, the latent structure features are extracted, the decoder is used to reconstruct the input gene expression matrix according to the latent structure features, and the multi-task prediction module at the output end predicts the infection stage probability distribution and pathogen load through parallel classification branch and regression branch;

[0012] Step 4. Train the constructed model:

[0013] First, minimize the reconstruction loss and regularization term of the graph variational autoencoder model; then train only the multi-task prediction module on the basis of freezing the parameters of the encoder and the decoder; finally, unfreeze all the model parameters and perform end-to-end joint optimization;

[0014] Step 5. Perform feature attribution analysis on the input genes according to the task output of the graph variational autoencoder model, and calculate the importance score of each gene;

[0015] Step 6. Generate a sorted key gene candidate set according to the importance score.

[0016] As a preferred embodiment of the present application, the preprocessing step in step 1 is: performing differential expression analysis on the gene expression data, screening differential expression genes meeting the threshold condition, and then performing variance stabilization transformation and standardization processing to construct the gene expression matrix X.

[0017] As a preferred embodiment of the present invention, the gene interaction network described in step 2 is constructed based on the STRING database, and a sparse adjacency matrix is ​​generated using interaction relationships with a confidence level ≥ 700.

[0018] As a preferred embodiment of the present invention, the expression for the encoder in step 3 is:

[0019] H (1) =ReLU((M⊙XW (1) )+b (1) );

[0020] μ=(M⊙H (1) W μ )+b μ ;

[0021] logσ 2 =(M⊙H (1) W σ )+b σ ;

[0022] Where X is the input gene expression matrix; W (1) W μ W σ b are trainable weight parameters; (1) ,b μ b σ σ represents the bias term; M⊙(•) represents the sparse mask matrix constraint; ReLU(•) is the nonlinear activation function, μ represents the mean of the latent distribution, and σ represents the standard deviation of the latent distribution;

[0023] Potential structural characteristics z=μ+σ⊙ε, ε—N(0,I).

[0024] As a preferred embodiment of the present invention, the expression for the decoder in step 3 is:

[0025] ;

[0026] in, W represents the reconstruction matrix output by the decoder. (dec) b are the trainable weight parameters for the decoding layer. (dec) For trainable bias terms of the decoding layer, A symmetric normalized sparse adjacency matrix is ​​defined as: A is a sparse adjacency matrix. identity matrix It is a degree matrix.

[0027] As a preferred embodiment of the present invention, when training the constructed model in step 4, the expression for minimizing the reconstruction loss and regularization term of the graph variational autoencoder model is as follows:

[0028] L vae =L rec +βD kl ;

[0029] Among them, L vae L represents the overall loss function of the graph variational autoencoder. rec D is the reconstruction error of the input and reconstruction matrix. kl Here, β is the potential distribution regularization term, and β is the loss balance coefficient.

[0030] The optimization objective for training the multi-task prediction module is:

[0031] L task = λ1L class + λ2L reg ;

[0032] Among them, L task L represents the overall loss function of the multi-task learning module. class For cross-entropy loss in infection phase classification, L reg λ1 and λ2 are the mean square error of the pathogen load regression, and the loss balance coefficients.

[0033] Minimize the overall loss function for end-to-end joint optimization:

[0034] L total =L rec +λ1L class +λ2L reg +βD Kl ;

[0035] Among them, L total This is the total loss function.

[0036] As a preferred embodiment of the present invention, step 5 employs a multi-dimensional interpretability scoring module for feature attribution analysis. This multi-dimensional interpretability scoring module includes input gradient method, ensemble gradient method, sparse weight norm method, and stability selection method. The gradient sensitivity score of the gene is calculated based on the input gradient method. Calculate the ensemble gradient score of genes based on the ensemble gradient method. The weight norm score of genes is calculated based on the sparse weight norm method. The stability frequency Freq of genes was calculated based on the stability selection method. j Finally, the overall score is calculated.

[0037] As a further preferred embodiment of the present invention, the expression for the comprehensive score is: .

[0038] Advantages and beneficial effects of the present invention:

[0039] (1) This invention proposes a graph variational autoencoder (Graph-VAE) model that integrates gene interaction network constraints. By introducing prior knowledge of biological networks to construct a sparse connection structure, it jointly completes multi-task learning in a low-dimensional embedding space. Furthermore, it combines various interpretability analysis strategies (including gradient attribution, integrated gradient and stability selection) to enable it to efficiently screen pathogenic factor genes that play a key regulatory role in the interaction between maize and Fusarium graminearum while ensuring prediction performance.

[0040] (2) The key pathogenic factor screening method proposed in this invention integrates prior information from gene networks to enhance the biological rationality of feature extraction: it introduces maize gene interaction networks from public databases such as STRING, and constructs a sparse adjacency matrix to achieve graph structure constraints between genes, enabling the model to retain the regulatory relationships between genes during the feature extraction stage. Compared with traditional fully connected neural networks, this method can effectively avoid the noise propagation problem caused by random connections, making the obtained potential features more consistent with the real biological network structure.

[0041] (3) Employing a graph variational autoencoder structure to achieve nonlinear feature compression and latent space reconstruction: This invention uses a graph convolutional variational encoder (Graph-VAE) to reduce the dimensionality of high-dimensional gene expression data and model latent features, which can adaptively capture complex nonlinear gene regulation patterns and achieve accurate characterization of dynamic changes in gene expression. Compared with traditional PCA or autoencoder methods, this model has stronger nonlinear fitting ability and higher noise robustness.

[0042] (4) The multi-task learning framework enhances the model's ability to model pathogen infection dynamics: This invention designs a joint classification and regression task at the decoder output. The classification task is used to distinguish infection stages (e.g., 0h, 12h, 24h), and the regression task is used to predict pathogen load, i.e., the relative infection intensity of Fusarium graminearum in the host. This multi-task mechanism enables the model to learn the representation of "whether infected" and "infection degree" simultaneously, thereby obtaining a more accurate description of the disease process; it solves the technical problem that existing gene screening methods have learning conflicts in different tasks (e.g., time point classification and pathogen load prediction), making it difficult to achieve multi-objective balance optimization.

[0043] (5) Achieving model interpretability and automatically screening key pathogenic factor genes: Based on the output of the Graph-VAE decoder, this invention combines a feature attribution mechanism (including gradient sensitivity analysis and graph convolution weight contribution calculation) to calculate the influence intensity of each gene on the infection phenotype prediction results, thereby obtaining gene importance scores; by ranking importance and combining stability screening, key pathogenic genes that play a dominant role in the infection process can be automatically identified, thus achieving interpretability of deep learning results; solving the technical problem that existing deep learning models do not consider the biological interaction between genes, and the model structure is inconsistent with the biological network, resulting in a lack of biological interpretability of the screening results.

[0044] (6) Robust and reliable results with cross-omics scalability: This invention introduces Bootstrap resampling and stability statistics methods to comprehensively analyze the results of multiple model trainings, thereby reducing fluctuations caused by random initialization and improving the stability and reproducibility of screening results; it solves the technical problem that existing gene screening methods lack a stable and interpretable analysis mechanism, making it difficult to accurately locate key regulatory genes. In addition, this framework can be extended to multi-omics data such as transcriptomics, proteomics, and metabolomics to achieve cross-level identification of disease regulatory factors.

[0045] (7) It has broad application value and promotion potential: The method proposed in this invention can be applied to the screening of key pathogenic genes in the maize-Fusarium graminearum interaction system, and can also be extended to the study of disease resistance mechanisms in other crop-pathogen interaction systems. This method provides a new and efficient computational tool for crop disease control, disease-resistant variety breeding and molecular breeding, and has high scientific research and industrial application value. Attached Figure Description

[0046] Figure 1 This is a flowchart of the key pathogenic factor screening method based on deep embedding and gene network constraints of the present invention;

[0047] Figure 2 The graph-VAE model constructed for this invention;

[0048] Figure 3 This refers to the classification performance metrics of the Graph-VAE model in this invention.

[0049] Figure 4 Receiver operating characteristic (ROC) curve of the Graph-VAE model in a classification task during the period of maize pathogen infection;

[0050] Figure 5 The results of gene ontology GO and KEGG pathway enrichment analysis were performed on the set of 50 key genes ranked by importance in the model output. Detailed Implementation

[0051] To enable those skilled in the art to better understand the technical solutions and advantages of the present invention, the present application will be described in detail below with reference to the accompanying drawings, but this is not intended to limit the scope of protection of the present invention.

[0052] like Figure 1 , Figure 2 As shown, this invention provides a method for screening key pathogenic factors based on deep embedding and gene network constraints, the method comprising the following steps:

[0053] Step 1. Data Collection and Preprocessing:

[0054] Gene expression data (transcriptome expression data) of maize before and after Fusarium graminearum infection were obtained from the GEO database, number GSE147508. The original data files were downloaded in the format of count matrix (gene × sample) and sample information file. Nine wild-type B73 maize samples were selected from them, which were three groups of samples infected by Fusarium graminearum at 0 hours, 12 hours and 24 hours.

[0055] Screening for highly variable genes: Differential expression analysis was performed on samples from different infection stages to screen for differentially expressed genes that met the threshold conditions, and high-variable genes were retained by variance screening; then variance stabilization transformation (log2 transformation) and Z-score standardization were performed to eliminate batch effects and ensure data comparability.

[0056] Specifically, the expression variance of each gene in the nine samples is calculated, and the top 1000 genes with the highest variance are selected as model inputs. Based on this, an input gene expression matrix X∈R is constructed. nxp Where n is the number of samples and p is the number of genes; this matrix X serves as the input feature data for the subsequent Graph-VAE model.

[0057] Define the labels by defining two labels based on the sample metadata file:

[0058] The classification labels (infection period) are F0, F12, and F24, which are three sampling time points after Fusarium graminearum infected corn, corresponding to pre-infection, early infection, and mid-to-late infection; in the model, they are coded as 0, 1, and 2, respectively.

[0059] Pathogen load: To reflect the dynamic process of pathogen accumulation in maize tissues over time, this invention introduces a continuous indicator—pathogen load—to characterize the degree of pathogen infection. Since this dataset does not provide absolute abundance measurements of pathogens, the pathogen load is approximated based on infection time.

[0060] .

[0061] This assignment does not represent three discrete categories, but rather represents continuous sampled values ​​of pathogen load at different time points, i.e., the numerical approximation of pathogen infection intensity on the time axis.

[0062] During training, the model treats this variable as a continuous quantity for regression optimization, learning its smooth change pattern in the latent space, thereby fitting the nonlinear growth trend of pathogen infection degree over time.

[0063] The dataset was divided into training set (70%), validation set (15%), and test set (15%) according to sample number. Stratified sampling was used to ensure consistency of time distribution.

[0064] After the above steps, a standardized gene expression matrix file is generated:

[0065] GSE174508_expr_ZmID_final.csv: Rows correspond to genes, columns correspond to samples; GSE174508_sample_label.csv: Includes infection period and pathogen_load.

[0066] Step 2. Gene Network Construction:

[0067] To introduce biological network constraints into the model, this invention utilizes the STRING v11.5 database (protein interaction database) to construct a maize gene interaction network. Specifically, using the input genes selected above as nodes, interactions with a confidence level ≥700 in the STRING database are extracted to generate a sparse adjacency matrix A∈{0,1}. pxp Among them, A i,j =1 indicates that gene g i With g j Interactions exist. To ensure the model structure conforms to the topological characteristics of biological networks, an identity matrix I is added to the sparse adjacency matrix to form a sparse mask matrix M = A + I; where the identity matrix I is defined as:

[0068]

[0069] This design adds self-connections to each node in the sparse adjacency matrix, thereby preserving the gene's own expression characteristics during graph convolution propagation and avoiding the problem of over-smoothing of features caused by relying solely on the features of neighboring nodes.

[0070] The sparse adjacency matrix A generated based on the STRING database is saved in MatrixMarket format (filename: GSE174508_adjacency_from_STRING.mtx) and the corresponding gene sequence text (GSE174508_gene_list.csv) to guide the Graph-VAE model in implementing a sparse connection structure with gene network constraints in the coding layer.

[0071] Step 3. Construct a graph variational autoencoder (Graph-VAE) model constrained by a gene network:

[0072] This invention uses GSE17450 maize transcriptome data from the GEO database as an example and combines STRING gene interaction network information to construct a gene network-constrained graph variational autoencoder (Graph-VAE) model.

[0073] like Figure 2 As shown, the Graph-VAE model proposed in this invention mainly consists of three parts: an encoder, a decoder, and a multi-task prediction module. The Graph-VAE model also includes an input layer, where the input is a standardized gene expression matrix X∈R. nxp , where n is the number of samples and p is the number of genes;

[0074] In the encoder (encoder layer) computation process of Graph-VAE, this invention employs a sparse mask matrix M to restrict weight connections in the form of element-wise multiplication (the input layer is multiplied by the gene interaction sparse mask matrix M to constrain the connection relationships between neurons, allowing only gene pairs with biological interactions to transmit information), thereby achieving structured feature learning constrained by the gene network and sparse graph structure propagation. The encoding layer calculation formula is as follows:

[0075] H (1) =ReLU((M⊙XW (1) )+b (1) )

[0076] μ=(M⊙H (1) W μ )+b μ

[0077] logσ 2 =(M⊙H (1) W σ )+b σ

[0078] Where, X∈R nxp Input gene expression matrix; W (1) W μ Wσ b are trainable weight parameters; (1) ,b μ b σ σ is the bias term; M⊙(•) is the sparse mask matrix constraint, which retains only the connections that exist in the interaction network; ReLU(•) is the nonlinear activation function, μ represents the mean of the latent distribution, and σ represents the variance of the latent distribution.

[0079] The model obtains a latent variable z through reparameterized sampling, which is used to describe the potential structural features of maize gene expression:

[0080] z = μ + σ⊙ε, ε — N(0, I)

[0081] This design ensures that the feature learning process is consistent with the biological network topology, enabling structured latent variable modeling under gene interaction constraints.

[0082] The decoder (Dncoder) part of the Graph-VAE in this invention is used to reconstruct the input gene expression matrix based on the latent variable z, realizing the inverse mapping of the feature space. The decoder also adopts a two-layer graph convolution structure, and its calculation expression is as follows:

[0083]

[0084] in, W represents the reconstruction matrix output by the decoder. (dec) b are the trainable weight parameters for the decoding layer. (dec) For trainable bias terms of the decoding layer, A symmetric normalized sparse adjacency matrix is ​​defined as: A is a sparse adjacency matrix. identity matrix For a degree matrix, its diagonal elements are defined as:

[0085] in, Representative gene g i With g j The degree matrix, Representative gene g i With g j The sparse adjacency matrix, Representative gene g i With g j The identity matrix;

[0086] Through the above normalized adjacency operation, the decoder can balance the feature propagation weights of high-connectivity nodes and low-connectivity nodes during the reconstruction process, which not only ensures numerical stability but also allows the model to retain the structural consistency of biological networks in the latent space.

[0087] Finally, the reconstruction matrix output by the decoder The reconstruction error L between the input gene expression matrix X and the input gene expression matrix X rec Defined as:

[0088]

[0089] This loss term measures the model's ability to fit the original gene expression distribution and is optimized together with the classification loss, regression loss, and regularization term during training.

[0090] Furthermore, in order to simultaneously identify the infection stage of maize and quantitatively predict the pathogen load, this invention designs a multi-task prediction module at the decoder output, and adopts a multi-task learning structure that combines classification and regression optimization, so that the model can learn discrete and continuous biological features in the latent space at the same time, thereby improving the biological rationality and interpretability of the latent representation.

[0091] In this embodiment, the input of the multi-task prediction module is the reconstructed gene expression matrix output by the decoder, which contains the feature expression information of the sample after reconstruction in the latent space.

[0092] Based on this, the Graph-VAE model sets up two parallel branches:

[0093] (1) Classification branch: used to identify the infection stage of the sample, and output the probability distribution of the three stages: before infection (F0), 12 hours after infection (F12), and 24 hours after infection (F24). The calculation formula is as follows:

[0094]

[0095] Among them, W c b c These are the classification layer parameters in the multi-task prediction module. ∈R nx3 This represents the prediction result in the probability phase.

[0096] (2) Regression branch: Used to predict pathogen load, i.e., the quantitative influence of pathogens at different infection stages. Pathogen load, as a continuous indicator of infection intensity, can be approximately set to 0, 1, and 2 based on the infection time, corresponding to the three stages F0, F12, and F24. The regression branch adopts a linear mapping structure:

[0097]

[0098] Among them, W r b r These are the regression layer parameters in the multi-task prediction module. ∈R nx1 This represents the predicted pathogen load.

[0099] To optimize both classification and regression tasks simultaneously, this invention defines the total loss function (comprehensive loss function) as follows:

[0100] L total =L rec +λ1L class +λ2L reg +βD Kl

[0101] in, The mean square error between the input and the reconstructed matrix;

[0102] : Cross-entropy loss for classifying infection phases; where, This represents the true category label of the i-th sample, using one-hot encoding, to indicate the infection stage (0 h, 12 h, or 24 h) of the sample. This represents the predicted probability distribution of the infection stage for the i-th sample by the model, and is the output of Softmax.

[0103] : represents the mean squared error of pathogen load prediction; where, The true pathogen load value of the representative sample; This represents the pathogen load value predicted by the model.

[0104] D Kl λ1, λ2, and β are the potential distribution regularization terms; λ1, λ2, and β are the loss balance coefficients.

[0105] Through this multi-task joint optimization structure, the model can simultaneously learn the features of infection stages (discrete changes) and pathogen load (continuous changes) based on shared latent representations, thereby capturing the dynamic changes in maize during Fusarium graminearum infection. The classification task provides discrete label constraints, while the regression task provides continuous trend information. The two tasks are mutually regularized, effectively enhancing the model's generalization ability and improving the biological interpretability of latent features.

[0106] Step 4. Train the constructed model.

[0107] To ensure the stability of the model in latent space learning, task optimization, and overall convergence, this invention adopts a three-stage training strategy, which includes a Graph-VAE pre-training stage, a task training stage, and a joint fine-tuning stage, to progressively optimize model performance and improve the biological interpretability of the results. The specific steps are as follows:

[0108] Step 4.1. Graph-VAE pre-training:

[0109] This stage aims to minimize the reconstruction loss and regularization term of the graph variational autoencoder model, enabling the model to learn a stable and structured latent space representation. The training objective is:

[0110] L vae =L rec +βD kl

[0111] Among them, L vae L represents the overall loss function of a graph variational autoencoder (Graph-VAE). rec D is the mean square error (reconstruction error) of the input and reconstruction matrix. kl Let β be the Kullback-Leibler divergence between the latent distribution and the standard normal distribution, i.e., the regularization term of the latent distribution, and β be the annealing weight, also known as the loss balance coefficient. In the early stages of training, a Kl annealing strategy is used, gradually increasing β from 0 to β... max (typically β) max =0.1), to prevent the occurrence of latent variable collapse.

[0112] The optimizer uses AdamW (learning rate 1×10⁻⁶). -3 Weight decay 1×10 -4 The training rounds are 50–100, and the early stopping strategy uses the validation set reconstruction error as the monitoring metric (patience value 15 rounds).

[0113] Through this pre-training stage, the model can learn the potential feature distribution that reflects the structure of the gene network, providing a stable initial representation for downstream tasks.

[0114] Step 4.2. Task Training:

[0115] In this stage, based on the frozen encoder and decoder parameters, only the classification and regression task heads (multi-task prediction module) are trained to learn the mapping relationship between latent features and infection stage and pathogen load.

[0116] The optimization objective is:

[0117] L task = λ1L class + λ2L reg

[0118] Among them, L task L represents the comprehensive loss function of the multi-task learning module, used to simultaneously optimize the infection stage classification task and the pathogen burden regression task. class For cross-entropy loss in infection phase classification, L reg λ1 and λ2 represent the mean squared error of the pathogen burden regression, and λ1 and λ2 are the task weight coefficients, also known as loss balance coefficients. The optimizer uses Adam (learning rate 5×10⁻⁶).-4 The training run can be 20–50 epochs. By fine-tuning the task head on a stable latent space, the accuracy and generalization ability of the model on classification and regression tasks can be effectively improved.

[0119] Step 4.3. Joint fine-tuning

[0120] To further improve overall performance and multi-task collaboration, all model parameters are unfrozen in this stage, and end-to-end joint optimization is performed to minimize the comprehensive loss function:

[0121] L total =L rec +λ1L class +λ2L reg +βD Kl

[0122] The optimizer still uses AdamW (learning rate 1×10⁻⁶). -4 The training run consists of 50–150 epochs, with an early stopping strategy based on the comprehensive loss from the validation set. During training, the task loss weights are dynamically balanced to ensure the model maintains reconstruction accuracy while also considering performance in predicting infection stages and pathogen burden. Through this joint training phase, the model achieves synergistic optimization of latent features, reconstructed representations, and task outputs, resulting in optimal overall performance and biological interpretability.

[0123] Step 5. Based on the constructed multi-dimensional interpretability scoring module, a comprehensive evaluation of gene importance is conducted to form a unified gene importance score;

[0124] After the model training is completed and a stable latent representation is obtained, in order to identify genes that play a key regulatory role in the pathogen infection process, this invention designs a multi-dimensional interpretability scoring module. This multi-dimensional interpretability scoring module includes four scoring methods, which comprehensively evaluate gene importance through gradient contribution, weight contribution, and stability selection, ultimately forming a unified gene importance score. The specific steps are as follows:

[0125] Step 5.1. Input gradient method:

[0126] The input gradient method is used to measure the sensitivity of the model output to changes in the expression of each gene, reflecting the direct contribution of genes in the prediction task.

[0127] For the classification output y of the sample c (e.g., infection class probability in the Softmax output) and input genes Its importance is defined as: ;

[0128] Among them, S j For the first Gradient sensitivity score of each gene; x j Input gene expression matrix The Middle Expression values ​​of individual genes; y c This represents the probability value corresponding to category c in the model's predicted output. It is the gradient of the output with respect to the input gene.

[0129] The gradient is taken as an absolute value to eliminate the influence of sign, so that the score reflects the magnitude of the contribution rather than the direction.

[0130] Step 5.2. Integrated Gradient Method:

[0131] To address the gradient saturation problem that may occur with single-point gradient methods, this invention employs integrated gradients to further evaluate the path contribution of genes to the predicted output. For gene x... j Its integrated gradient score is defined as:

[0132]

[0133] Among them, IG j The integrated gradient score for the j-th gene; This represents the true expression value of the j-th gene in the input sample. is the reference expression value (e.g., zero vector or average expression) of the j-th gene in the baseline sample; F(•) represents the model prediction function (which can be a classification probability or regression output); a is the integral path parameter from 0 to 1; ) is from the baseline vector to the real input vector The interpolation vector.

[0134] Integrated gradients contribute more robustly and smoothly than single-point gradients by accumulating contributions through the integral path.

[0135] Step 5.3. Sparse Weight Norm Method:

[0136] To measure the "structural importance" of each gene in the model structure, this invention characterizes the degree of gene influence by calculating the weight norm of the task-related layer.

[0137] I j = ||W :,j (1) ||2

[0138] Among them, I j W represents the weight norm score of the j-th gene; :,j (1) For the weight matrix and genes The corresponding column (all weights mapping from latent features to output); ||.||2 is the L2 norm (Euclidean mode).

[0139] If the model has sparse constraints (such as L1 / L0), the metric will naturally exhibit sparsity, improving the interpretability of importance screening.

[0140] Step 5.4. Stability Selection Method:

[0141] To improve the robustness of the screening results, this invention employs a stability selection method to mitigate fluctuations caused by data noise and random model initialization. The specific steps are as follows:

[0142] Step 5.4.1. Repeat the training of the model N times (usually N=100), each time constructing a training set through different random initializations or bootstrap resampling;

[0143] Step 5.4.2. Calculate the gradient score, ensemble gradient score, and weight norm score in each training iteration;

[0144] Step 5.4.3. Record the set of genes that rank the Top-K (e.g., K=100) in each training session;

[0145] Step 5.4.4. Calculate the frequency of each gene appearing in the Top-K list during N repeated training iterations:

[0146] Freq j

[0147] Among them, Freq j ∈[0,1] represents the stability frequency of gene j;

[0148] Step 5.5. Overall Score and Key Gene Screening:

[0149] Calculate the overall score based on the scores from the four scoring methods mentioned above:

[0150]

[0151] Among them, Score j The overall score is calculated (the higher the score, the more important it is).

[0152] In this invention, the combination of three types of contribution and stability frequency effectively improves the reliability of the screening and avoids the influence of single indicator bias on the results.

[0153] Step 6. Generate a high-confidence key gene candidate set based on the comprehensive score;

[0154] In this invention, a comprehensive score is selected. j The top M genes (e.g., M=50 or M=100) are considered as a candidate set of key regulatory genes in maize when it is infected by Fusarium graminearum.

[0155] Step 7. Result Verification:

[0156] This invention performs GO and KEGG enrichment analyses on candidate genes to verify their enrichment in pathways such as defense response, hormone signaling, and cell wall synthesis, and outputs a final list of key genes with high confidence.

[0157] In this embodiment, the data source is the GEO database GSE174508, which contains transcriptome data of maize B73 strain infected at 0h, 12h, and 24h, totaling 9 samples.

[0158] Model settings: latent dimension d=16, β=2, Dropout=0.3, learning rate 1e-3, optimizer AdamW.

[0159] Training phase: Graph-VAE pre-training for 100 rounds; downstream task training for 100 rounds; joint fine-tuning for 150 rounds.

[0160] Results: Five-fold cross-validation was used to evaluate the performance of the model on the classification task of infection period (0h, 12h, 24h). The results are as follows: Figure 3 As shown in the figure, the Graph-VAE model provided by this invention has good overall classification performance and strong predictive ability.

[0161] To verify the effectiveness of the model of this invention in the classification task of maize pathogen infection stage, the receiver operating characteristic (ROC) curve and area under the curve (AUC) were calculated for the test set samples. The results are shown in Figure 4. The ROC curve of the Graph-VAE model of this invention generally shows a shape close to the upper left corner, and the AUC reaches 0.80, indicating that the model maintains high sensitivity and specificity at different thresholds.

[0162] To verify the biological relevance of key maize genes screened by the deep learning model in the Fusarium graminearum infection process, gene ontology GO and KEGG pathway enrichment analyses were performed on the set of 50 key genes ranked by importance output by the model. The results are as follows: Figure 5 As shown. Analysis results indicate that the key genes screened were significantly enriched in multiple functional categories closely related to plant-pathogen interactions (corrected p-value <0.05), including:

[0163] (1) Biological Process: microtubule depolymerization, regulation of DNA-templated transcription, mitotic cell cycle, and diterpenoid metabolic process.

[0164] (2) Cellular component: microtubule;

[0165] (3) Molecular Function: heme binding and structural constituent of cytoskeleton.

[0166] (4) KEGG pathway: motor proteins.

[0167] The above enrichment results have clear biological significance for plant disease resistance. Among them, diterpenoid metabolism is the core pathway for maize to synthesize antifungal phytoprotectants; microtubule dynamic regulation and cytoskeleton remodeling are involved in limiting the spread of pathogens in tissues and mediating the intracellular transport of defense-related substances; the enrichment of transcriptional regulation and cell cycle-related genes reflects the immune response and resource redistribution mechanism initiated by the host under pathogen infection; and heme-binding proteins (such as the cytochrome P450 family) are widely involved in defense responses such as phytoprotectant modification and reactive oxygen species bursts.

[0168] In summary, the key genes identified in this invention are significantly enriched in biological pathways known to be highly associated with Fusarium graminearum resistance, which functionally verifies the important role of these genes in maize's response to Fusarium graminearum infection, fully supporting the practicality and technical effectiveness of the present invention.

[0169] The above-described specific examples are for illustrative purposes only and are not intended to limit the scope of the invention. Those skilled in the art can make various simple deductions, modifications, or substitutions based on the principles of this invention. Therefore, the scope of protection of this invention should be determined by the scope of the claims.

Claims

1. A method for screening key pathogenic factors based on deep embedding and gene network constraints, characterized in that, The method includes the following steps: Step 1. Obtain gene expression data of crops under pathogen infection conditions, and obtain gene expression matrix X after preprocessing; Step 2. Construct a gene interaction network corresponding to the gene expression matrix in Step 1 based on a public protein interaction database, generate a sparse adjacency matrix A, and form a sparse mask matrix M by adding it to the identity matrix I; Step 3. Construct a graph variational autoencoder model constrained by a gene network. The input of the graph variational autoencoder model is the preprocessed gene expression matrix X. The encoder uses a sparse mask matrix M to restrict the weight connections in the form of element-wise multiplication and extract latent structural features. The decoder is used to reconstruct the input gene expression matrix based on the latent structural features. The multi-task prediction module at the output end predicts the probability distribution of the infection stage and the pathogen load through parallel classification and regression branches. Step 4. Train the constructed model: First, minimize the reconstruction loss and regularization term of the graph variational autoencoder model; then, train only the multi-task prediction module while freezing the encoder and decoder parameters; finally, unfreeze all model parameters and perform end-to-end joint optimization. Step 5. Perform feature attribution analysis on the input genes based on the task output of the graph variational autoencoder model, and calculate the importance score of each gene; Step 6. Generate and sort a candidate set of key genes based on importance scores; The encoder expression in step 3 is: A (1) =ReLU((M⊙XW (1) )+b (1) ); μ=(M⊙H (1) W μ )+b μ ; logσ 2 =(M⊙H (1) W σ )+b σ ; Where X is the input gene expression matrix; W (1) W μ W σ b are trainable weight parameters; (1) ,b μ b σ σ represents the bias term; M⊙(⋅) is the sparse mask matrix constraint; ReLU(⋅) is the nonlinear activation function, μ represents the mean of the latent distribution, and σ represents the standard deviation of the latent distribution; The potential structural feature is z = μ + σ⊙ε, ε ∼ N(0, I); The expression for the decoder in step 3 is: ； in, W represents the reconstruction matrix output by the decoder. (dec) b are the trainable weight parameters for the decoding layer. (dec) For trainable bias terms of the decoding layer, A symmetric normalized sparse adjacency matrix is ​​defined as: ; A It is a sparse adjacency matrix. identity matrix It is a degree matrix.

2. The method for screening key pathogenic factors based on deep embedding and gene network constraints according to claim 1, characterized in that, The preprocessing step in step 1 is as follows: perform differential expression analysis on gene expression data, screen differentially expressed genes that meet the threshold conditions, and then perform variance stabilization transformation and standardization to construct the gene expression matrix X.

3. The method for screening key pathogenic factors based on deep embedding and gene network constraints according to claim 2, characterized in that, The gene interaction network described in step 2 is constructed based on the STRING database, and a sparse adjacency matrix is ​​generated using interaction relationships with a confidence level of ≥700.

4. The method for screening key pathogenic factors based on deep embedding and gene network constraints according to claim 3, characterized in that, When training the constructed model in step 4, the expression for minimizing the reconstruction loss and regularization term of the graph variational autoencoder model is as follows: L vae =L rec +βD kl ; Among them, L vae L represents the overall loss function of the graph variational autoencoder. rec D is the reconstruction error of the input and reconstruction matrix. kl is the potential distribution regularization term, and is the loss balance coefficient; The optimization objective for training the multi-task prediction module is: L task = λ1L class + λ2L reg ; Among them, L task The comprehensive loss function representing the multi-task learning module, 𝐿 class Cross-entropy loss for classification of infection phases, 𝐿 reg λ1 and λ2 are the mean square error of the pathogen load regression, and the loss balance coefficients. Minimize the overall loss function for end-to-end joint optimization: L total =L rec +λ1L class +λ2L reg +βD KL ; Among them, L total This is the total loss function.

5. The method for screening key pathogenic factors based on deep embedding and gene network constraints according to claim 4, characterized in that, Step 5 employs a multi-dimensional interpretability scoring module for feature attribution analysis. This module includes input gradient method, ensemble gradient method, sparse weight norm method, and stability selection method. The gradient sensitivity score of the gene is calculated based on the input gradient method. Calculate the ensemble gradient score of genes based on the ensemble gradient method. The weight norm score of genes is calculated based on the sparse weight norm method. Calculate gene stability frequencies based on stability selection methods Freq j Finally, the overall score is calculated.

6. The method for screening key pathogenic factors based on deep embedding and gene network constraints according to claim 5, characterized in that, The expression for the comprehensive score is: 。

Images