Method for constructing gene regulatory network based on multi-stage collaborative score and application thereof

Abstract

The invention relates to the technical field of bioinformatics, in particular to a method for constructing a gene regulatory network based on multi-stage collaborative scoring and application thereof, and the method comprises the following steps: firstly, constructing an RNA sequencing library at a plurality of time points and carrying out data processing to obtain the expression quantity of each gene at different time points; then, performing logarithmic transformation on the expression quantity of each pair of genes, and calculating a Pearson's correlation coefficient so as to quantify a linear relation between the genes; meanwhile, expression changes of adjacent time points are evaluated, the similarity of change amplitudes is calculated, and finally a comprehensive score is obtained. And carrying out average calculation on the scores to obtain the multi-stage expression cooperation score of the gene pair. And constructing a gene regulation network based on the scores. Complexity and functionality in organisms are effectively reflected, an innovative thought is provided for research of a gene regulatory network, and the method particularly has wide application potential in research of plants on an environmental stress mechanism.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and in particular to a method for constructing a gene regulatory network based on multi-stage synergy scores and an application thereof. Background Art

[0002] With the rapid development of genomics and transcriptomics, RNA sequencing (RNA-seq) technology has become an important tool for studying gene expression regulation. RNA-seq can not only provide high-throughput gene expression information, but also reveal complex transcriptional regulation phenomena between genes. In recent years, people have paid more and more attention to how to effectively construct gene regulatory networks from large amounts of RNA-seq data to deeply understand the interactions between genes and their functions in biological processes.

[0003] Traditional gene regulatory network construction methods are often based on linear correlations between genes. For example, the Pearson correlation coefficient is used to measure the similarity between gene expression patterns. The limitation of this method is that it assumes that the interactions between genes are linear, which may miss some nonlinear and complex interaction effects. For example, the regulation of a gene is not only affected by its own expression level, but may also be regulated by the expression of other genes. This complexity is difficult to capture through simple correlation analysis. In addition, when processing transcriptome data at multiple time points, traditional methods often fail to fully consider the dynamic changes of time series, resulting in inaccurate biological significance of network construction. Therefore, it is urgent to develop new methods to improve the construction accuracy and biological interpretation ability of gene regulatory networks.

[0004] The present invention proposes a new method based on multi-stage expression synergy score, which combines transcriptome data at multiple time points, comprehensively considers the linear relationship between genes and their expression changes at different time points, so as to more comprehensively reveal the complex interactions between genes. This method can not only improve the accuracy of network construction, but also provide new perspectives and tools for the study of plant adaptation mechanisms under environmental stress. In addition, there is still a large gap in the field of plant biology, especially in the study of gene regulation under environmental stress. The dynamic changes of gene regulatory networks of many plants in the face of environmental changes have not been fully understood. The method based on multi-stage expression synergy score proposed in the present invention can provide an effective solution for research in this field, promote in-depth exploration of plant gene regulatory networks, and thus provide a theoretical basis and technical support for molecular breeding to improve crop stress resistance and yield. Summary of the invention

[0005] The purpose of the present invention is to address the limitations of the gene regulatory network construction methods in the prior art, and to propose a method for constructing a gene regulatory network based on multi-stage synergy scores and its application. This method comprehensively considers the expression levels of genes and their changes at different time points, accurately captures the complex interactions between genes, and thus improves the accuracy and reliability of network construction.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] A method for constructing a gene regulatory network based on multi-stage synergy scores comprises the following steps:

[0008] Step 1: RNA-seq library construction and data processing:

[0009] Taking cotton leaf tissue under cold stress as an example, RNA-seq library construction was performed at multiple time points to obtain the expression information of each gene at different time points. The specific data processing process includes first, using trimgalore software to perform quality control on the original sequencing data to remove low-quality sequences and connector contamination; then, using Tophat2 software to align the quality-controlled data to the cotton reference genome to achieve accurate sequence alignment. Finally, StringTie software was used to calculate gene expression and obtain the transcript abundance of each gene, which was presented in the form of TPM (Transcripts Per Million) value.

[0010] Step 2: Logarithmic transformation of expression and correlation calculation:

[0011] 1. Perform logarithmic transformation on the expression of each pair of genes at each time point:

