Gene network construction method based on two-stage Granger causality and application

Abstract

The invention relates to the technical field of bioinformatics, in particular to a method for constructing a gene network based on two-stage Granger causality and application, and the method comprises the following steps: firstly, carrying out global Granger causality test based on gene expression time sequence data, and constructing an initial candidate network containing all potential regulation and control relationships; and then, refining the candidate network by adopting a dynamic pruning strategy based on event driving. And before the target gene expression has significant fluctuation, detecting whether the candidate regulatory gene has synchronous significant expression change or not, and verifying the causal driving force of the candidate regulatory gene. And only the regulation and control relationship passing the verification is reserved, so that wrong causal association is effectively removed, and optimization of the network structure is realized. The network constructed by the invention better conforms to a real biological regulation rule, and can better reflect the dynamics and directivity of a regulation relationship in a biological system; the method is especially suitable for the research fields of plant stress response, development regulation and control, system biology and the like, and has wide application prospects and practical values.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, and in particular to a method and application for constructing gene networks based on two-stage Granger causality. Background Technology

[0002] Gene Regulatory Network (GRN) analysis is a crucial tool for studying gene-gene interactions and their regulatory mechanisms in biological systems. It reveals functional changes in cells under different physiological states, signal transduction pathways, and key regulatory relationships in the formation of complex traits, thus providing a theoretical basis for understanding the mechanisms of organismal development, environmental adaptation, and response. In recent years, the rapid development of high-throughput sequencing technology has enabled us to acquire high-precision gene expression data at multiple time points, making time-series-based gene regulatory network inference possible and providing a key tool for elucidating the dynamic regulatory behavior of biological systems.

[0003] Traditional methods for constructing gene regulatory networks often rely on correlations between gene expressions, such as using Pearson correlation coefficients to infer gene associations. However, correlation-based networks often struggle to distinguish between direct and indirect regulation, tending to generate highly dense and directionless networks that fail to reflect the dynamic process of gene regulation over time. Furthermore, correlation does not imply causation, thus failing to reveal the direction of regulation or elucidate potential lag relationships between genes. To improve causal inference capabilities, some studies have employed Granger causality analysis, using time-lag models to determine whether changes in the expression of one gene can predict the future state of another. However, traditional Granger causality methods are susceptible to noise and parameter redundancy in small-sample, high-dimensional gene expression data, leading to a high number of false positive connections and poor network stability. Simultaneously, these methods often treat the entire time series as a single continuous process, failing to reflect the dynamic reconstruction characteristics of regulatory networks across different states when organisms are under stress or experience physiological events. Therefore, there is an urgent need to develop new methods to improve the reliability of causal inference and capture the dynamic changes of gene regulatory networks under different states.

[0004] This invention proposes a gene network construction method and its application based on two-stage Granger causality. The method first identifies all potential regulatory directions through global Granger causality analysis, and then combines this with an event-driven dynamic pruning strategy to filter indirect and unstable causal relationships, resulting in a sparser, well-defined, and biologically interpretable regulatory network. This method can improve the reliability of causal inference under time-limited conditions and effectively capture the dynamic reconstruction characteristics of the regulatory network during environmental stimuli or physiological events. The technology proposed in this invention can provide a new analytical tool for studying plant stress response mechanisms, help reveal the functional roles of key regulatory genes at different stages, and provide a theoretical basis and technical support for crop stress resistance improvement and molecular breeding. Summary of the Invention

[0005] The purpose of this invention is to address the shortcomings of existing technologies by proposing a two-stage Granger causality-based gene network construction method and its application. This method can construct gene regulatory networks with topological structures that are more consistent with biological characteristics from short-time-series gene expression data.

[0006] To achieve the above objectives, the present invention adopts the following technical solution:

[0007] A method for constructing gene networks based on two-stage Granger causality includes the following steps:

[0008] Step 1: Acquire and process gene expression data from multiple time points to obtain the expression level matrix of each gene at different time points;

[0009] Preferably, in step 1, acquiring and processing gene expression data includes: constructing an RNA-seq library and performing high-throughput sequencing; using data quality control software to perform quality control on the raw sequencing data; using alignment software to align the quality-controlled data to a reference genome; and using gene expression calculation software to calculate the expression level of each gene in each sample.

[0010] Here, we take cotton leaf tissue under cold stress as an example to obtain gene expression time-series data:

[0011] 1. Sample collection and RNA extraction:

