Transcriptome data-based plant stress mitigation key gene tracing method and system

Abstract

The application relates to the technical field of biological information, and discloses a plant stress relief key gene tracing method and system based on transcriptome data. The method comprises the following steps: performing high-throughput transcriptome sequencing on plant leaves of four treatment groups, obtaining clean transcriptome data through quality filtering; obtaining a whole genome expression matrix through transcript assembly and redundancy removal processing, performing intersection screening on stress induction up-regulated genes and relief agent silencing down-regulated genes, and obtaining a stress induction-relief agent silencing candidate gene set; constructing a stress induction co-expression topology network by taking the candidate gene set as a silencing annotation basis, calculating a weighted neighbor silencing rate to obtain a topology feature matrix; training a selective silencing probability prediction model through a logistic regression model by taking the topology feature matrix and the candidate gene set, and performing channel enrichment calculation by taking the silencing probability as a weight to screen and obtain a core tracing gene set. The application improves the biological credibility of key gene tracing results and the cross-scene adaptability of the method.

Description

Technical Field

[0001] This application relates to the field of bioinformatics, and in particular to a method and system for tracing the origins of key genes for plant stress relief based on transcriptome data. Background Technology

[0002] Nitric oxide, as an important gaseous signaling molecule in plants, is widely involved in physiological processes such as plant germination, mitochondrial function regulation, hormone signal transduction, and stress response. Exogenous application of the nitric oxide donor sodium nitroprusside can effectively alleviate the toxicity caused by cadmium stress in plants, including alleviating symptoms such as chlorosis, wilting, photosynthetic inhibition, and cell death. Cadmium, as a highly toxic heavy metal, harms normal plant physiological functions through multiple pathways, including inducing changes in membrane composition, disrupting the photosystem electron transport chain, affecting stomatal closure, and inhibiting nutrient absorption. Kentucky bluegrass, due to its well-developed root system and strong cadmium absorption and accumulation capacity, possesses potential phytoremediation capabilities. Transcriptome analysis platforms based on high-throughput RNA sequencing technology provide genome-wide gene expression data support for studying the above-mentioned mitigation mechanisms.

[0003] However, existing methods for tracing key genes in plant stress relief based on transcriptome data have significant shortcomings. These methods typically only compare the stress group and the control group, and the screening of differentially expressed genes is based solely on the statistical test values ​​and fold changes of single pairwise comparisons. This fails to capture gene expression information from both the stress-induced signal and the silencing signal of the alleviating agent, resulting in a large number of logically contradictory genes mixed in with the selected candidate genes. Functional enrichment analysis assigns equal weight to all differentially expressed genes, causing numerous but functionally low metabolic pathways to dominate the enrichment results. The final candidate gene list is merely a post-hoc statistical description of observed expression changes and cannot reveal the intrinsic patterns of the selective silencing of specific gene subsets by alleviating agents. Summary of the Invention

[0004] This application provides a method and system for tracing key genes for plant stress relief based on transcriptome data. It solves the problems in the prior art where transcriptome tracing methods cannot quantitatively model the selective silencing patterns of alleviating agents and the candidate gene screening is limited to post-hoc statistical description. This improves the biological reliability of key gene tracing results and the cross-scenario adaptability of the method.

[0005] Firstly, this application provides a method for tracing the origins of key genes for plant stress relief based on transcriptome data, the method comprising: Step S1: High-throughput transcriptome sequencing was performed on the plant leaves of the control group, stress treatment group, alleviator treatment group and combined treatment group. After quality filtering, clean transcriptome data of each treatment group were obtained. Step S2: The clean transcriptome data of each treatment group are processed by transcript assembly and redundancy removal to obtain a whole genome expression matrix; the intersection of stress-induced upregulated genes and reliever-silenced downregulated genes is screened to obtain a stress-induced-relieving agent silencing candidate gene set; Step S3: Using the stress-induced / alleviating agent silencing candidate gene set as the silencing labeling basis, construct a co-expression topology network for the stress-induced upregulated genes to obtain a stress-induced co-expression topology network; calculate the weighted neighbor silencing rate for each node gene in the network to obtain a topology feature matrix, where the weighted neighbor silencing rate is the weighted proportion of genes belonging to the candidate gene set among the direct neighbors of the node gene to all direct neighbors; Step S4: Train the topological feature matrix and the candidate gene set using a logistic regression model to obtain a selective silencing probability prediction model; calculate the silencing probability score matrix by using the prediction model to obtain the stress-induced upregulated genes, and screen to obtain a high-confidence silencing prediction gene set; calculate the high-confidence silencing prediction gene set by using the silencing probability as a weight through pathway enrichment to obtain a priority pathway ranking list for remission agent selective silencing, and screen to obtain a core source gene set.

[0006] Secondly, this application provides a system for tracing key genes for plant stress relief based on transcriptome data, the system comprising: The sequencing module is used to perform high-throughput transcriptome sequencing on plant leaves from the control group, stress treatment group, alleviator treatment group, and combined treatment group. After quality filtering, clean transcriptome data for each treatment group are obtained. The screening module is used to assemble and deredundate the clean transcriptome data of each treatment group to obtain a whole genome expression matrix; and to perform intersection screening of stress-induced upregulated genes and reliever-silenced downregulated genes to obtain a stress-induced-relieving agent silencing candidate gene set. The construction module is used to construct a co-expression topology network for the stress-induced upregulated genes based on the stress-induced silencing candidate gene set as the silencing label, thereby obtaining a stress-induced co-expression topology network; and to calculate the weighted neighbor silencing rate for each node gene in the network to obtain a topology feature matrix, wherein the weighted neighbor silencing rate is the weighted proportion of genes belonging to the candidate gene set among the direct neighbors of the node gene to all direct neighbors; The calculation module is used to train the topological feature matrix and the candidate gene set using a logistic regression model to obtain a selective silencing probability prediction model; to extrapolate the stress-induced upregulated genes using the prediction model to obtain a silencing probability score matrix, and to screen to obtain a high-confidence silencing prediction gene set; and to calculate the high-confidence silencing prediction gene set by using the silencing probability as a weight through pathway enrichment to obtain a priority pathway ranking list for remission agent selective silencing, and to screen to obtain a core source gene set.