[0012] L ij =log2(E ij +∈)

[0013] Among them, L ij represents the logarithm-transformed expression of gene i at time point j, E ij is the original expression, ∈ is a small constant (such as 1×10 -6 ), used to avoid zero values ​​in logarithmic operations.

[0014] 2. Calculate the Pearson correlation coefficient:

[0015]

[0016] Among them, r ij represents the linear correlation between gene i and gene j, L ik and L jk are the logarithmic expression levels of gene i and gene j at all time points k, and is the corresponding average expression level.

[0017] Step 3: Expression change assessment:

[0018] 1. Calculate the difference in expression change of each pair of genes:

[0019] ΔE ij (t)=(E i (t)-E i (t-1))-(E j (t)-E j (t-1))

[0020] Among them, ΔE ij (t) represents the difference in expression between gene i and gene j at adjacent time points, E i (t) represents the expression level of gene i at time point t, E i (t-1) represents the expression level of gene i at time point t-1, E j (t) represents the expression level of gene j at time point t, E j (t-1) represents the expression level of gene j at time point t-1.

[0021] 2. Calculate the similarity of the change magnitude:

[0022]

[0023] Among them, S ij (t) represents the similarity of the expression change amplitude of gene i and gene j at time point t, ΔE i (t) represents the expression change of gene i at time point t, ΔE j (t) represents the expression change of gene j at time point t, and ∈ is a small constant (such as 1×10 -6 ) to avoid division by zero in calculations.

[0024] Step 4: Calculate the multi-stage expression synergy score:

[0025]

[0026] Among them, MECN Score represents the multi-stage expression synergy score, n represents the total number of time points, C k represents the correlation score at time point k, S k Indicates the similarity of the change amplitude at time point k.

[0027] A regulatory network was constructed based on the multi-stage expression synergy scores among genes for subsequent analysis.

[0028] Preferably, in step one, RNA-seq library construction and data processing include: RNA-seq library construction and sequencing; quality control of raw sequencing data using trim_galore software; alignment of data to a reference genome using Tophat2 software; and calculation of gene expression using StringTie software to obtain transcript abundance of each gene.

[0029] Preferably, in step 4, the filtering threshold of the multi-stage expression synergy score is set to 0.6 to ensure the stability and accuracy of the constructed network, thereby improving the credibility of the gene regulatory network analysis.

[0030] The present invention also provides an application of the above method in plants, specifically: the plant is cotton, and a gene regulatory network is constructed by analyzing the time-series RNA-seq data of cotton leaf tissue under cold stress conditions to accurately evaluate the interaction between genes.

[0031] By adopting the above technical scheme: taking the transcriptome data of cotton leaf tissue at six time points under cold stress as an example, the performance of the constructed multi-stage expression synergy score network is better than the correlation network, which is specifically reflected in the standardized clustering coefficient of the network of 4.4435 and the large and small world characteristics of 2.0585, which are higher than the 1.9440 and 1.2228 of the traditional network, reflecting the close interaction between genes and efficient information transfer. In addition, the modularity value is 0.3184, which is significantly higher than the 0.1872 of the traditional network, showing a stronger functional module recognition ability. The present invention not only improves the accuracy and reliability of gene regulatory network construction, but also provides strong data support for studying the response mechanism of plants under cold stress. In general, this method effectively reflects the complexity and functionality in organisms, provides innovative ideas for the study of gene regulatory networks, and has a wide range of application potential in the study of plant response to environmental stress mechanisms, helps to reveal the interaction between genes, and provides important support for in-depth research in related fields.

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

[0033] 1. The present invention can more comprehensively reflect the expression changes of genes at different time points through the calculation of multi-stage expression synergy scores, thereby accurately capturing the interactions between genes, improving the accuracy of network construction, and providing an important reference for subsequent functional studies.

[0034] 2. Taking the transcriptome data of cotton leaves at six time points under cold stress as an example, the multi-stage expression synergy score network constructed by the present invention is superior to the traditional correlation network in performance. Specifically, the standardized clustering coefficient of the network is 4.4435, which is significantly higher than the 1.9440 of the traditional network, which shows that the network of this method can better capture the close interactions between genes; the big and small world characteristics are 2.0585, compared with 1.2228 of the traditional network, showing that information transfer is more efficient; the modularity value is 0.3184, which is significantly higher than the 0.1872 of the traditional network, indicating that the network of this method is stronger in identifying functional modules. These results show that this method has more advantages in capturing complex gene interactions and provides more reliable data support for understanding the adaptation mechanism of plants in adversity. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 A schematic diagram of a regulatory network based on multi-stage expression synergy scores constructed for the present invention;