[0012] Sufficient leaf samples were collected at multiple pre-set time points (0h, 2h, 4h, 6h, 12h, 24h) and immediately flash-frozen in liquid nitrogen to inhibit RNA degradation. Subsequently, total RNA was extracted using the Trizol method, and the integrity and purity of the RNA were detected by gel electrophoresis and spectrophotometry.

[0013] 2. RNA-seq library construction and high-throughput sequencing:

[0014] High-quality total RNA was used to enrich mRNA with magnetic beads and then fragmented. Using the fragmented mRNA as a template, the first strand of cDNA was synthesized via reverse transcription, followed by the synthesis of the second strand to form double-stranded cDNA. The double-stranded cDNA underwent end repair, A-tailing, and sequencing adapter ligation to construct an RNA-seq library suitable for high-throughput sequencing. After the library passed quality control, sequencing was performed using high-throughput sequencing platforms such as Illumina to generate raw sequencing data (Raw Reads).

[0015] 3. Raw data quality control and comparison:

[0016] Quality control: The raw sequencing data was processed using the trim_galore software to remove sequencing adapter sequences and low-quality bases (bases with a Phred value below 20) to ensure the accuracy of the analysis data.

[0017] Sequence alignment: Using Tophat2 alignment software, quality-controlled short reads (CleanReads) were aligned to the cotton reference genome to determine their origin position in the genome.

[0018] 4. Gene expression quantification:

[0019] Using StringTie software, transcripts were assembled based on sequence alignment results, and the expression abundance of each gene was calculated. Finally, the expression levels were normalized to TPM (Transcripts Per Million) values ​​and organized into a gene expression matrix M. raw The rows of this matrix represent genes, the columns represent different time points, and the values ​​in the matrix are the TPM values ​​of the corresponding genes at specific time points.

[0020] Step 2: Based on the expression matrix, construct a candidate network using a global Granger causality test, including:

[0021] 1. Data preprocessing:

[0022] First, the TPM values ​​were transformed using log2(TPM+1), and then the time series of each gene was normalized using Z-score to eliminate dimensional differences.

[0023]

[0024] Where x is the original expression level, μ is the mean of the gene expression time series, and σ is its standard deviation. The processed result is the gene expression matrix M used for analysis.

[0025] 2. Global Granger causality test:

[0026] For any pair of gene time series X(t) and Y(t) in matrix M, a vector autoregression (VAR) model is constructed to test whether X is a Granger cause of Y. The lag order of the model is set to p = 1. Compare the following two models:

[0027] Unrestricted Model:

[0028] Y(t)=c1+α1Y(t-1)+β1X(t-1)+ε1(t)

[0029] Where Y(t) and X(t) represent the expression levels of gene Y and gene X at time point t, respectively, c1 is a constant term, α1 is the coefficient of the lag term of gene Y, β1 is the coefficient of the lag term of gene X, and ε1(t) is the model error term.

[0030] Restricted Model:

[0031] Y(t)=c2+γ1Y(t-1)+ε2(t)

[0032] Where Y(t) and X(t) represent the expression levels of gene Y and gene X at time point t, respectively, c2 is a constant term, γ1 is the coefficient of the lag term of gene Y, and ε2(t) is the model error term.

[0033] The F-test is used to determine whether the coefficient β1 is significantly non-zero.

[0034]

[0035] Among them, RSS r and RSS u These are the sum of squared residuals for the constrained and unconstrained models, respectively, where T is the length of the time series and p is the lag order.

[0036] 3. Generate candidate networks:

[0037] Calculate the p-value based on the results of the F-test. Set a significance level threshold α. GC =0.05. Iterate through all gene pairs. If the p-value obtained from the test is less than 0.05, then gene X is considered a potential regulator of gene Y, and the edge X→Y is added to the candidate set. This ultimately yields a graph G. c =(V, E) c ), where V is the set of all genes, E c It is the set of all potential regulatory relationships that pass the Granger causality test.

[0038] Step 3: Perform event-based dynamic pruning on the candidate network to form a refined gene regulatory network, including:

[0039] 1. Target gene event detection:

[0040] On the time series Y(t) of the target gene Y, the Pelt algorithm based on the Akaike Information Criterion (AIC) penalty is used to detect change points, identify event points where the expression level changes significantly, and form the event point set S. Y =τ1, τ2, ... The algorithm finds the optimal point of change by minimizing a cost function:

[0041]