[0007] The technical solution provided in this application establishes a four-treatment group architecture: a control group, a stress treatment group, a remission agent treatment group, and a combined treatment group. This allows transcriptome sequencing data to simultaneously carry gene expression information in two dimensions: stress-induced signals and remission agent silencing signals. Based on this, the intersection of stress-induced upregulated genes and remission agent-silencing downregulated genes is screened to obtain a stress-induced-remission agent silencing candidate gene set. This dual expression feature constraint is established from the data acquisition stage, eliminating logically contradictory genes introduced by the single-comparison differential expression screening method in existing technologies, thus improving the biological reliability of the candidate gene set. In the stress-induced co-expression topology network construction stage, this invention uses the stress-induced-remission agent silencing candidate gene set as the basis for silencing annotation. A weighted neighbor silencing rate is calculated for each node gene in the network. The weighted neighbor silencing rate directly incorporates the member information of the candidate gene set into the topological feature calculation, quantitatively characterizing the local density of remission agent-silenced genes in the network neighborhood of the node gene. This establishes a quantitative correlation between network topological location and selective silencing behavior, a feature completely absent in existing transcriptome tracing methods.

[0008] The weighted neighbor silencing rate, along with node degree, betweenness centrality, and clustering coefficient, are used as input features of the logistic regression model to train a selective silencing probability prediction model. This model outputs the probability value of each stress-induced upregulated gene being selectively silenced by the alleviating agent. This transforms the post-hoc statistical description of observed expression changes in existing technologies into a mechanistic prediction of selective silencing patterns. The resulting silencing probability score matrix is ​​used as the weight for pathway enrichment calculation, making the priority pathway ranking list for alleviating agent selective silencing dominated by the pathway attribution of genes with high silencing probabilities, rather than by the number of genes. This changes the location of the core source gene set from quantity-driven to probability intensity-driven. Finally, the co-expression network construction threshold and silencing probability admission threshold are linked and corrected using qRT-PCR validation data to form a closed-loop source tracing mechanism. This allows the method of this invention to have parameter adaptive correction capabilities under different plant species, different stress types, and different alleviating agent conditions. Attached Figure Description

[0009] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0010] Figure 1 This is a schematic diagram of an embodiment of the plant stress relief key gene tracing method based on transcriptome data in this application. Figure 2 This is a schematic diagram showing the distribution of the number of differentially expressed genes in each of the four pairs of comparisons in this application embodiment. Detailed Implementation

[0011] This application provides a method and system for tracing the origins of key genes for plant stress relief based on transcriptome data. The terms "first," "second," "third," "fourth," etc. (if present) in the specification, claims, and accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in a sequence other than that illustrated or described herein. Furthermore, the terms "comprising" or "having" and any variations thereof are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0012] For ease of understanding, the specific process of the embodiments of this application is described below. Please refer to [link / reference]. Figure 1 One embodiment of the method for tracing the origins of key genes for plant stress relief based on transcriptome data in this application includes: Step S1: High-throughput transcriptome sequencing was performed on the plant leaves of the control group, stress treatment group, alleviator treatment group and combined treatment group. After quality filtering, clean transcriptome data of each treatment group were obtained. Step S2: The clean transcriptome data of each treatment group are processed by transcript assembly and redundancy removal to obtain the whole genome expression matrix; the intersection of stress-induced upregulated genes and reliever-silenced downregulated genes is screened to obtain the stress-induced-relieving agent silencing candidate gene set; Step S3: Using the stress-induced-relief silencing candidate gene set as the silencing labeling basis, construct a co-expression topology network for stress-induced upregulated genes to obtain a stress-induced co-expression topology network; calculate the weighted neighbor silencing rate for each node gene in the network to obtain the topology feature matrix. The weighted neighbor silencing rate is the weighted proportion of genes belonging to the candidate gene set among the direct neighbors of a node gene to all direct neighbors. Step S4: Train the topological feature matrix and candidate gene set using a logistic regression model to obtain a selective silencing probability prediction model; calculate the silencing probability score matrix by using the prediction model to infer stress-induced upregulated genes, and screen to obtain a high-confidence silencing prediction gene set; calculate the priority pathway ranking list of the selective silencing of alleviators by using the silencing probability as the weight of the high-confidence silencing prediction gene set through pathway enrichment calculation, and screen to obtain the core source gene set.

[0013] It is understood that the executing entity of this application can be a system for tracing key genes for plant stress relief based on transcriptome data, or it can be a terminal or a server; no specific limitation is made here. This application's embodiments use a server as an example for illustration.

[0014] Specifically, in the four-treatment group design, the control group was given normal culture medium, the stress treatment group was given cadmium chloride solution, the relief agent treatment group was given sodium nitroprusside solution, and the combined treatment group was given both cadmium chloride and sodium nitroprusside. The four-treatment group setting allows subsequent differential expression analysis to simultaneously capture gene expression changes in two dimensions: cadmium stress-induced signal and nitric oxide relief silencing signal. In the existing technology, only the stress group and the control group are compared, which cannot simultaneously obtain data in the above two dimensions, and therefore cannot identify genes with dual expression characteristics.

[0015] The stress-induced-relief silencing candidate gene set is obtained by taking the intersection of upregulated genes obtained from comparing the control group with the stress treatment group and downregulated genes obtained from comparing the stress treatment group with the combined treatment group. The biological basis for this intersection operation is that genes that simultaneously meet the two conditions of being upregulated by cadmium stress and significantly downregulated after nitric oxide treatment are the direct molecular targets of nitric oxide selective silencing. Existing differential expression screening methods with a single comparison group output logically contradictory genes and true target genes in the same way. This invention eliminates logically contradictory genes from the source by using dual expression feature constraints.

[0016] The weighted neighbor silencing rate is a newly defined topological feature in this invention. This feature is not defined or applied in existing transcriptome source literature. Its essential difference from the three classic topological features in existing technologies—node degree, betweenness centrality, and clustering coefficient—lies in that these only describe the pure structural position of nodes in the network and do not carry any biological information related to allergen silencing behavior. In contrast, the weighted neighbor silencing rate directly incorporates the membership information of the stress-induced allergen silencing candidate gene set into the topological feature calculation, quantitatively characterizing the local density of nitric oxide-silenced genes in the network neighborhood of a node gene, and establishing a quantitative correlation between network topological position and selective silencing behavior. This correlation is completely absent in the existing technical framework. Experimental verification shows that after introducing the weighted neighbor silencing rate into the logistic regression model, the area under the receiver operating characteristic (AUC) curve on the test set is improved by at least 0.05 compared to the baseline model using only node degree, betweenness centrality, and clustering coefficient, demonstrating that the weighted neighbor silencing rate makes an independent technical contribution to the prediction of selective silencing probability.