[0036] Figure 2 Schematic diagram of the regulatory network constructed for the traditional correlation method;

[0037] Figure 3 It is a schematic diagram for comparing the network performance indicators of the method of the present invention and the traditional correlation method. DETAILED DESCRIPTION

[0038] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with 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 protection scope of the present invention. The embodiments described in the present invention are only a part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments in the present invention, all other embodiments obtained by ordinary technicians in the field without making creative work are within the scope of protection of the present invention.

[0039] A method for constructing a gene regulatory network based on multi-stage synergy scores comprises the following steps:

[0040] (1) RNA-seq library construction and data processing

[0041] Three-week-old cotton seedlings were selected and placed in a low temperature environment of 4°C for cold stress treatment. To study the temporal expression changes of genes, leaf tissues were collected at six time points (0h, 2h, 4h, 6h, 12h and 24h). At least three biological replicates were prepared for each time point to ensure the reliability of the experimental results. Leaf samples were quickly frozen in liquid nitrogen and then stored at -80°C until RNA was extracted.

[0042] RNA was extracted using Omega Plant RNA Kit (Omega Bio-tek, Cat. no. R6827-01) according to the instructions of the kit. During the extraction process, the sample was first ground into a fine powder in liquid nitrogen, and then the lysis buffer was added and column purification was performed to remove proteins, DNA and other impurities. After extraction, the concentration and purity of RNA were detected by NanoDrop spectrometer, and the integrity of RNA was observed by agarose gel electrophoresis to ensure that the A260 / A280 value of the sample was in the range of 1.8 to 2.0, indicating that the RNA quality was good.

[0043] Next, an RNA sequencing library was constructed using the Illumina TruSeq RNA Sample Preparation Kit (NEB, Cat. no. E7530). First, mRNA enrichment was performed, and the mRNA in the total RNA was enriched by magnetic beads bound to poly-T chains. The enriched mRNA was fragmented to break it into small fragments of 200-300 bp in size, and then reverse transcribed to synthesize cDNA. Random primers were used to synthesize the second strand of cDNA, and the library was constructed by end repair and adapter connection. Subsequently, PCR amplification of the library was performed. After amplification, the size distribution and concentration of the library were detected using an Agilent 2100 bioanalyzer to ensure that it met the sequencing requirements.

[0044] After the library was constructed, high-throughput double-end 150bp sequencing was performed using the Illumina NovaSeq 6000 platform. The raw data (raw reads) generated by sequencing were first processed by the quality control tool trimgalore to remove low-quality reads and adapter sequences. The quality control steps included: filtering bases with a Phred quality value lower than 20, and removing the parts of the sequence containing adapter contamination. After quality control, the processed reads were aligned to the cotton reference genome (ZJU-v2.1 version) using Tophat2 software. In the alignment step, a certain number of mismatches were allowed to improve the sensitivity and specificity of the alignment. The aligned data were used to calculate the expression level of each gene using StringTie software, and transcript abundance (TPM) was used as a standardized indicator of expression. After the expression data was sorted, subsequent downstream analysis could be performed.

[0045] (2) Logarithmic transformation of expression and correlation calculation

[0046] 1. Logarithmic transformation of expression

[0047] In order to avoid errors caused by large differences in expression levels, the expression levels of each pair of genes were first logarithmically transformed:

[0048] Lij =log2(E ij +∈)

[0049] Among them, L ij represents the logarithm-transformed expression of gene i at time point j, E ij is the original expression, ∈ is a small constant (such as 1×10 -6 ) is used to avoid zero values ​​in logarithmic operations. The purpose of logarithmic transformation is to compress the expression in a large range so that the expression is more consistent with the normal distribution, which is convenient for subsequent correlation calculations.

[0050] 2. Pearson correlation coefficient calculation