[0042] Where τ is the location of the change point to be detected, T is the total length of the time series, and y i Let i be the observation value of the time series at time point i. To evaluate data point y i The cost function is the degree of fit between the model and the segment it belongs to. βf(τ) is a penalty term used to penalize model complexity, based on the Akaike Information Criterion, where β is the penalty coefficient used to balance the model's goodness of fit and complexity.

[0043] 2. Validation of regulatory gene behavior:

[0044] For each event point τ of the target gene Y k ∈S Y To verify the predicted regulatory gene X at a time point (τ) before the event occurred. k -1) Whether significant changes in activity are observed:

[0045] |X(τ k -1)-X(τ k -2)|>θ change

[0046] Where X represents a potential regulatory gene, τ k This is a significant event point indicating a change in the expression of the target gene y. θ change It is a threshold for regulating gene changes dynamically determined based on data distribution, and its value is the 90th percentile of the magnitude of change in the expression level of all genes at all time points. If for set S... Y There exists at least one event point τ. k If the above conditions are met, then the candidate edge e = (X → Y) is considered to have been dynamically verified and is an effective control relationship.

[0047] 3. Generate a refined network:

[0048] All edges e that have passed dynamic validation are retained to form the final refined gene regulatory network G. f =(V, E) f ), where E fIt is a collection of all effective regulatory relationships, which can be used for subsequent analysis of biological mechanisms.

[0049] This invention also provides an application of the above-described method in the study of plant gene regulatory networks.

[0050] Preferably, the plant is cotton, and the gene regulatory network is constructed using time-series RNA-seq data of cotton leaf tissues under cold stress to reveal the interaction relationships between genes in cotton under cold stress conditions.

[0051] By employing the above technical solution, taking transcriptome data from cotton leaf tissue at six time points under cold stress as an example, the regulatory network constructed using the method of this invention significantly outperforms the traditional Pearson correlation network in key topological indicators. Specifically, in terms of network heterogeneity, the heterogeneity value of the method of this invention is 0.6375, higher than the 0.6250 of the traditional method, indicating that the differences in node degree distribution in the network are more obvious and more consistent with the characteristics of real biological networks; in terms of scale-free network goodness of fit (R0.05), the method of this invention has a heterogeneity value of 0.6375, which is higher than the traditional method's 0.6250, indicating that the differences in node degree distribution in the network are more obvious and more consistent with the characteristics of real biological networks; 2 In terms of network density, the method of this invention achieves a density of 0.0404, higher than the 0.0330 of the traditional method, and is closer to a scale-free topology. Regarding network density, the method of this invention has a network density of only 0.1693, significantly lower than the 0.3427 of the traditional method, reducing network complexity by approximately 50% and effectively eliminating a large number of indirect or redundant connections. In terms of network information transmission efficiency, the method of this invention has an average path length of 1.7394, shorter than the 1.8163 of the traditional method, indicating fewer propagation steps in the network and a more pronounced "small-world" characteristic, which better aligns with the regulatory characteristics of organisms rapidly responding to external stimuli. In summary, the method of this invention not only significantly improves the accuracy and reliability of gene regulatory network construction but also provides data support for the study of plant response mechanisms under cold stress, and offers a reliable technical tool for gene function analysis, identification of key regulatory factors, and molecular breeding research.

[0052] Compared with the prior art, the present invention has the following beneficial effects:

[0053] 1. This invention effectively reduces false positive connections caused by noise and indirect regulation by integrating global Granger causality testing and dynamic expression change screening. Taking cotton cold stress data as an example, the network density constructed by this invention is only 0.1693, which is about 50% lower than the 0.3427 of the traditional Pearson correlation network, successfully eliminating a large number of redundant connections; at the same time, the scale-free network fit (R0.05) is significantly improved. 2 The value increased from 0.0330 to 0.0404, indicating that the network topology is closer to that of a real biological network, thus significantly improving the recognition quality of core regulatory genes.

[0054] 2. By employing a dynamic pruning strategy, this invention can capture the most significant regulatory events in temporal expression, enabling the network to better reflect dynamic regulatory patterns. Taking cotton under cold stress as an example, the average path length of the network constructed by this method is 1.7394, shorter than the 1.8163 of the traditional method, and the signal transmission efficiency is improved by approximately 4.24%. This better reflects the characteristics of organisms responding rapidly under adversity, thereby enhancing the ability to analyze stress response mechanisms.

[0055] 3. The method of this invention is designed for time-series data and can capture dynamic patterns of gene regulatory relationships changing over time. It is particularly suitable for analyzing biological processes with stage-specific characteristics, such as stress responses and developmental processes. This method can extract significant regulatory events with limited time points, providing a powerful technical means for understanding the dynamic response mechanisms of biological systems. Attached Figure Description