[0017] In the closed-loop correction mechanism, the squared Pearson correlation coefficient between the results of real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR) verification and the expression levels of RNA sequencing (RNA-Seq) is used. As the triggering basis for parameter correction, The current parameter combination is considered valid if the value is not lower than 0.80. Tighten the threshold for constructing co-expression network edges and the threshold for the silencing probability admission when the value is between 0.65 and 0.80. When the threshold is below 0.65, the above thresholds are further tightened, and the false discovery rate threshold for differential expression analysis is simultaneously adjusted to 0.01. The basis for setting the above three thresholds is as follows: A correlation of at least 0.80 corresponds to a strong correlation criterion between the transcriptome and the experimental validation results, which is consistent with the generally accepted standard for transcriptome validation experiments in this field. A correlation below 0.65 corresponds to a moderately weak correlation, indicating that the biological reliability of the screening results under the current parameter combination is insufficient, and the screening rigor needs to be improved by tightening the parameters.

[0018] In one specific embodiment, step S1 includes: At least three biological replicates were set up for each of the control group, stress treatment group, alleviator treatment group and combined treatment group. Plant leaves were collected from each treatment group, flash-frozen in liquid nitrogen and then stored at -80℃ to obtain frozen leaf samples from each treatment group. Total RNA was extracted from frozen leaf samples of each treatment group using TRIzol reagent. The integrity of the extracted RNA was tested. RNA samples with an RNA integrity index of not less than 7 were processed for library construction to obtain sequencing libraries for each treatment group. The sequencing libraries of each treatment group were input into a high-throughput sequencing platform for paired-end sequencing to obtain the raw sequencing data of each treatment group. The CutAdapter tool was used to remove adapter sequences and filter low-quality bases from the raw sequencing data of each treatment group. The filtering quality threshold was set to be qualified if the Phred value was not lower than 20, the Q20 base ratio was not lower than 97%, and the Q30 base ratio was not lower than 94%. This yielded clean transcriptome data for each treatment group.

[0019] Specifically, at least three biological replicates were used to ensure the statistical validity of the negative binomial distribution model in the subsequent differential expression analysis using the DESeq software. After collection, plant leaves from each treatment group were flash-frozen in liquid nitrogen to instantly terminate the activity of RNA-degrading enzymes in the cellular ribonucleic acid (RNA), preventing alterations to transcriptome information after collection. The leaves were then stored at -80°C for unified extraction. RNA extraction was performed using TRIzol reagent, a single-phase lysis buffer containing phenol and guanidine isothiocyanate, used to simultaneously extract total RNA from plant tissues. After extraction, the RNA integrity was tested using the RNA integrity index. Only samples with an RNA integrity index of at least 7 were allowed into the library construction process; samples with an index below this threshold were excluded due to unacceptable RNA degradation, ensuring the consistency of RNA integrity used in library construction.

[0020] After library construction, the sequencing libraries of each treatment group were input into a high-throughput sequencing platform for paired-end sequencing. Paired-end sequencing refers to sequencing both ends of the same deoxyribonucleic acid (DNA) fragment separately, which can obtain more complete sequence information compared to single-end sequencing. After sequencing, the raw sequencing data of each treatment group were obtained. The raw sequencing data contained adapter sequences and low-quality bases introduced during library construction. The CutAdapter adapter sequence removal and quality filtering tool was used to remove adapter sequences and filter low-quality bases from the raw sequencing data. The filtering quality threshold was set to a Phred value of not less than 20. The Phred value is the logarithmic transformation value of the base identification error probability. A Phred value of 20 corresponds to a base identification error rate of 1%. After filtering, the data were subjected to quality control indicators. The qualified criteria were that the proportion of Q20 bases was not less than 97% and the proportion of Q30 bases was not less than 94%. Q20 and Q30 refer to the proportion of bases with Phred values ​​of not less than 20 and not less than 30, respectively. Data that met the above criteria were output as clean transcriptome data for each treatment group.

[0021] In one specific embodiment, step S2 includes: The clean transcriptome data from each treatment group were input into the Trinity assembly tool for de novo transcript assembly to obtain candidate transcript sequences. The candidate transcript sequences were then processed for sequence similarity redundancy removal using the CD-HIT tool. The sequence similarity threshold was set to 0.95, and the longest sequence was retained as the representative sequence to obtain single gene sequence files. Single gene sequence files were input into the RSEM tool for expression level estimation to obtain the expected read values ​​of each gene in each treatment group. The whole genome expression matrix was constructed using the expected read values. Based on the DESeq tool, four pairs of differential expression analyses were performed on the whole genome expression matrix. The original test values ​​were corrected by the Benjamini-Hochberg method. With the error detection rate after correction being less than 0.05 as the threshold, the stress-induced upregulated gene sets were obtained by comparing the control group with the stress treatment group, the relief agent response gene sets were compared with the control group with the relief agent treatment group, and the relief agent silencing downregulated gene sets were compared with the stress treatment group and the combined treatment group. Intersection screening was performed on the stress-induced upregulated gene set and the remission agent-silenced downregulated gene set to obtain the stress-induced-remission agent-silenced candidate gene set.

[0022] Specifically, clean transcriptome data from each treatment group were input into the Trinity de novo transcriptome assembly software for de novo transcript assembly. Trinity sequentially ran three modules: the linear sequence extension module Inchworm, the graph construction module Chrysalis, and the transcript parsing module Butterfly. The Inchworm module greedily assembled short reads from the clean data into continuous sequences. The Chrysalis module clustered overlapping continuous sequences to construct a de Bruijn graph. The de Bruijn graph is a directed graph constructed with k-mers as nodes and the overlap relationship between adjacent k-mers as edges. A k-mer refers to a subsequence of length k bases extracted from the sequencing read. The Butterfly module traversed the de Bruijn graph and output candidate transcript sequences. A large number of highly similar redundant sequences exist among the candidate transcript sequences. The CD-HIT sequence clustering and redundancy removal software was used to remove redundancy from the candidate transcript sequences. The sequence similarity threshold was set to 0.95, meaning transcripts with a sequence similarity of at least 0.95 were grouped into the same single gene (Unigene), and the longest sequence was retained as the representative sequence, resulting in a single-gene sequence file. The single-gene sequence file was then input into RSEM, an RNA-Seq expression estimation software based on the expectation-maximization algorithm. RSEM used the single-gene sequence file as a reference, comparing the clean transcriptome data of each treatment group with the reference sequence, counting the number of reads aligned to each single gene, and outputting the expected read value for each single gene in each sample. The expected read value is the estimated number of reads obtained by RSEM after probabilistically assigning multiple aligned reads based on the expectation-maximization algorithm. The expectation-maximization algorithm calculates the probability of multiple aligned reads belonging to each single gene in the expectation step and updates the expression estimate of each single gene in the maximization step, alternating between the two steps until convergence. A whole-genome expression matrix is ​​constructed using the expected read values ​​of all single genes in all samples, with rows representing single genes and columns representing samples.