[0051] For each pair of genes (gene i and gene j), linear correlation analysis was performed on the logarithmic transformed expression at different time points, and the Pearson correlation coefficient was used to measure the linear relationship between genes. The formula is:

[0052]

[0053] Among them, r ij represents the linear correlation between gene i and gene j, L ik and L jk are the logarithmic expression levels of gene i and gene j at all time points k, and is the corresponding average expression. The Pearson correlation coefficient here is used to quantify the similarity of the expression patterns of gene i and gene j at each time point. The value range is [-1, 1]. Positive values ​​indicate positive correlation, and negative values ​​indicate negative correlation.

[0054] (3) Expression change assessment

[0055] 1. Calculate the difference in expression change of each pair of genes

[0056] For each pair of genes i and j, the difference in the magnitude of their expression changes at adjacent time points is first calculated using the following formula:

[0057] ΔE ij (t)=(E i (t)-E i (t-1))-(E j (t)-E j (t-1))

[0058] Among them, ΔE ij (t) represents the difference in expression between gene i and gene j at adjacent time points, E i (t) represents the expression level of gene i at time point t, E i (t-1) represents the expression level of gene i at time point t-1, E j(t) represents the expression level of gene j at time point t, E j (t-1) represents the expression level of gene j at time point t-1. The calculation of the difference in the expression change amplitude can reflect the expression change trend of genes i and j at adjacent time points.

[0059] 2. Calculate the similarity of the change magnitude

[0060] The similarity of the expression change amplitudes of gene i and gene j at time point t is calculated by the following formula:

[0061]

[0062] Among them, S ij (t) represents the similarity of the expression change amplitude of gene i and gene j at time point t, ΔE i (t) represents the expression change of gene i at time point t, ΔE j (t) represents the expression change of gene j at time point t, and ε is a small constant (such as 1×10 -6 ) is used to avoid division by zero in calculations. The closer the similarity of the change amplitude is to 1, the more similar the expression change patterns of the two genes are.

[0063] (4) Calculation of multi-stage expression synergy scores

[0064] Based on the expression change patterns at multiple time points, the multi-stage expression synergy score (MECN Score) can be calculated, and the formula is as follows:

[0065]

[0066] Among them, MECN Score represents the multi-stage expression synergy score, n represents the total number of time points, C k represents the correlation score at time point k, S k It represents the similarity of the change amplitude at time point k. By averaging the correlation scores and change amplitude similarity at each time point, the co-expression score of the gene pair over the entire time interval can be obtained.

[0067] (5) Construction of gene regulatory network

[0068] When constructing a gene regulatory network, potential gene pairs are first screened based on the MECN Score. The threshold is set to 0.6, that is, when the MECN Score of a gene pair is greater than 0.6, the gene pair is considered to have a synergistic regulatory relationship. The screened gene pairs will be used as candidate edges of the gene regulatory network and the network will be constructed. The results are shown in Figure 2. Figure 1As shown, through this method, the constructed gene regulatory network can capture the synergistic expression relationship of gene pairs in multiple stages, which helps to better understand the gene regulatory mechanism in complex biological processes.

[0069] (6) Comparison of the network of this method and the correlation network

[0070] In order to verify the effectiveness of the proposed gene regulatory network construction method, this method was compared with the traditional correlation network-based method ( Figure 2-3 ), mainly from the following network characteristics for comparative analysis:

[0071] 1. Calculate the standardized clustering coefficient

[0072] The standardized clustering coefficient measures the degree of clustering of nodes in the network and reflects the connection between local nodes in the network. Its calculation formula is as follows:

[0073]

[0074] Among them, C norm is the standardized clustering coefficient, n is the total number of nodes in the network, e i is the actual number of edges between neighbors of node i, k i is the degree of node i, that is, the number of edges connected to the node, and m is the total number of edges in the network.

[0075] The larger the standardized clustering coefficient, the more nodes in the network tend to form small groups with high connections with their neighbors, that is, the network has stronger local aggregation.

[0076] 2. Computing big and small world characteristics

[0077] The small-world characteristic reflects the overall structural characteristics of the network, and its calculation formula is as follows:

[0078]

[0079] Among them, S small-world is the coefficient of the big and small world characteristics, C norm is the standardized clustering coefficient of the network, that is, the ratio of the actual clustering coefficient of the network to the clustering coefficient of the random network, L norm It is the standardized average path length of the network, that is, the ratio of the actual average path length of the network to the average path length of the random network.

