EMD-based analysis of complexity with dissociated EEG amplitude and frequency information: a data-driven robust tool -for Autism diagnosis- compared to multi-scale entropy approach

1.   Introduction

Colorectal cancer (CRC), also named as colon cancer, is among top 3 common cancer diseases worldwide. With the fast development of powerful bioinformatics tools like microarray and high-throughput sequencing technique, the researchers are able to have a peek into the gene expression level of such disease. Omics data-based modeling and analysis does not only provide opportunity of systematically uncover the molecular mechanism underlying CRC development but also facilitate designing more effective treatments and prognosis.

Many studies focus on exploring the gene expression profile to select differentially expressed genes and connect their known gene function with the disease. In such analysis, genes are viewed as independent molecules and their connections are often ignored. As such, the conventional differential gene expression studies fall short on investigating key genes at systems level. Network-based methods have been used to analyze high-throughput transcriptomic data of cancer patients ^[1]. The prevalent methods of reconstructing gene networks from cross-sectional clinical transcriptomic data mainly include correlational methods, such as PCC-based methods ^[2], mutual information ^[3], regression methods ^[4] or machining-learning methods ^[5], which can only infer correlations or associations between genes. However, these network reconstruction methods cannot infer directions of the edges in the network which is important for detecting regulatory relationship between genes underlying disease development.

To reconstruct gene regulatory networks (GRNs) from gene expression data, many methods have also been developed in the past years ^[6]. The typical methods include Boolean network methods ^[7], ordinary differential equation (ODE) methods ^[8,9], Bayesian network methods ^[10]. These methods usually require time-course gene expression data ^[11] for inferring directed GRN. However, the temporal information is always not available for the clinical transcriptomic data of patients, which restricted to apply these methods.

On the other hand, many studies employed static network approach to bridge the connection between molecular profile and biological functions ^[1]. However, cancer is a dynamic evolving disease and the intrinsic molecular interactions should also have dynamic changing patterns. Therefore, constructing dynamic regulatory networks that are rewiring at different stages is a rational way to understand the mechanism underlying cancer development.

In this paper, to detect dynamic rewiring of gene regulation network along different disease stage of colorectal cancer, we reconstruct the stage-specific gene regulatory networks from a combined colorectal cancer dataset. We first used a supervised approach to select genes for network construction and analysis by integrating genes’ prognostic significance, functional association and correlation relationship and dynamic trends. An initial gene network was constructed using mutual information, which was then refined by employing a dynamic system model-based optimization approach to reconstruct the dynamic GRNs based on different colorectal cancer stages. The dynamic changes of the GRNs were analyzed to reveal important genes and their regulations during CRC development.

2.   Materials and methods

2.1.   Data collection and integration

In this article, four datasets of different sources are collected from both local hospital and gene expression omnibus (GEO) repository. GSE12945 is collected from Germany Max Planck Institute ^[12]; GSE14333 is collected from Australia Melbourne Royal Hospital ^[13] and GSE17536 is collected from United States Moffitt Cancer Center ^[14]. We also used a dataset of a Chinese cohort in our previous study, where 67 colorectal samples are processed through microarray analysis ^[15]. Both the gene expression data and clinical information of patients were collected.

After close examination of the datasets, we found that the sample number of different cancer grade is unequally distributed among the four datasets, in which GSE12945 does not include grade A and grade D samples while GSE14333 does not include grade D samples. The best way to solve the inequality of sample number is to combine them together as a largely merged dataset, where the sample number of each grade is more balanced. Using selected variables including GSE IDs, overall survival months, living status, gender, cancer grade, age at diagnosis and data source, we combined the 4 different datasets together, resulting in 529 samples with valid phenotype information. Detailed information can be found in Table 1. We used ‘removeBatchEffect’ in limma R package to normalize the merged data and to reduce the between-experiments batch effects, which is used in the downstream analysis.

2.2.   Mutual information estimation