[0056] Figure 1 A comparison graph of the degree distribution of the networks is constructed for the three methods (the method of this invention, Granger causal network, and Pearson correlation network);

[0057] Figure 2 A visualization of the gene regulatory network constructed using the two-stage method described in this invention;

[0058] Figure 3 A visualization of a gene regulatory network constructed using only global Granger causality tests (without dynamic pruning);

[0059] Figure 4 A visualization of a gene co-expression network constructed using the traditional Pearson correlation coefficient method;

[0060] Figure 5 A bar chart showing the quantitative comparison of the three networks across five key topological metrics. Detailed Implementation

[0061] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings, so that those skilled in the art can better understand the advantages and features of the present invention, thereby making a clearer definition of the scope of protection of the present invention. The embodiments described in this invention are only some embodiments of the present invention, not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0062] Example:

[0063] A method for constructing a gene regulatory network based on two-stage Granger causality includes the following steps:

[0064] (1) Data preprocessing

[0065] The transcriptome sequencing dataset used in this invention was obtained from the GSA database of the National Gene Bank (CNGB), accession number CRA019230. This experiment used cotton (Gossypium hirsutum) seedlings as material to study the dynamic changes in gene expression in their leaves under 4℃ low-temperature stress. Sampling was conducted at six time points: 0h, 2h, 4h, 6h, 12h, and 24h, with three biological replicates at each time point. The extracted RNA was sequenced using the Illumina high-throughput sequencing platform, yielding high-quality RNA-seq data. The data preprocessing workflow is as follows:

[0066] 1. Raw data quality control and comparison:

[0067] First, the raw sequencing data (Raw Reads) underwent rigorous quality control using the trim_galore software. This step mainly involved removing adapter sequences, filtering out low-quality bases with a quality score below 20, and discarding sequences shorter than 50 bp after processing to ensure the accuracy of subsequent analyses. The high-quality sequences (Clean Reads) that passed quality control were used for the next alignment analysis. The Tophat2 alignment tool was used to align the Clean Reads to the cotton reference genome. Tophat2 can effectively identify splicing events across exons and generate BAM (Binary Alignment Map) files containing complete alignment information, providing a foundation for subsequent transcript assembly and expression quantification.

[0068] 2. Quantification of gene expression:

[0069] The BAM files generated from Tophat2 alignments were processed using StringTie software. First, StringTie assembled transcripts based on the alignment results, constructing all expressed transcripts for each sample. Next, based on the assembly, StringTie estimated the expression abundance of each transcript. To eliminate the influence of sequencing depth and gene length on expression levels, TPM (Transcripts Per Million) was used as the standardized unit of expression. Finally, the TPM values ​​of each gene in all samples were summarized to construct a complete gene expression matrix. Each row of this matrix represents a gene, and each column represents a sample; the values ​​in the matrix represent the TPM expression level of the corresponding gene in the corresponding sample.

[0070] 3. Data transformation and normalization:

[0071] For subsequent analysis, the TPM expression matrix underwent a two-step processing. First, a log2(TPM+1) transformation was performed to mitigate heteroscedasticity and make the data distribution closer to a normal distribution. Second, Z-score normalization was applied to the time-series data of each gene to eliminate differences in expression dimensions between different genes.

[0072]

[0073] Where x is the expression level of the gene at a certain time point, and μ and σ are the mean and standard deviation of the gene expression level at all time points, respectively.

[0074] (2) Network Construction

[0075] 1. Construct a global causal candidate network:

[0076] All potential causal relationships were identified using the Granger causality test. For any pair of gene time series X(t) and Y(t), the following two vector autoregressive (VAR) models (lag order p = 1) were compared:

[0077] Unrestricted model:

[0078] Y(t)=c1+α1Y(t-1)+β1X(t-1)+ε1(t)

[0079] Where Y(t) and X(t) represent the expression levels of gene Y and gene X at time point t, respectively, c1 is a constant term, α1 is the coefficient of the lag term of gene Y, β1 is the coefficient of the lag term of gene X, and ε1(t) is the model error term.

[0080] Restricted model:

[0081] Y(t)=c2+γ1Y(t-1)+ε2(t)

[0082] Where Y(t) and X(t) represent the expression levels of gene Y and gene X at time point t, respectively, c2 is a constant term, γ1 is the coefficient of the lag term of gene Y, and ε2(t) is the model error term.