[0023] The whole-genome expression matrix was input into the DESeq tool for pairwise differential expression analysis of four groups. The pairwise comparisons of the four groups were: control group vs. stress treatment group, control group vs. relief agent treatment group, relief agent treatment group vs. combined treatment group, and stress treatment group vs. combined treatment group. DESeq used a negative binomial distribution model to perform statistical tests on the expression differences of each gene between the two groups. The negative binomial distribution model can handle the overdispersion problem that is common in transcriptome data. Overdispersion means that the expression variation of the same gene among repeated samples exceeds the range that the Poisson distribution can describe. The Poisson distribution assumes that the mean and variance are equal, but the variance of gene expression in transcriptome data is usually greater than the mean. The negative binomial distribution models the additional variation that exceeds the Poisson distribution by introducing a discrete parameter. The original test values ​​were corrected to false discovery rates (FDR) using the Benjamini-Hochberg multiple test correction method (hereinafter referred to as BH correction). BH correction is obtained by sorting all the original p-values ​​of the tests in ascending order, multiplying them by the total number of tests, and dividing by the correction coefficient of each ranking. In multiple test correction, the expected false discovery rate is controlled, and a false discovery rate of less than 0.05 is used as the threshold for significant differential expression. The control group and the stress treatment group output the stress-induced upregulated gene set, the control group and the relief agent treatment group output the relief agent response differential gene set, the relief agent treatment group and the combined treatment group output the relief agent and cadmium combined treatment response differential gene set, and the stress treatment group and the combined treatment group output the relief agent silenced downregulated gene set. The relief agent response differential gene set and the relief agent and cadmium combined treatment response differential gene set are retained as background reference data in the whole genome expression matrix for expression pattern consistency verification during subsequent closed-loop correction. The intersection of the stress-induced upregulated gene set and the alleviator-silenced downregulated gene set is taken. The intersection operation retains genes that appear in both gene sets, resulting in a stress-induced-alleviator-silenced candidate gene set. Genes in this candidate gene set simultaneously meet the two conditions of being upregulated by cadmium stress and significantly downregulated after nitric oxide treatment.

[0024] Figure 2The distribution of upregulated and downregulated genes was shown in pairwise comparisons between the control group and the stress treatment group, the control group and the alleviating agent group, the alleviating agent group and the combined treatment group, and the stress treatment group and the combined treatment group. The results were obtained from differential expression analysis of the whole genome expression matrix using the DESeq tool. Specifically, the control group and the stress treatment group had 7022 upregulated genes and 1033 downregulated genes, the control group and the alleviating agent group had 7463 upregulated genes and 2146 downregulated genes, the alleviating agent group and the combined treatment group had 381 upregulated genes and 2101 downregulated genes, and the stress treatment group and the combined treatment group had 466 upregulated genes and 2118 downregulated genes. The stress-induced upregulated gene set and the alleviating agent silenced downregulated gene set were screened to obtain the stress-induced-alleviating agent silenced candidate gene set.

[0025] In one specific embodiment, step S3 involves constructing a co-expression topology network for stress-induced upregulated genes, using a set of stress-induced silencing candidate genes as the basis for silencing annotation. This includes: The expression quantum matrix of stress-induced upregulated genes in all samples was extracted from the whole genome expression matrix; Based on the expression quantum matrix, Pearson correlation coefficients were calculated for each pair of stress-induced upregulated genes to obtain the inter-gene correlation coefficient matrix. Edge construction was performed based on the correlation coefficient matrix between genes. The edge connection condition was that the absolute value of the correlation coefficient was not lower than the initial threshold of 0.90. The stress-induced upregulated genes were used as nodes, and the gene pairs that met the edge connection condition were used as edges to obtain the stress-induced co-expression topology network.

[0026] Specifically, the expression levels of stress-induced upregulated genes in all samples are extracted from the whole-genome expression matrix to form an expression quantum matrix. The expression quantum matrix lists stress-induced upregulated genes as the main genes and all samples as the target samples. Each value in the matrix represents the expected read value of the corresponding gene in its respective sample. Based on the expression quantum matrix, Pearson correlation coefficients are calculated pairwise for each stress-induced upregulated gene. The Pearson correlation coefficient measures the linear correlation between the expression levels of two genes in all samples. For genes i and j, using their expected read values ​​in all samples as input, the mean of each gene is calculated, and the sum of the products of their differences is divided by the product of the square roots of their respective differences, resulting in a correlation coefficient ranging from -1 to 1. The closer the absolute value of the correlation coefficient is to 1, the more consistent the expression trends of the two genes under different treatment conditions. After pairwise calculations for all stress-induced upregulated genes, a gene correlation coefficient matrix is ​​obtained.

[0027] The stress-induced co-expression topology network was constructed based on the inter-gene correlation coefficient matrix. An edge connection condition was set where the absolute value of the correlation coefficient was not lower than an initial threshold of 0.90. Specifically, for any two stress-induced upregulated genes i and j, an undirected edge was established between gene i and gene j nodes if the absolute value of their Pearson correlation coefficient was not lower than 0.90. The edge weight was assigned the absolute value of the correlation coefficient. Gene pairs with correlation coefficients lower than 0.90 were not connected. Using all stress-induced upregulated genes as the node set and gene pairs satisfying the above edge connection condition as the edge set, the stress-induced co-expression topology network was constructed. The weight of each edge in the network reflects the similarity of the expression trends of the corresponding two genes under cadmium stress; a higher weight indicates a closer co-expression relationship between the two genes.

[0028] In one specific embodiment, step S3 involves calculating the weighted neighbor silencing rate for each node gene in the stress-induced co-expression topology network, including: The node degree of each node gene in the stress-induced co-expression topology network is calculated to obtain the set of direct neighbor nodes and the corresponding edge weights of each node gene. The set of neighboring silent genes for each node gene is obtained by performing intersection statistical processing on the set of direct neighbor nodes and the set of candidate genes for stress-induced-relief silencing. The weighted neighbor silencing rate of each node gene is obtained by dividing the sum of the edge weights corresponding to each gene in the set of silent neighbor genes by the sum of the edge weights corresponding to all genes in the set of direct neighbor nodes. When the node degree of a node gene is 0, the weighted neighbor silencing rate is assigned a value of 0.