Mutual information can be used to measure the nonlinear correlation relationship between two genes, which is a bilateral relation measurement. We use $\text{I}\left(\text{A}, \text{B} \right)$ to denote the mutual information between gene A and gene B. If $\text{I}\left(\text{A}, \text{B} \right) = 0$, then two gene expression patterns are uncorrelated, while larger $\text{I}\left(\text{A}, \text{B} \right)$ stands for closer pattern similarity and potentially implies close biological relationship.

In order to construct a mutual information network, we first calculate the mutual information value for each pair of genes, then a threshold is set to construct the connection network. In this step, the estimation method of mutual information plays a pivotal role in the performance of the ultimate mutual information network. Currently, two methods are commonly used to estimate the mutual information. First is to increase the accuracy of the estimation of joint probability density function from probability density estimation, such as histogram method and kernel density method ^[16]. The second way is to avoid the estimation of joint probability density function, such as K-nearest neighbor (KNN) method ^[17].

In terms of gene expression data, kernel density estimation is a better way because of its advantage over the high-dimensional dataset and the intrinsic normal distribution pattern of gene expression values. Define ${{X}_{1}}, {{X}_{2}}, \cdots, {{X}_{n}}$ as independent identically distributed random variables with joint probability density function $f\left(x \right)$, then the kernel density estimation for function $f\left(x \right)$ is

where $K\left(· \right)$ is called kernel function, and h is a pre-defined positive number, often referred as bandwidth funcion.

From the calculation of mutual information above, we know that the estimation for mutual information $I\left(X, Y \right)$ is $\hat{I}\left(X, Y \right)$, and it reads

2.3.   Network sparsification

Based on the above MI estimation of each pair of genes, we apply data process inequality (DPI) algorithm ^[18] to eliminate indirect associations between genes. In this way, a sparse network could be constructed, which is more suitable for the following regulatory network estimation.

Given the following gene interactions: $X\to Y\to Z$, where $Z$ is correlated with Y but not with X, then the DPI algorithm considers that

That is, if (X, Y) and (Y, Z) are correlated directly but (X, Z) is indirectedly correlated through Y, then $I\left(X, Z \right)\le I\left(X, Y \right)$, and $I\left(X, Z \right)\le I\left(Y, Z \right)$. Using the DPI algorithm, we could exclude those indirect correlations, and retain those estimated to be directly associated gene pairs. We use ${{e}_{ij}}$ to denote the edges in the MI network after network sparsification using DPI algorithm. ${{e}_{ij}} = 1$ if there is a direct association between gene i and gene j, otherwise ${{e}_{ij}} = 0.$

2.4.   Dynamic GRNs reconstruction

Assume ${{x}_{i}}\left(t \right)$ to be the expression level of gene i at time t, then the ordinary differetial equation (ODE) of transcriptionl kinetics would be

where ${{\alpha }_{i}}$ is the mRNA turnover rate, ${{R}_{i}}$ is the set of regulators of gene i and ${{\beta }_{ij}}$ is the regularotry strength from gene j to gene i. We incorporate ${{e}_{ij}}$ derived from the MI network into the above ODE model to exclude indirectly-correlated gene pairs from the network model for the purpose of reducing computational cost in the following parameter estimation steps. If the time-course data are provided, this equation can be utilized to construct a dynamic gene regulation networks (GRN), and the coefficients can be estimated by a linear regression method. The coefficients can also be used to demonstrate the gene interaction in a GRN.

Many available cancer datasets for research use is stage-static data, in which all the case with the same cancer progression grade will be counted without time-course data. Thus, it makes people hard to examine the dynamic gene evolving process. In order to overcome such obstacle, we employed a dynamic cascaded methods (DCM) algorithm ^[19] (Text S1) to reconstruct the dynamic gene networks from sample-based transcriptional data.

According to the intra-stage steady-rate assumption and the continuity assumption ^[19], we are able to build up a model that connects two consecutive stages in terms of gene profile.