[0083] Use the F-test to determine whether the coefficient β1 of X(t-1) is significantly non-zero:

[0084]

[0085] Among them, RSS r and RSS u These are the sums of squared residuals for the constrained and unconstrained models, respectively.

[0086] T is the length of the time series, and p is the lag order. If the p-value of the F-test is less than 0.05, then the edge X→Y is added to the candidate network.

[0087] 2. Event-based dynamic pruning:

[0088] First, for each gene in the expression matrix M, calculate the magnitude of its expression change over all adjacent time points, i.e., |M i (t)-M i (t-1)|. The calculated magnitudes of change for all genes at all time points are aggregated into a set. The 90th quantile of this set is calculated and defined as the global change threshold θ. change This threshold represents the minimum standard for a significant change in expression levels.

[0089] Subsequently, target gene event points were identified. For candidate network G... c For each edge X→Y in the algorithm, the time series Y(t) of the target gene Y is extracted, and the Pelt algorithm is used to identify the set of event points S in which the expression level of Y undergoes significant mutation. Y =τ1,τ2,...,τ m The algorithm finds the optimal point of change by minimizing a cost function:

[0090]

[0091] Where τ is the location of the change point to be detected, T is the total length of the time series, and y i Let i be the observation value of the time series at time point i. To evaluate data point y i The cost function is the degree of fit between the model and the segment it belongs to. βf(τ) is a penalty term used to penalize model complexity, based on the Akaike Information Criterion, where β is the penalty coefficient used to balance the model's goodness of fit and complexity.

[0092] Next, the synchronization of regulatory genes was verified. This involved iterating through each event point τ of the target gene Y. k ∈S Y The test was conducted at the previous time point (τ). k -1) Does the change in the expression level of regulatory gene X meet the following conditions?

[0093] |X(τ k -1)-X(τ k -2)|>θ change

[0094] Where X represents a potential regulatory gene, τ k This is a significant event point indicating a change in the expression of the target gene Y. θ changeIt is a threshold for regulating gene changes that is dynamically determined based on data distribution. Its value is the 90th percentile of the magnitude of change in the expression level of all genes at all time points.

[0095] Then, pruning and retention are performed. For an edge X→Y, if at least one of the event points of the target gene Y passes the synchronicity verification, the edge is retained. Otherwise, it is pruned from the network.

[0096] Finally, all edges that pass the dynamic pruning verification described above are aggregated to form the final gene regulatory network G. f ( Figure 2 ).

[0097] 3. Granger causal network:

[0098] This method, as a control, constructs the network using only the global Granger causality test (P-value < 0.05) without dynamic pruning. This network represents a relatively broad set of causal relationships that have not been dynamically validated. Figure 3 ).

[0099] 4. Pearson Related Network:

[0100] This method represents a traditional co-expression-based network construction strategy. The Pearson correlation coefficient (PCC) of the expression profiles between all gene pairs (X,Y) is calculated:

[0101]

[0102] Where, r XY The Pearson correlation coefficient is the coefficient between gene X and gene Y, where T is the total number of time points, and X is the number of time points. i and Y i These represent the expression levels of gene X and gene Y at time point i, respectively. and These represent the average expression levels of genes X and Y at all time points, respectively. Figure 4 ).

[0103] (3) Network Comparison

[0104] To verify the effectiveness of the proposed gene regulatory network construction method, the following network topology indices were calculated and compared. Figure 1 , Figure 5 ):

[0105] 1. Network heterogeneity: This indicator reflects the degree of unevenness in the distribution of node degrees within a network. Networks with higher heterogeneity typically contain a few highly connected core nodes and a large number of low-degree connected peripheral nodes, which is more consistent with the characteristics of key genes playing a core regulatory role in real biological networks.

[0106]

[0107] Where H represents network heterogeneity, <k>This represents the average degree of all nodes in the network. <k 2 > represents the second moment of nodal degree.

[0108] 2. Scale-Free Network Fit (R) 2 This metric measures whether the degree distribution of a network follows a power-law distribution, a typical characteristic of scale-free networks. In biological networks, scale-free properties indicate robustness. The goodness-of-fit R is calculated by performing a linear regression on the logarithmic coordinates of the degree distribution. 2 The closer the value is to 1, the more the network conforms to the scale-free property:

[0109] log(P(k))=-γlog(k)+c

[0110] Where P(k) is the probability of a node of degree k appearing in the network, k is the node degree, γ is the power law exponent, and c is a constant.