[0080] The larger the small-world characteristic, the stronger the local connectivity (high clustering coefficient) and the shorter the global connectivity (low average path length) of the network, indicating the efficient information transmission capability of the network.

[0081] 3. Computational modularity

[0082] Modularity is used to measure the community structure of the network and reflects the distribution of genes in different functional modules in the network. The calculation formula is as follows:

[0083]

[0084] Among them, Q is the modularity value of the network, A ij Indicates whether there is an edge between nodes i and j. If so, A ij =1, otherwise A ij =0, k i and k j They represent the degree of nodes i and j, that is, the number of edges connected to the nodes, m is the total number of edges in the network, δ(c i ,c j ) is the Kronecker function. If nodes i and j belong to the same community, then δ(c i ,c j )=1, otherwise δ(c i ,c j )=0.

[0085] The larger the value of modularity, the more obvious the community structure in the network and the higher the degree of clustering of genes in different modules. Usually, the value of Q is between 0 and 1, and Q>0.3 indicates that the network has a significant community structure.

[0086] In summary, the present invention proposes a method for constructing a gene regulatory network based on multi-stage synergy scores, which can effectively capture the expression changes and mutual regulatory relationships of genes at multiple time points. Compared with the traditional correlation network-based method, this method shows significant advantages in multiple network characteristics, such as higher clustering coefficient, small-world characteristics and modularity, indicating that the constructed gene regulatory network has stronger biological interpretation capabilities.

[0087] The description and practice disclosed in the present invention are easy to think and understand for ordinary technicians in the technical field, and several improvements and modifications can be made without departing from the principles of the present invention. Therefore, modifications or improvements made without departing from the spirit of the present invention should also be regarded as the protection scope of the present invention.

Claims

1. A method for constructing a gene regulatory network based on multi-stage synergy scores, characterized in that: The following steps are involved: Step 1: construct RNA sequencing libraries at multiple time points and perform data processing to obtain the expression level of each gene at different time points; Step 2: Logarithmically transform the expression of each pair of genes and calculate their Pearson coefficient to quantify the linear relationship between genes: L ij =log2(E ij +∈) Among them, L ij represents the logarithm-transformed expression of gene i at time point j, E ij is the original expression, ∈ is a small constant used to avoid zero values ​​in logarithmic operations; Among them, r ij represents the linear correlation between gene i and gene j, L ik and L jk are the logarithmic expression levels of gene i and gene j at all time points k, and is the corresponding average expression; Step 3: Evaluate the expression changes at adjacent time points and calculate the similarity of the change amplitudes to obtain a comprehensive score: Among them, S ij (t) represents the similarity of the expression change amplitude of gene i and gene j at time point t, ΔE i (t) represents the expression change of gene i at time point t, ΔE j (t) represents the expression change of gene j at time point t, ∈ is a small constant used to avoid division by zero in calculation; Step 4: Average these scores to obtain the multi-stage expression synergy score of the gene pair: Among them, MECN Score represents the multi-stage expression synergy score, n represents the total number of time points, C k represents the correlation score at time point k, S k Indicates the similarity of the change amplitude at time point k; The gene regulatory network was constructed based on the above scores.

2. A method for constructing a gene regulatory network based on multi-stage synergy scores according to claim 1, characterized in that: In step 1, RNA-seq library construction and data processing include: RNA-seq library construction and sequencing; data quality control using trim_galore software; alignment of RNA-seq data to the cotton reference genome using Tophat2 software; and calculation of gene expression using StringTie software.

3. The method for constructing a gene regulatory network based on multi-stage synergy scores according to claim 1, characterized in that: In step 4, the filtering threshold of the multi-stage expression synergy score is set to 0.6 to retain significant gene interactions and ensure that the constructed gene regulatory network can effectively capture key interaction relationships.

4. Use of the method according to any one of claims 1 to 3 in plants.

5. The use according to claim 4, characterized in that: The plant is cotton, and the method is used to construct a gene regulatory network based on the time-series RNA-seq data of cotton leaf tissue under cold stress. The network can reveal the interaction between genes of cotton under cold stress conditions.

Images