where $\bar{x}_{i}^{\left(s \right)}$ and $\bar{x}_{i}^{\left(s-1 \right)}$ stand for the mean expressions of gene i in stage s and s-1. In addition, the coefficients in this formula is closely related to the previous ODE formula:

where ${{L}^{\left(s \right)}}$ is the time length of stage s, ${{\alpha }_{i}}$ is the degradation rate, and $\beta _{ij}^{\left(s \right)}$ is the regulatory coefficient.

In order to describe the dynamic pattern within a stage, we should define the $\lambda$ fraction of a stage as the proportional interpolation between the earliest and latest time points of the stage. According to the previous two assumptions, gene expression should be of linear trend change continuously, hence the largest and smallest gene expression should be on the either end of a gene expression stage. $t_{1}^{\left(s \right)}$ and $t_{{{N}^{\left(s \right)}}}^{\left(s \right)}$ denotes the earliest and the latest time points of stage s. The time of the $\lambda$ fraction in stage s can be expressed as $t_{\lambda }^{\left(s \right)} = \text{ }\!\!~\!\!\text{ }t_{1}^{\left(s \right)}+\lambda ·\left(t_{{{N}^{\left(s \right)}}}^{\left(s \right)}-t_{1}^{\left(s \right)} \right)$. After certain transformations, the dynamic model equation can be written as:

In this equation, gene expressions at the same ($\lambda$) fraction of two consecutive stages are connected. It describes the inter-stage dynamics of the GRN. And the intra-stage dynamical GRN can be described by the former ODE equation.

The details of the DCM algorithm are described in Text S1. All the steps listed above are performed in R (version 3.3.3) using packages including parallel, annotate, hgu133a.db, igraph, and most importantly, glmnet package that allows to fast implementation of elasticnet regression in selecting appropriate penalty coefficient and equation coefficients.

3.   Results

3.1.   Description of the merged datasets

In the merged dataset, there are about an equal number of male and female samples, counting for 54.4 and 45.6% respectively. The median age in those patients is 66 years old, corresponding to the higher hazard rate for old people to suffer from colorectal cancer. There are also young patients, who are only 19 years old, while the oldest patients are in their 90s. When it comes to the number of samples for each cancer grade, Grade B has the largest group of samples, while Grade D has the lowest group of samples, which is only 15. This phenomenon is also reasonable considering the terrible prognosis and lack of proper treatment for colorectal cancer in its late days. Throughout the whole survey, there are in total 262 recorded death cases and the median survival time length by day unit is 1161.1.

3.2.   Gene set selection and enrichment

To reconstruct dynamic GRNs during CRC progression, we focus on cancer progression-related genes that should be functionally important and associated to the survival of patients. So, we used Cox PH model in this step to select the survival-associated genes for downstream analysis. The log-rank test p-value was calculated to assess the significance of the association between a given gene and survival of patients. We selected 1868 genes with p-value less than 0.01 for the downstream analysis, to reduce probability of Type I error.

We next used Gene Ontology to perform gene set enrichment analysis ^[20]. The most common method for GO enrichment analysis is the Fisher Exact Test, where contingency tables are made for each GO terms about the number of genes have or don’t have the annotation and whether or not the genes are considered significant. 1868 survival-associated genes (log-rank p value < 0.01) were subject to gene set enrichment analysis. ‘topGO’ R package was used to perform the GO enrichment analysis and the top 10 GO functions are selected in Table 2, ranked by Kolmogorov-Smirnov test with more conservative ‘elim’ algorithm, which will render a more confident enrichment analysis result. From Table 2 we can see that the most enriched GO terms are closely related to the regulation of ubiquitin-protein, cell division, DNA elongation and replication and protein transportation, which have been verified to be closely related to cancer progression. For example, the deregulated ubiquitin ligases (E3s) have been experimentally validated to play important roles in tumor progression and cancer drug resistance ^[21]. For another example, protein trafficking has been shown to facilitate invasion and metastasis of CRC cells, and impaired transportation of member receptors represents important mechanisms of resistance to CRC-targeted therapy ^[22]. Furthermore, cell division and DNA elongation and replication are well-known biological processes associated with progression of many cancers including CRC ^[23].