[0029] Specifically, node degree calculation involves counting the number of directly connected neighboring nodes for each gene in the stress-induced co-expression topology network. For each gene g in the network, the rows corresponding to gene g in the inter-gene correlation coefficient matrix are traversed, and all genes with an absolute correlation coefficient of not less than 0.90 are extracted as the set of direct neighboring nodes of gene g. At the same time, the edge weights between gene g and each direct neighboring node are recorded, and the edge weights are the absolute values ​​of the corresponding Pearson correlation coefficients. The node degree is the number of genes in the set of direct neighboring nodes. The intersection statistics are performed based on the set of direct neighboring nodes of each gene and the set of candidate genes for stress-induced and alleviating agent silencing. The intersection statistics determine whether each gene in the set of direct neighboring nodes also appears in the set of candidate genes for stress-induced and alleviating agent silencing. Genes that appear simultaneously constitute the set of neighboring silent genes of that gene. The set of neighboring silent genes reflects the group of genes in the cadmium stress co-expression network that are directly connected to the gene in the node and are simultaneously selectively silenced by nitric oxide.

[0030] The weighted neighbor silencing rate is calculated by dividing the sum of the edge weights corresponding to each gene in the neighbor-silenced gene set by the sum of the edge weights corresponding to all genes in the direct neighbor node set by the sum of the edge weights. The numerator is the sum of the absolute values ​​of the Pearson correlation coefficients between the node gene g and each gene in its neighbor-silenced gene set, and the denominator is the sum of the absolute values ​​of the Pearson correlation coefficients between the node gene g and all its direct neighbors. The weighted neighbor silencing rate ranges from 0 to 1. A higher value indicates a higher proportion of co-expression intensity of genes silenced by nitric oxide among the direct neighbors of the node gene g; that is, genes with closer co-expression relationships with the node gene in its network neighborhood are more likely to be selectively silenced by nitric oxide. When the node degree of a gene is 0, meaning the gene has no direct neighbors in the network, the denominator is 0, and the weighted neighbor silencing rate is directly assigned a value of 0.

[0031] In one specific embodiment, in step S3, the weighted neighbor silencing rate is calculated for each node gene in the stress-induced co-expression topology network to obtain a topology feature matrix, including: Betweenness centrality calculation was performed on each node gene in the stress-induced co-expression topological network to obtain the ratio of the frequency of each node gene in the shortest path of all node pairs in the network. Clustering coefficients were calculated for each node gene in the stress-induced co-expression topology network to obtain the ratio of the actual number of edges between the direct neighbors of each node gene to the maximum possible number of edges. The node degree, betweenness centrality, clustering coefficient, and weighted neighbor silencing rate of each node gene are integrated into a four-dimensional feature vector. The topological feature matrix is ​​constructed with stress-induced upregulated genes as rows and the four-dimensional feature vector as columns.

[0032] Specifically, betweenness centrality measures the degree to which a node gene acts as an information transmission hub in a stress-induced co-expression topological network. This invention uses unweighted shortest path to calculate betweenness centrality. The unweighted shortest path refers to the path between two nodes with the fewest edges, regardless of edge weights. For each node gene g in the network, all nodes in the network other than node gene g are considered as unordered node pairs. Let the number of nodes other than node gene g in the network be N-1. Then the total number of all unordered node pairs is N-1 multiplied by N-2 divided by two. When there is no connected path between two nodes, this node pair is not included in the total number of unordered node pairs in the denominator. The number of paths passing through node gene g in the unweighted shortest path between all unordered node pairs included in the denominator is counted. This number is divided by the total number of unordered node pairs included in the denominator. The resulting ratio is the betweenness centrality of node gene g. The higher the betweenness centrality value, the stronger the path transmission function of the node gene in the network, and the more likely it is to be located at the intersection of multiple co-expression relationships.

[0033] The clustering coefficient measures the tightness of the connections between the direct neighbors of a gene node. For a gene node g, its set of direct neighbors consists of all genes whose absolute Pearson correlation coefficient with gene g is not less than 0.90. The actual number of connections between any two neighboring nodes in the set of direct neighbors is counted, and the absolute Pearson correlation coefficient between any two neighboring nodes is also not less than 0.90. Let the number of direct neighbors be k. The theoretical maximum number of edges that can exist in the set of direct neighbors is k multiplied by k minus one divided by two. The actual number of edges divided by the maximum possible number of edges is the clustering coefficient. When the number of direct neighbors k is less than 2, the maximum possible number of edges is zero, and the clustering coefficient is directly assigned to zero. The value of the clustering coefficient is between 0 and 1. The higher the value, the tighter the co-expression relationship between the direct neighbors of the gene node.

[0034] After calculating betweenness centrality and clustering coefficients, the four types of eigenvalues ​​of all stress-induced upregulated genes in the stress-induced co-expression topology network were subjected to min-maximum normalization. The minimum and maximum values ​​were statistically obtained within the range of all stress-induced upregulated genes. Specifically, each eigenvalue was subtracted from its minimum value among all stress-induced upregulated genes, and then divided by the difference between its maximum and minimum values ​​among all stress-induced upregulated genes. After normalization, all four types of eigenvalues ​​were mapped to between 0 and 1, eliminating the need for non-negative integer node degrees and compatibility issues with betweenness centrality and clustering. The coefficients, weighted neighbor silence rate, and other features are decimals between 0 and 1, representing dimensional differences. After normalization, a variance inflation factor test is performed on the four features to rule out multicollinearity. The variance inflation factor measures the degree of linear correlation between features; a higher variance inflation factor indicates a stronger linear correlation between that feature and other features. Features with a variance inflation factor exceeding 10 are considered to have severe multicollinearity. Principal component analysis (PCA) is used to reduce the dimensionality of features with severe multicollinearity. PCA linearly groups the multicollinear features. The principal component scores are combined to form mutually orthogonal features. After replacing the original collinear features with the principal component scores, the variance inflation factor of all features is re-examined until the variance inflation factor of all features does not exceed 10. After passing the test, the features are arranged in a fixed order of node degree, betweenness centrality, clustering coefficient, and weighted neighbor silencing rate to form a four-dimensional feature vector. The weighted neighbor silencing rate is calculated by weighting the absolute value of the Pearson correlation coefficient between node gene g and its direct neighbors. Unlike the three features of node degree, betweenness centrality, and clustering coefficient, which only describe the pure structural position of nodes, the weighted neighbor silencing rate directly introduces the member information of the stress-induced-relief silencing candidate gene set into the topological feature calculation, establishing a quantitative correlation between network topological position and selective silencing behavior. There is no definition and application of this feature in the existing technology. The four-dimensional feature vectors of all node genes are arranged in rows with stress-induced upregulated genes as rows and normalized fixed-order four-dimensional feature vectors as columns to construct a topological feature matrix. The number of rows in the topological feature matrix is ​​the total number of stress-induced upregulated genes, and the number of columns is fixed at four.