[0111] 3. Network Density: This metric measures the ratio of the actual number of edges in a network to the theoretically maximum number of edges, reflecting the sparsity of the network. For a network containing N nodes and E edges:

[0112]

[0113] Where D is the network density, E is the number of actual edges in the network, and N is the total number of nodes in the network.

[0114] The lower the network density, the sparser the network, the fewer redundant connections, and the more precise the control relationships.

[0115] 4. Average path length: This refers to the average of the shortest path lengths between any two nodes in a network. It measures the average efficiency of information propagation in the network.

[0116]

[0117] Where L is the average path length, N is the total number of nodes in the network, and d(v i ,v j ) is node v i and v j The shortest distance between them. A shorter average path length is often associated with "small-world" network characteristics, indicating that the network can respond to signals quickly.

[0118] 4. Average Clustering Coefficient: This metric measures the degree of clustering of nodes in a network, that is, the tendency for a node's neighbors to also be neighbors. It reflects the modularity of the network. Local Clustering Coefficient of Nodes:

[0119]

[0120] Among them, C i Let k be the local clustering coefficient of node i. i E is the degree of node i. i It is its k i The number of actual edges between each neighboring node. The average clustering coefficient of the entire network is the sum of the values ​​of all nodes C. i The average value.

[0121] In summary, this invention proposes a method for constructing gene networks based on two-stage Granger causality and its application. Compared with traditional correlation-based methods, this method exhibits significant advantages in several network properties, such as higher network heterogeneity, better scale-free network fit, lower network density, and shorter average path length, indicating that the constructed gene regulatory network has stronger biological explanatory power.

[0122] The descriptions and practices disclosed in this invention are readily apparent and understandable to those skilled in the art, and various modifications and refinements can be made without departing from the principles of this invention. Therefore, any modifications or improvements made without departing from the spirit of this invention should also be considered within the scope of protection of this invention.< / k>

Claims

1. A method for constructing gene networks based on two-stage Granger causality, characterized in that, Includes the following steps: Step 1: Acquire and process gene expression data from multiple time points to obtain the expression level matrix of each gene at different time points; Step 2: Based on the expression matrix, construct a candidate network using the global Granger causality test, including: for any pair of gene time series X(t) and Y(t), using the residual sum of squares (RSS) of the restricted and unrestricted models. r and RSS u Perform an F-test to determine the significance of the causal relationship; Among them, RSS r and RSS u These are the residual sums of squares for the constrained model and the unconstrained model, respectively, where T is the length of the time series and p is the lag order; all causal relationships X→Y with P values ​​less than a preset threshold in the F-test are selected to form the candidate network. Step 3: Perform event-based dynamic pruning on the candidate network to form a refined gene regulation network, including: for each edge X→Y in the candidate network, firstly, on the time series Y(t) of the target gene Y, use a change point detection algorithm to identify event points τ where the expression level changes significantly; then, for each identified event point τ... k To verify that the regulatory gene X was at a time point τ before the event occurred. k-1 The magnitude of change in expression level ΔX(τ) k -1)=|X(τ k -1)-X(τ k -2) Is it greater than the preset significance threshold θ? change Finally, only edges that satisfy the verification conditions at least once are retained to form the final gene regulation network.

2. The method for constructing gene networks based on two-stage Granger causality according to claim 1, characterized in that, In step 1, acquiring and processing gene expression data includes: constructing an RNA-seq library and performing high-throughput sequencing; using data quality control software to perform quality control on the raw sequencing data; using alignment software to align the quality-controlled data to a reference genome; and using gene expression calculation software to calculate the expression level of each gene in each sample.

3. The method for constructing gene networks based on two-stage Granger causality according to claim 1, characterized in that, In step three, the change point detection algorithm finds the optimal change point τ by minimizing a cost function, wherein the cost function is: Where τ is the location of the change point to be detected, T is the total length of the time series, and y i Let i be the observation value of the time series at time point i. To evaluate data point y i The cost function is the degree of fit between the model and the segment it belongs to. βf(τ) is a penalty term used to penalize the complexity of the model, based on the Akaike Information Criterion, where β is the penalty coefficient used to balance the goodness of fit and complexity of the model.

4. An application of the method described in any one of claims 1-3 in the study of plant gene regulatory networks.

5. The application according to claim 4, characterized in that, The plant in question is cotton. The study used time-series RNA-seq data from cotton leaf tissues under cold stress to construct a gene regulatory network, in order to reveal the interactions between genes in cotton under cold stress conditions.

Images