To select genes with known direct interaction or indirect associations and narrow the range to select both statistically and biologically significant genes in colorectal cancer, we used STRING database (https://string-db.org/) to construct a PPI network using 1868 genes identified in the previous step as input. The PPI link tsv file is downloaded and 249 unique genes are selected in consideration of potential biological meanings.

Because of the multiple Affymetrix IDs for the same gene symbol, we here have more than 249 genes selected out of the merged gene expression matrix. Eventually, a matrix used to plot heatmap is selected with 539 genes as rows and 529 patients as columns. The heatmap (Figure 1) shows that these genes are largely correlated with living status, which is shown in the first row of column side colors (green for live and red for death). Many living patients are associated with the higher expression level of these 539 selected genes.

3.3.   Mutual information network

Using 249 genes selected from the previous STRING PPI network, we calculated mutual information for all pairs of genes and used the above-mentioned DPI algorithm to exclude indirectly associated gene pairs. In order to demonstrate the network structure, the Fruchterman-Reingold algorithm is used to automatically display the network as implemented in igraph R package. Because of the relatively large number of genes in our network, it’s inappropriate to show them all together without a certain difference, hence we modify the text size, edge width and node size in the plot to help better visualize the weight of each gene. The text size is set to be proportionate to the ratio of node degree over maximum node degree of the network; edge width is set to be proportionate to the ratio of log transformed edge weight (which is also the mutual information value between gene pairs) over the log transformed maximum edge weight; and the node size is set to be the ratio of the sum of the mutual information value between one gene and all the other genes, over the maximum number of it.

The visualization of the constructed mutual information network is shown in Figure 2. We also examined the degree distribution (Figure 3a) and cumulative degree distribution (Figure 3b) for the mutual information network. Most genes (about 45%) have only one connection, verifying the sparsity of the constructed MI network using ARACNE method and DPI algorithm. High-degree genes in gene networks are often considered as hub genes or potential key genes which may possess important biological functions. In our analysis, we found that gene CCT7 does not only has the largest text size but also it has the node size (Figure 2), indicating that CCT7 has both the largest connection with other genes and the sum of mutual information value. Importantly, CCT7 has been experimentally validated to be closely associated with some cancers, such as endometrial carcinoma ^[24]. This generates hypothesis of CRC as a key gene in CRC for experimental tests.

3.4.   Reconstruction of dynamic GRNs

The mutual information network enables us to have a big picture of the gene association within CRC, however, it cannot provide more information about the association direction or the level of such association. To investigate the gene regulatory roadmap underlying CRC, we next infer the directed regulatory networks for each stage of CRC.

A smaller gene expression matrix (with 127 Affymetrix probes and 529 patients) is extracted from the original matrix (with 22277 Affymetrix probes and 529 patients). The genes are selected based on the following criteria: Survival association significance assessed (genes with log-rank test p-value less than 0.01), PPI interaction information (249 genes in String database) and mutual information network information (genes with the sum of mutual information > 1). The column number of the gene expression matrix corresponds with the phenotype dataset (which has 529 patient entries). Sample-based gene-profiling data of different colorectal cancer grade is created with 61 Grade A, 276 Grade B, 177 Grade C and 15 Grade D patients. The difference in the cancer grade represents the stepwise cancer progression process from early stage to late stage.

To determine the gene expression evolving direction, a gene evolving trend analysis was conducted. The overall connecting error is used to help to solve this problem. It is defined as the L1-norm of all the individual connecting errors:

where $t_{1}^{\left(s \right)}$ and $t_{{{N}^{\left(s-1 \right)}}}^{\left(s-1 \right)}$ are the starting times of stage s and the ending time of stage s-1 respectively, $x\left(t_{1}^{\left(s \right)} \right)$ and $x\left(t_{{{N}^{\left(s-1 \right)}}}^{\left(s-1 \right)} \right)$ are the corresponding gene expression level of one gene. According to the continuity assumption, $x\left(t_{1}^{\left(s \right)} \right)$ and $x\left(t_{{{N}^{\left(s-1 \right)}}}^{\left(s-1 \right)} \right)$ should be the same. However, in practice, there could be error at the connecting point of two stages. By summing up the overall connecting error, we can minimize it in order to get the shortest path from grade A to grade D. Then, this path is used to determine the gene evolving trend within each stage. Figure S2 illustrates the possible combinations between each stage.

The gene trend analysis is realized in R by using igraph package. First, the grade separated gene-profiling data is generated as large list object, in which every list in the root is the expression profile, and four sublists contain different gene expressions in different cancer stages. Then, quantile function is used to select genes according to a different $\lambda \text{ }\!\!~\!\!\; \text{ fraction}$ in a stage. $\lambda$ is selected to be 0, 0.5, and 1 for interpolation. As for extrapolation, ${{\lambda }_{head}}$ and ${{\lambda }_{tail}}$ is set to be both at 4 for better controlling the linearity of the model equations. After that, a bootstrapping method is used to exploit the limited samples. We heuristically set the bootstrap sample size to half of the selected genes, and a total of 6150 bootstrap groups are generated for all the genes, and a total of 30750 corresponding model equations are obtained. The model coefficients in Eq (7) are then estimated by a regularized regression method. In this study we used elastic-net method to guarantee both sparsity and robustness of the estimation. We formulated the following objective function with a penalty to estimate ${{b}_{ij}}^{\left(s \right)}$ for each s,

where ${\mathsf{α}}_{ela}$ is an elasticnet mixing parameter (between 0 and 1), and θ is the penalty coefficient. ${{P}_{1}} = \underset{i}{\overset{{}}{\mathop \sum }}\, {{\left(a_{i}^{\left(s-1, s \right)} \right)}^{2}}+~\underset{ij}{\overset{{}}{\mathop \sum }}\, {{\left(~{{b}_{ij}}^{\left(s-1, ~s \right)} \right)}^{2}}+~~\underset{ij}{\overset{{}}{\mathop \sum }}\, {{\left(~{{b}_{ij}}^{\left(s \right)} \right)}^{2}}~$, and ${{P}_{2}} = \underset{i}{\overset{{}}{\mathop \sum }}\, \left| a_{i}^{\left(s-1, s \right)} \right|+~\underset{ij}{\overset{{}}{\mathop \sum }}\, \left| ~{{b}_{ij}}^{\left(s-1, ~s \right)} \right|+~~\underset{ij}{\overset{{}}{\mathop \sum }}\, \left| ~{{b}_{ij}}^{\left(s \right)} \right|~$.

Cross validation strategy is also used to select the optimal θ. After getting the ${{b}_{ij}}^{\left(s \right)}$ matrix, we then convert them to $\beta _{ij}^{\left(s \right)}$ matrix according to previous formula Eq (6). Then a gene regulatory network (GRN) is constructed according to these $\beta _{ij}^{\left(s \right)}$. The network construction process is repeated 50 times, and ultimately 50 networks are obtained. The confidence of a connection is calculated as occurring frequency among 50 networks, and network connectivity ${{P}_{ER}}$ is defined as the number of connections reserved in a network.

A sequence of ${\rm{\alpha}}_{ela}$ is selected in our model to optimize the sparse gene regulatory network. Confidence-Connectivity plot (C-C plot) is used to demonstrate the relationship between the threshold confidence level of a network and the network connectivity. From the Figure S3, we found that the GRNs of Grade A and C have a similar curve, while Grade B is a little higher in the network confidence level as the network connectivity increases. Grade D curve has a similar descending curve to Grade A and C in the beginning, but then diverts from all the other 3 stages and descends more slowly. Among 15006 possible connections (with direction) for a set of 123 genes in colorectal cancer, most connections have relatively low confidence, in other words, low frequency of occurrence. The connection confidence decreases following an exponential trend as the network connectivity increases. By setting a cutoff value of 30% confidence level, there is no more than 40% connectivity in either grade of the GRN (2615, 4999, 3157 and 4542 connections in Grade A, B, C, and D respectively).

Four stage-specific gene regulatory network (GRN) are reconstructed with the connectivity setting of 0.02. Four GRNs are then used by Cytoscape (http://www.cytoscape.org) and shown in Figure 4. In these figures, the node size is set to be proportionate to the sum of the products of the confidence and strength over all the direct incoming and outgoing connections. This product stands for the activity of a gene in the network. Whether the strength is for upregulation and downregulation of the target gene is marked with either an arrow (→) or a stop (┤). The level of the orange color of an edge is proportionate to the confidence level of the connection, and the width of the line is set corresponding to the strength of this regulation. The circular layout is ordered by the degree of the nodes.

3.5.   Analysis of dynamic GRNs

In the following, we discuss the biological implications of the reconstructed dynamic GRNs for CRC at each stage.

During CRC stage A (Figure 4a), gene RPS15 has a big positive regulation over RPL37, and the former is negatively regulated by SFSWAP, MAPK4 and other two genes. One known connection with high confidence level is the negative regulation from MYH2 to ACTN2.

During stage B (Figure 4b), RAPGEF1, ACTN2, MYH2, TGFA, MAPK4, SMG7, CHERP and DYRK1B are potential genes with biological significance. On literature review, RAPGEF1 gene have its first intron under somatic demethylation of a relaxed-criterion CpG island in 40% of colon cancers ^[25]. ACTN2, as isoform of ACTN4, whose amplifications has been verified to be associated with worse outcomes of cancer patients ^[26], is also believe to be critical to the outcome of CRC patients. MYH2 was found to be highly up-regulated (with 63-fold change) in CRC patients with tumors of high intrinsic COX-2 expression. TGFA is of the same family of TGF-beta, whose dysregulation has been shown to be closely related to colorectal cancer ^[27]. CHERP was found in Colorectal Cancer Atlas ^[28] with functions of negative regulation on cell proliferation.

In CRC stage C (Figure 4c), there is an interesting pattern of negative regulation from gene BAK1 to several other genes such as RAPGEF1, MYH2. Although BAK1 has relatively small value of product sum (corresponding to gene activity), it has a broad regulation function on multiple genes. Another interesting pattern to notice is that genes with low activity might have a strong regulation effect on other genes, such as RPL37. This is reasonable because important regulatory genes are not always active while their effects could be tremendous.

Figure 4d shows the stage D-specific regulatory network. The previous observation reoccurred that two low-expressed genes RPL37 and NPM1 have strong influence on the expression of each other as well as being regulated by many other active genes. From this pattern we hypothesize that those low-expressed but influential genes might be the downstream gene regulated by multiple less influential but high-expressed genes. The additive effect of the upstream genes might largely affect the downstream biological process.

4.   Discussion

In this study, we employed an integrated approach to analyze clinical transcriptomic data of CRC patients. We are interested at the directional gene regulatory networks for each stage and their dynamic changes during CRC development. Driven by this question, we collected and merged 4 datasets of CRC containing both gene expression data and clinical information. To infer GRNs from the large-scale clinical gene expression dataset, the survival-associated genes were selected as input for network construction, and an initial mutual information network was constructed to reduce the computational cost of the parameter estimation of the dynamical system model of GRN. By introducing a new regularized optimization method into the DCM algorithm, directional GRNs were inferred for each stage of CRC.

The DCM algorithm ^[19] employed in this paper relies on two assumptions: The intra-stage steady-rate assumption and the continuity assumption. The intra-stage steady-rate assumption or linear-dynamic assumption, models the dynamic gene expression profile with a linear trend within each stage of a process. The continuity assumption assumes that the gene expression trajectory is continuous both within each disease stage and at the connecting point of two adjacent stages. Obviously, these two assumptions are simplification of the real gene expression in patients. In reality, it’s largely plausible that the gene expression dynamics should be nonlinear, and the gene mutations or other genetic variations may result in discrete transition of gene expression trajectory. Therefore, the linear-dynamic assumption and the continuity assumption need to be generalized or waived to be more realistic.

In addition, the number of samples in each stage (stage A to D) is uneven in our study. The minimal sample number in that study is 13, which is comparative to 15 in our paper. Although the previous study has shown that the DCM algorithm is not sensitive to the sample size (Figure 3B in ^[19]), the uneven sample number could be considered for improving the DCM algorithm. For example, this consideration may be addressed by assigning different penalties or weights for different groups according to the number of samples.

In the current study, we focus on GRNs rewiring during CRC development from stage A to stage D. We did not consider normal samples of CRC since they are not available in the collected dataset, which is another limitation of our study. From the viewpoint of dynamic theory, disease progression can be considered a state transition, from a normal state, to pre-disease state and then to disease state. So, integration of normal samples into dynamic network analysis would reveal some important insights for cancer initiation and occurrence, which will be studied in our future study.

Inferring gene regulatory network from high throughput data is a fundamental question in systems biology yet remains as a challenging task. Analyzing GRNs and their dynamic rewiring are instrumental to identifying master regulator genes or casual drivers in human disease at systems level. The existing methods for inferring GRNs largely depend on time-course gene expression data. However, large samples of time-course data are rarely available in clinical situations since longitudinal surveys are often challenging to conduct. In contrast, cross-sectional studies of a large population based on high-throughput molecular omics data are more prevalent due to their relative feasibility.

The lack of temporal information in clinical transcriptomic data leads to a major challenge for inferring GRN and its translation to precision medicine. Therefore, how to decode the temporal information underlying disease progression is an interesting question and may be key to address the above challenge. The sample similarity-based approach in evolution and genetics studies ^[29], for instance, phylogenetic trees based on microarray data ^[30] and genetic linkage maps based on genetic markers ^[31], have shown developed and shown success in recovering evolutionary dynamics. Importantly, some studies have shown that it is possible to measure cancer progression by using gene expression profiles to estimate “distances” between tumors ^[32]. Another study used graph-based minimum spanning tree to recover biological progression from microarray samples ^[33], which can be applied for reconstructing temporal trajectory of differentiation, development, cell cycle or disease progression.

In summary, in this study, we analyzed the dynamic changes of GRNs along different stages of CRC by using an integrated approach. We combined several clinical transcriptomic datasets of CRC and reconstructed GRNs of CRC across different stages by integrating gene selection, mutual information-based network sparcification and a previous DCM algorithm. Our study demonstrated that the dynamic rewiring of gene regulations could be inferred by reconstructing dynamic GRNs based on clinical transcriptomic profiling. Our study revealed several key genes underlying CRC development, which may generate new hypothesis for further experimental test. We anticipate new methods in future studies to quantitatively estimate disease progression based on the clinical cross-sectional transcriptomic data, which may be promising for developing dynamical system-based GRN inference approach.

Data availability

The gene expression datasets as well as clinical information of the patients were downloaded from the NCBI GEO database, as described in the main text.

Acknowledgements

X Sun was supported by grants from the National Key R & D Program of China (2018YFC0910500), the National Natural Science Foundation of China (11871070), the Guangdong Basic and Applied Basic Research Foundation (2020B1515020047), Guangdong Key Field R & D Plan (2019B020228001）and the Opening Project of Guangdong Province Key Laboratory of Computational Science at the Sun Yat-sen University (2018003).

Conflict of Interests

The authors declare that they have no conflict of interest.

Author contributions

X.S. designed research; A.D. performed research; A.D. and X.S. analyzed data and wrote the paper.