[0035] In one specific embodiment, step S4 includes: The topological feature matrix and the stress-induced-relief silencing candidate gene set were divided into training set and test set in a 7:3 ratio, and the ratio of genes within the candidate gene set to genes outside the candidate gene set was maintained at 1:1 in both the training set and the test set. The training set is input into the logistic regression model for parameter training. The regression coefficients are solved by maximum likelihood estimation. The test set is input into the trained logistic regression model for validation. The validation pass condition is that the area under the receiver operating characteristic curve is not less than 0.75 and the F1 score is not less than 0.70. The selective silence probability prediction model is obtained. The topological feature matrix of stress-induced upregulated genes is input into the selective silencing probability prediction model for extrapolation to obtain the silencing probability score matrix. The silencing probability score matrix is ​​then filtered with a silencing probability of not less than 0.70 to obtain a high-confidence silencing prediction gene set. The high-confidence silencing prediction gene set was mapped to KEGG database pathways. Enrichment scores were calculated for each pathway using the silencing probability of each gene in the high-confidence silencing prediction gene set as weights. After the pathway hypergeometric test value was corrected by the Benjamini-Hochberg method, a significant enrichment threshold of less than 0.05 was used. Significantly enriched pathways were sorted from high to low enrichment scores to obtain a priority list of pathways for remission agent selective silencing. Genes with a silencing probability of not less than 0.80, at least two occurrences of the pathway, and a weighted neighbor silencing rate of not less than 0.50 in the top 20% of the priority pathways for remission agent selective silencing were selected to obtain the core source gene set.

[0036] Specifically, the topological feature matrix and the stress-induced / alleviator silencing candidate gene set are divided into training and testing sets in a 7:3 ratio. The division method is to randomly sort all stress-induced upregulated genes and then split them in a 7:3 ratio. The training set contains 70% of the genes and their corresponding four-dimensional feature vectors and candidate gene set member labels, while the testing set contains the remaining 30% of the genes and their corresponding four-dimensional feature vectors and candidate gene set member labels. The candidate gene set member labels are binary labels. Genes belonging to the stress-induced / alleviator silencing candidate gene set are labeled as 1, and those not belonging are labeled as 0. During the division, the genes labeled as 1 and 0 are independently split in a 7:3 ratio and then merged to ensure that the ratio of the number of genes labeled as 1 to the number of genes labeled as 0 in both the training and testing sets is maintained at 1:1, avoiding sample imbalance from interfering with model training.

[0037] The four-dimensional feature vectors and corresponding binary labels of all genes in the training set are input into the logistic regression model for parameter training. The logistic regression model uses the sigmoid function to map the linear combination of the four types of features to probability values ​​between 0 and 1. The input to the sigmoid function is the weighted sum of the four types of feature values ​​and their corresponding regression coefficients, plus the intercept term. The regression coefficients are obtained through maximum likelihood estimation. Maximum likelihood estimation finds the combination of regression coefficients on the training set that maximizes the product of the probabilities of all sample observation labels. The closer the probability value of a sample labeled 1 is to 1, and the closer the probability value of a sample labeled 0 is to 0, the larger the likelihood function value. After obtaining the regression coefficients through iterative optimization, the four-dimensional feature vectors of all genes in the test set are input into the trained logistic regression model to obtain the values ​​for each gene in the test set. The predicted gene silencing probability is binarized into a predicted label using a classification threshold of 0.50. The area under the receiver operating characteristic (ROC) curve and the F1 score are calculated. The ROC curve area is obtained by plotting the true positive rate and false positive rate under different classification thresholds and then calculating the area under the curve. The F1 score is the harmonic mean of precision and recall. Precision is the proportion of true labels with a predicted label of 1 that are also true labels with a predicted label of 1. Recall is the proportion of true labels with a predicted label of 1 that are also true labels with a predicted label of 1. The model is considered validated when the ROC curve area is not less than 0.75 and the F1 score is not less than 0.70, and a selective silencing probability prediction model is obtained. If the model fails validation, the process returns to step S3, the threshold for constructing the co-expression network edge is adjusted from 0.90 to 0.85, and the process is repeated.

[0038] The four-dimensional feature vectors of all stress-induced upregulated genes are input into the selective silencing probability prediction model. The model outputs a silencing probability value between 0 and 1 for each gene. This probability value reflects the predicted likelihood that the gene will be selectively silenced by the alleviating agent. The silencing probability values ​​of all stress-induced upregulated genes are arranged by gene number to form a silencing probability score matrix. The number of rows in the silencing probability score matrix is ​​the total number of stress-induced upregulated genes, and the number of columns is one column. Genes with a silencing probability value of not less than 0.70 in the silencing probability score matrix are screened to obtain a high-confidence silencing prediction gene set. All genes in the high-confidence silencing prediction gene set were mapped to KEGG database pathways, which are a set of known gene metabolism and signal transduction pathways. An enrichment score was calculated for each pathway. The enrichment score is the arithmetic mean of the silencing probabilities of all genes from the high-confidence silencing prediction gene set in that pathway, multiplied by the pathway's gene coverage, and then multiplied by the negative logarithm of the pathway's hypergeometric test statistical significance. The pathway gene coverage is the number of genes from the high-confidence silencing prediction gene set in that pathway divided by the total number of genes in the high-confidence silencing prediction gene set. The hypergeometric test statistical significance is the original test value obtained by calculating the enrichment degree of genes in the high-confidence silencing prediction gene set using the hypergeometric distribution. The original test value is then subjected to Benjamini-Hoc... The Hochberg method corrects for the false discovery rate. The Benjamini-Hochberg method obtains the corrected false discovery rate by sorting all original test values ​​of all pathways in ascending order, multiplying them by the total number of pathways, and dividing by the correction coefficient of each ranking. A significant enrichment threshold of less than 0.05 is used. Pathways that meet the significant enrichment threshold are sorted in descending order of enrichment score to obtain a list of pathways with priority for selective silencing by remission agents. Among all genes in the top 20% of pathways in the list of pathways with priority for selective silencing by remission agents, genes that simultaneously meet the following three conditions are selected: silencing probability not less than 0.80, at least 2 pathways appearing in the top 20% of pathways, and weighted neighbor silencing rate not less than 0.50. This results in the core source gene set.

[0039] The above describes the method for tracing key genes for plant stress relief based on transcriptome data in the embodiments of this application. The following describes the system for tracing key genes for plant stress relief based on transcriptome data in the embodiments of this application. One embodiment of the system for tracing key genes for plant stress relief based on transcriptome data in the embodiments of this application includes: The sequencing module is used to perform high-throughput transcriptome sequencing on plant leaves from the control group, stress treatment group, alleviator treatment group, and combined treatment group. After quality filtering, clean transcriptome data for each treatment group are obtained. The screening module is used to assemble and deredundate the clean transcriptome data of each treatment group to obtain a whole genome expression matrix; and to perform intersection screening of stress-induced upregulated genes and reliever-silenced downregulated genes to obtain a stress-induced-relieving agent silencing candidate gene set. The construction module is used to construct a co-expression topology network for the stress-induced upregulated genes based on the stress-induced silencing candidate gene set as the silencing label, thereby obtaining a stress-induced co-expression topology network; and to calculate the weighted neighbor silencing rate for each node gene in the network to obtain a topology feature matrix, wherein the weighted neighbor silencing rate is the weighted proportion of genes belonging to the candidate gene set among the direct neighbors of the node gene to all direct neighbors; The calculation module is used to train the topological feature matrix and the candidate gene set using a logistic regression model to obtain a selective silencing probability prediction model; to extrapolate the stress-induced upregulated genes using the prediction model to obtain a silencing probability score matrix, and to screen to obtain a high-confidence silencing prediction gene set; and to calculate the high-confidence silencing prediction gene set by using the silencing probability as a weight through pathway enrichment to obtain a priority pathway ranking list for remission agent selective silencing, and to screen to obtain a core source gene set.

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

Claims

1. A method for tracing the origins of key genes for plant stress relief based on transcriptome data, characterized in that, The method includes: Step S1: High-throughput transcriptome sequencing was performed on the plant leaves of the control group, stress treatment group, alleviator treatment group and combined treatment group. After quality filtering, clean transcriptome data of each treatment group were obtained. Step S2: The clean transcriptome data of each treatment group are processed by transcript assembly and redundancy removal to obtain a whole genome expression matrix; the intersection of stress-induced upregulated genes and reliever-silenced downregulated genes is screened to obtain a stress-induced-relieving agent silencing candidate gene set; Step S3: Using the stress-induced / alleviating agent silencing candidate gene set as the silencing labeling basis, construct a co-expression topology network for the stress-induced upregulated genes to obtain a stress-induced co-expression topology network; calculate the weighted neighbor silencing rate for each node gene in the network to obtain a topology feature matrix, where the weighted neighbor silencing rate is the weighted proportion of genes belonging to the candidate gene set among the direct neighbors of the node gene to all direct neighbors; Step S4: Train the topological feature matrix and the candidate gene set using a logistic regression model to obtain a selective silencing probability prediction model; calculate the silencing probability score matrix by using the prediction model to obtain the stress-induced upregulated genes, and screen to obtain a high-confidence silencing prediction gene set; calculate the high-confidence silencing prediction gene set by using the silencing probability as a weight through pathway enrichment to obtain a priority pathway ranking list for remission agent selective silencing, and screen to obtain a core source gene set.

2. The method for tracing key genes for plant stress relief based on transcriptome data according to claim 1, characterized in that, Step S1 includes: At least three biological replicates were set up for each of the control group, stress treatment group, alleviator treatment group and combined treatment group. Plant leaves were collected from each treatment group, flash-frozen in liquid nitrogen and then stored at -80℃ to obtain frozen leaf samples from each treatment group. Total RNA was extracted from the frozen leaf samples of each treatment group using TRIzol reagent. The integrity of the extracted RNA was tested. RNA samples with an RNA integrity index of not less than 7 were processed for library construction to obtain sequencing libraries for each treatment group. The sequencing libraries of each treatment group were input into a high-throughput sequencing platform for paired-end sequencing to obtain the raw sequencing data of each treatment group. The CutAdapter tool was used to remove adapter sequences and filter low-quality bases from the raw sequencing data of each treatment group. The filtering quality threshold was set to a Phred value of not less than 20, a Q20 base ratio of not less than 97%, and a Q30 base ratio of not less than 94% to be considered qualified, thus obtaining clean transcriptome data for each treatment group.

3. The method for tracing key genes for plant stress relief based on transcriptome data according to claim 1, characterized in that, Step S2 includes: The clean transcriptome data of each treatment group were input into the Trinity assembly tool for de novo transcript assembly to obtain candidate transcript sequences. The candidate transcript sequences were then subjected to sequence similarity deduplication based on the CD-HIT tool. The sequence similarity threshold was set to 0.95, and the longest sequence was retained as the representative sequence to obtain a single gene sequence file. The single gene sequence files are input into the RSEM tool for expression level estimation to obtain the expected read values ​​of each gene in each treatment group. The whole genome expression matrix is ​​constructed based on the expected read values. The genome-wide expression matrix was analyzed using the DESeq tool, and four pairs of differential expression analyses were performed. The original test values ​​were corrected using the Benjamini-Hochberg method. With a false positive rate of less than 0.05 after correction as the threshold, the stress-induced upregulated gene sets were obtained by comparing the control group with the stress treatment group, the remission agent response gene sets were compared with the control group with the remission agent treatment group, and the remission agent silencing downregulated gene sets were compared with the stress treatment group with the combined treatment group. The stress-induced upregulated gene set and the alleviator-silenced downregulated gene set are subjected to intersection screening to obtain a stress-induced-alleviator-silenced candidate gene set.

4. The method for tracing key genes for plant stress relief based on transcriptome data according to claim 1, characterized in that, In step S3, using the stress-induced / alleviating agent silencing candidate gene set as the silencing labeling basis, a co-expression topology network is constructed for the stress-induced upregulated genes, including: Extract the expression quantum matrix of the stress-induced upregulated gene in all samples from the whole genome expression matrix; Based on the expression quantum matrix, the Pearson correlation coefficients of the stress-induced upregulated genes were calculated pairwise to obtain the intergene correlation coefficient matrix. Based on the correlation coefficient matrix between genes, edge construction is performed. The edge connection condition is that the absolute value of the correlation coefficient is not lower than the initial threshold of 0.

90. The stress-induced upregulated genes are used as nodes, and the gene pairs that meet the edge connection condition are used as edges to obtain the stress-induced co-expression topology network.

5. The method for tracing key genes for plant stress relief based on transcriptome data according to claim 4, characterized in that, In step S3, the weighted neighbor silencing rate is calculated for each node gene in the stress-induced co-expression topology network, including: The node degree of each node gene in the stress-induced co-expression topology network is calculated to obtain the set of direct neighbor nodes and the corresponding edge weights of each node gene. Based on the intersection statistics of the set of direct neighbor nodes and the set of candidate genes for stress-inducing-relieving agents, the set of neighbor silent genes for each node gene is obtained. The weighted neighbor silencing rate of each node gene is obtained by dividing the sum of the edge weights corresponding to each gene in the set of silent neighbor genes by the sum of the edge weights corresponding to all genes in the set of direct neighbor nodes; when the node degree of a node gene is 0, the weighted neighbor silencing rate is assigned a value of 0.

6. The method for tracing the origins of key plant stress relief genes based on transcriptome data according to claim 5, characterized in that, In step S3, the weighted neighbor silencing rate is calculated for each node gene in the stress-induced co-expression topology network to obtain a topology feature matrix, including: Betweenness centrality calculation is performed on each node gene in the stress-induced co-expression topology network to obtain the ratio of the frequency of each node gene in the shortest path of all node pairs in the network. Clustering coefficients are calculated for each node gene in the stress-induced co-expression topology network to obtain the ratio of the actual number of edges between the direct neighbors of each node gene to the maximum possible number of edges. The node degree, betweenness centrality, clustering coefficient, and weighted neighbor silencing rate of each node gene are integrated into a four-dimensional feature vector. The stress-induced upregulated genes are used as rows and the four-dimensional feature vector is used as columns to construct a topological feature matrix.

7. The method for tracing key genes for plant stress relief based on transcriptome data according to claim 1, characterized in that, Step S4 includes: The topological feature matrix and the stress-inducing-relieving agent silencing candidate gene set are divided into a training set and a test set in a 7:3 ratio, and the ratio of genes within the candidate gene set to genes outside the candidate gene set is maintained at 1:1 in both the training set and the test set. The training set is input into the logistic regression model for parameter training. The regression coefficients are solved by maximum likelihood estimation. The test set is input into the trained logistic regression model for validation. The validation pass condition is that the area under the receiver operating characteristic curve is not less than 0.75 and the F1 score is not less than 0.

70. The selective silence probability prediction model is obtained. The topological feature matrix of the stress-induced upregulated genes is input into the selective silencing probability prediction model for extrapolation to obtain a silencing probability score matrix. The silencing probability score matrix is ​​then filtered with a silencing probability of not less than 0.70 to obtain a high-confidence silencing prediction gene set. The high-confidence silencing prediction gene set was mapped to KEGG database pathways. Enrichment scores were calculated for each pathway using the silencing probability of each gene in the high-confidence silencing prediction gene set as weights. After correction using the Benjamini-Hochberg method, a significant enrichment threshold of less than 0.05 was used. Significantly enriched pathways were sorted from highest to lowest enrichment score to obtain a priority list of pathways for allergen-selective silencing. Genes with a silencing probability of at least 0.80, at least two occurrences, and a weighted neighbor silencing rate of at least 0.50 in the top 20% of pathways in the priority list were further screened to obtain the core source gene set.

8. A system for tracing the origins of key plant stress-relieving genes based on transcriptome data, characterized in that, For implementing the method for tracing key genes for plant stress relief based on transcriptome data as described in any one of claims 1-7, the system for tracing key genes for plant stress relief based on transcriptome data comprises: The sequencing module is used to perform high-throughput transcriptome sequencing on plant leaves from the control group, stress treatment group, alleviator treatment group, and combined treatment group. After quality filtering, clean transcriptome data for each treatment group are obtained. The screening module is used to assemble and deredundate the clean transcriptome data of each treatment group to obtain a whole genome expression matrix; and to perform intersection screening of stress-induced upregulated genes and reliever-silenced downregulated genes to obtain a stress-induced-relieving agent silencing candidate gene set. The construction module is used to construct a co-expression topology network for the stress-induced upregulated genes based on the stress-induced silencing candidate gene set as the silencing label, thereby obtaining a stress-induced co-expression topology network; and to calculate the weighted neighbor silencing rate for each node gene in the network to obtain a topology feature matrix, wherein the weighted neighbor silencing rate is the weighted proportion of genes belonging to the candidate gene set among the direct neighbors of the node gene to all direct neighbors; The calculation module is used to train the topological feature matrix and the candidate gene set using a logistic regression model to obtain a selective silencing probability prediction model; to extrapolate the stress-induced upregulated genes using the prediction model to obtain a silencing probability score matrix, and to screen to obtain a high-confidence silencing prediction gene set; and to calculate the high-confidence silencing prediction gene set by using the silencing probability as a weight through pathway enrichment to obtain a priority pathway ranking list for remission agent selective silencing, and to screen to obtain a core source gene set.

9. The system according to claim 8, characterized in that, High-throughput transcriptome sequencing was performed on plant leaves from the control group, stress treatment group, alleviator treatment group, and combined treatment group. After quality filtering, clean transcriptome data for each treatment group were obtained, including: At least three biological replicates were set up for each of the control group, stress treatment group, alleviator treatment group and combined treatment group. Plant leaves were collected from each treatment group, flash-frozen in liquid nitrogen and then stored at -80℃ to obtain frozen leaf samples from each treatment group. Total RNA was extracted from the frozen leaf samples of each treatment group using TRIzol reagent. The integrity of the extracted RNA was tested. RNA samples with an RNA integrity index of not less than 7 were processed for library construction to obtain sequencing libraries for each treatment group. The sequencing libraries of each treatment group were input into a high-throughput sequencing platform for paired-end sequencing to obtain the raw sequencing data of each treatment group. The CutAdapter tool was used to remove adapter sequences and filter low-quality bases from the raw sequencing data of each treatment group. The filtering quality threshold was set to a Phred value of not less than 20, a Q20 base ratio of not less than 97%, and a Q30 base ratio of not less than 94% to be considered qualified, thus obtaining clean transcriptome data for each treatment group.

10. The system according to claim 8, characterized in that, The clean transcriptome data of each treatment group were processed by transcript assembly and redundancy removal to obtain the whole genome expression matrix; Intersection screening of stress-induced upregulated genes and alleviator-silencing downregulated genes yielded a set of stress-induced-alleviator-silencing candidate genes, including: The clean transcriptome data of each treatment group were input into the Trinity assembly tool for de novo transcript assembly to obtain candidate transcript sequences. The candidate transcript sequences were then subjected to sequence similarity deduplication based on the CD-HIT tool. The sequence similarity threshold was set to 0.95, and the longest sequence was retained as the representative sequence to obtain a single gene sequence file. The single gene sequence files are input into the RSEM tool for expression level estimation to obtain the expected read values ​​of each gene in each treatment group. The whole genome expression matrix is ​​constructed based on the expected read values. The genome-wide expression matrix was analyzed using the DESeq tool, and four pairs of differential expression analyses were performed. The original test values ​​were corrected using the Benjamini-Hochberg method. With a false positive rate of less than 0.05 after correction as the threshold, the stress-induced upregulated gene sets were obtained by comparing the control group with the stress treatment group, the remission agent response gene sets were compared with the control group with the remission agent treatment group, and the remission agent silencing downregulated gene sets were compared with the stress treatment group with the combined treatment group. The stress-induced upregulated gene set and the alleviator-silenced downregulated gene set are subjected to intersection screening to obtain a stress-induced-alleviator-silenced candidate gene set.

Images