Gene expression is controlled by the combinatorial effects of regulatory factors from different biological subsystems such as general transcription factors (TFs), cellular growth factors and microRNAs. A subsystem’s gene expression may be controlled by its internal regulatory factors, exclusively, or by external subsystems, or by both. It is thus useful to distinguish the degree to which a subsystem is regulated internally or externally–e. g. , how non-conserved, species-specific TFs affect the expression of conserved, cross-species genes during evolution. We developed a computational method (DREISS, dreiss. gerteinlab. org) for analyzing the Dynamics of gene expression driven by Regulatory networks, both External and Internal based on State Space models. Given a subsystem, the “state” and “control” in the model refer to its own (internal) and another subsystem’s (external) gene expression levels. The state at a given time is determined by the state and control at a previous time. Because typical time-series data do not have enough samples to fully estimate the model’s parameters, DREISS uses dimensionality reduction, and identifies canonical temporal expression trajectories (e. g. , degradation, growth and oscillation) representing the regulatory effects emanating from various subsystems. To demonstrate capabilities of DREISS, we study the regulatory effects of evolutionarily conserved vs. divergent TFs across distant species. In particular, we applied DREISS to the time-series gene expression datasets of C. elegans and D. melanogaster during their embryonic development. We analyzed the expression dynamics of the conserved, orthologous genes (orthologs), seeing the degree to which these can be accounted for by orthologous (internal) versus species-specific (external) TFs. We found that between two species, the orthologs have matched, internally driven expression patterns but very different externally driven ones. This is particularly true for genes with evolutionarily ancient functions (e. g. the ribosomal proteins), in contrast to those with more recently evolved functions (e. g. , cell-cell communication). This suggests that despite striking morphological differences, some fundamental embryonic-developmental processes are still controlled by ancient regulatory systems.
Gene regulatory networks systematically control the gene expression dynamics. These networks are highly modular, and consist of various sub-networks. Each sub-network contains a number of regulatory factors representing a subsystem that drives specific gene regulatory functions [1,2]. The subsystems interact with one another, and work together to carry out the entire gene regulatory function. For example, the gene expression in embryogenesis is controlled by the combinatorial effects of various regulatory subsystems composed of complex evolutionary regulatory networks [3]. These regulatory subsystems drive very diverse developmental programs, from the highly conserved (e. g. DNA replication) to the species-specific (e. g. body segmentation). As such the orthologous genes that are evolutionary conserved genes across species can therefore be regulated by both orthologous and species-specific transcription factors (TFs) [4]. The orthologous TFs form an “internal” regulatory network, while the species-specific TFs form an “external” one. Unfortunately, existing experimental gene expression data cannot decouple the expression components that are driven by the different subsystems. Thus, computational methods are required to assess the contribution from each factor or subsystem from the gene expression data. In this study, we propose a novel computational method, DREISS—dynamics of gene expression driven by external and internal regulatory networks based on state space model. Using DREISS, we are able to identify temporal gene expression dynamic patterns for evolutionarily conserved genes during embryonic development, as driven by conserved and species-specific regulatory subsystems. These results advance our current understanding of gene regulatory networks during evolution, as well as the differentiation during development. Developmental gene regulatory networks control gene expression during the developmental processes. These particular regulatory networks have evolved, making it difficult to understand their regulatory mechanisms at the system level. Hence, one typically compares developmental gene expression across species to infer biological activities of developmental gene regulatory networks. For example, embryogenesis provides a platform to study the evolution of gene expression between different species. Recent work has showed that significant biological insight can be gained by cross-species comparisons of the expression profiles during embryogenesis for worms [5], flies [6], frogs [7] and several other vertebrates [8]. It was found that the orthologous genes have minimal temporal expression divergence during the phylotypic stage, a middle phase during the embryonic development across species within the same phylum. These patterns are often characterized as “hourglass” [9]. In addition, the conserved hourglass patterns were observed even within a single species while comparing the developmental gene expression data across distant species, such as worm and fly [10]; i. e. , the expression divergence among evolutionarily conserved genes become minimal during the phylotypic stage in both worm and fly. However, much less is known about how the orthologous genes in each species eventually contribute to their species-specific phenotypes due to the lack of appropriate computational approaches. Thus, we aim to use DREISS to discover the components of the orthologous gene expression during embryonic development driven by species-specific transcription factors. The state-space model has been widely used in engineering [11], and also in biology for the analysis of gene expression dynamics [12–14]. It models the dynamical system output as a function of both the current internal system state and the external input signal. A well-known example in engineering is the vehicle cruise control system where the system state can be the vehicle’s speed. Based on the road conditions, the cruise control requires various fuel amounts in order to keep the desired speed level. In biology, we can look at the transcription factors and microRNAs as internal and respectively external regulatory factors of the protein-coding genes expression (See more internal-external examples in S1 Table). Similarly, the state-space model can be applied for studying the expression of orthologous genes at different developmental stages using information regarding their expression (internal) and species-specific regulatory factors (external) at the current known developmental stage. Unlike earlier studies that calculate the expression correlation between individual genes, the state-space model predicts the temporal causal relationships at the system level; i. e. , the state at a time is determined by the state and external input at the previous time. The earlier work applied the state-space model to study the gene expression dynamics focusing on small-scale systems, and did not explore the analytic dynamic characteristics of the inferred state-space models. The complex and large-scale biological datasets, especially temporal gene expression data, are very noisy, and high dimensional (i. e. , the number of genes is much greater than the number of time samples), thereby preventing an accurate estimation of the state-space model’s parameters. The dimensionality reduction techniques have thus been used to project high-dimensional genes to low-dimensional meta-genes (i. e. , the selected features representing de-noised and systematic expression patterns [1,15,16]) as well as the principal dynamic patterns for those meta-genes [17,18]. Using DREISS, we are able to apply the dimensionality reduction to the gene expression data, and develop an effective state-space model for their meta-genes, and then identify a group of canonical temporal expression trajectories representing the dynamic patterns driven by the effective conserved and species-specific meta-gene regulatory networks according to the model’s analytic characteristics. These dynamic patterns reveal temporal gene expression components that are controlled by conserved or species-specific GRNs. DREISS is a general-purpose tool and can be used to study the gene regulatory effects from any different subsystems for a given group of genes. As an illustration, we applied DREISS to the gene expression data during embryonic development for two model organisms, worm (Caenorhabditis elegans) and fly (Drosophila melanogaster). In both species, we were able to identify the expression patterns of worm-fly orthologs driven by the conserved regulatory network consisting of the worm-fly orthologous TFs (i. e. , the conserved regulatory subsystems between two species), as well as the worm/fly-specific regulatory network consisting of non-orthologous TFs (i. e. , the species-specific regulatory subsystem). Our results reveal that, in addition to executing conserved developmental functions between worm and fly, their orthologous genes are also regulated by species-specific TFs to involve in species-specific developmental processes. In summary, DRIESS provides a framework to analyze both distantly and closely related species allowing for a better understanding of the gene regulatory mechanisms during development.
A gene regulatory network is made up of various subsystems [1,2]. These subsystems work together to execute regulatory functions. Given a group of N1 genes in a subsystem, defined as the internal gene set, Ω, their gene expression levels are not only controlled by internal interactions among Ω, but also affected by the regulatory factors from other subsystems outside Ω. We define an external gene set, Ψ consisting of those external regulatory factors. For example, we consider the worm-fly orthologous genes as internal set Ω. The worm-fly orthologous TFs from internal set Ω are the internal regulatory factors, and non-orthologous TFs such as worm- or fly- specific TFs are the external regulatory factors. Both the internal and external regulatory factors control gene expressions in dynamic ways (i. e. , their regulatory signals at the current time will affect gene expressions at subsequent times). Thus, the regulatory mechanisms for gene expressions form a control system. In this study, we used a state-space model (defined by linear first-order difference equations, Fig 2A) to formulate temporal gene expression dynamics for internal set Ω (comprising N1 genes) with external regulation from external set Ψ (comprising N2 genes) at time points 1,2, …, T as follows: Xt+1=AXt+BUt (1), where the vector Xt∈RN1×1, the “state”, includes N1 gene expression levels at time t in Ω, and the vector Ut∈RN2×1, the “input or control”, includes N2 gene expression levels at time t in Ψ. The system matrix A∈RN1×N1 captures internal causal interactions among genes in Ω (i. e. , the ith, jth element of A, Aij describes the contribution from the jth gene expression at time t to the ith gene expression at the next time t+1), which instantiates a gene regulatory network. The control matrix B∈RN1×N2 captures external causal regulations from the genes in Ψ to genes in Ω (i. e. , the ith, jth element of B, Bij describes the contribution from the jth gene expression in Ψ at time t to the ith gene expression in Ω at the next time t+1). R represents the real number domain. According to the state space model (1), the gene expression dynamics in Ω is determined by the system matrix A and the control matrix B. In particular, based on Eq 1, the state Xt can be expanded as follows: Xt=AXt−1+BUt−1=A (AXt−2+BUt−2) +BUt−1=A2Xt−2+ABUt−2+BUt−1=A3Xt−3+A2BUt−3+ABUt−2+BUt−1=⋯=At−1X1+At−2BU1+At−3BU2+⋯+ABUt−2+BUt−1=At−1X1⏟XtINT+∑k=1t−2AkBUt−1−k⏟XtINTER+BUt−1⏟XtEXT (2), where XtINT=At−1X1 is defined as the expression vector of the gene components driven only internally by genes in Ω. The rest terms ∑k=1t−2AkBUt−1−k+BUt−1 captures the expression expression vector of the gene components in Ω affected externally by the genes in Ψ. In particular, XtEXT=BUt−1 represents the expression vector of gene components in Ω driven purely by the genes in Ψ since it only involves B and U, and XtINTER=∑k=1t−2AkBUt−1−k captures the expression vector of gene components in Ω driven by the interactions between internal and external groups for involving A, B and U. In this paper, we mainly focus on the purely internal dynamics. As for the external-related terms, we should emphasize that any subdivision of the rest of the terms ∑k=1t−2AkBUt−1−k+BUt−1 is completely arbitrary. That is, although we subdivided it into a purely external term and an interaction term here, one could subdivide it multiple ways. That is, given a particular type of subdivision, each of the subdivided terms sums up a group of terms from AkBUt−1−k, k = 0,1, 2, …, t-2. For example, one can look at ∑k=2t−2AkBUt−1−k+ (ABUt−2+BUt−1), where ABUt−2 + BUt−1 shows the contribution from the inputs up to two time steps back to Xt. The temporal gene expression experiments normally have limited time samples (for example, there may only be a dozen time points), which are far less than the time samples needed to estimate the large matrices A and B when internal and external groups, Ω and Ψ are composed of hundreds or thousands of genes. One way to deal with lack of time samples is dimensionality reduction. Thus, we project high dimensional temporal gene expressions to much lower dimensional meta-gene expression levels using a dimensionality reduction technique (Fig 2B). Those meta-gene expression levels should capture original gene expression patterns, such as the ones having the greatest degree of co-variation. We calculate the meta-gene expression levels as follows: X˜t=WX*Xt; U˜t=WU*Ut (3), where X˜t∈RM1×1, the “meta-gene state” at time t, includes M1 (<< N1 and <T) meta-gene expression levels; i. e. , the first M1 elements of the tth row of the matrix whose columns are right-singular vectors of the matrix [X1 X2 ⋯ XT] in Ω by the singular value decomposition (SVD) [19]; the vector U˜t∈RM2×1, the “meta-gene input or control” at time t, includes M2 (<< N2 and <T) meta-gene expression levels; i. e. , the first M2 elements of the tth row of the matrix whose columns are right-singular vectors from SVD of the matrix [U1 U2 ⋯ UT] in Ψ; WX∈RN1×M1 is the linear projection matrix of SVD from M1 meta-gene expression space to N1 gene expression space in Ω, WU∈RN2×M2 is the linear projection matrix of SVD from M2 meta-gene expression space to N2 gene expression space in Ψ, and (.) * is a pseudo-inverse operation; i. e. , W*W = I, where I is the identity matrix. Next, we obtain the effective state-space model for meta-genes using linear projections WX and WU between genes and meta-genes as follows (Fig 2C). By replacing (1) using (3), we obtain that WXX˜t+1=AWXX˜t+BWUU˜t. (4) After multiplying the pseudo-inverse of WX, WX*∈RM1×N1 s. t. WX*WX=I where I is an identity matrix, at both sides of (4), we have that X˜t+1=WX*AWX⏟A˜X˜t+WX*BWU⏟B˜U˜t=A˜X˜t+B˜U˜t (5), where the effective meta-gene system matrix A˜=WX*AWX∈RM1×M1 captures internal causal interactions among meta-genes in Ω (i. e. , an element of A˜, A˜ij describes the contribution from the jth meta-gene expression at time t to ith meta-gene expression at time t+1), and the effective control matrix B˜=WX*BWU∈RM1×M2 captures external causal regulations from meta-genes of Ψ to meta-genes of Ω (i. e. , the ith, jth element of B˜, B˜ij describes the contribution from the jth meta-gene expression in Ψ at time t to ith meta-gene expression in Ω at time t+1). Eq 5 describes the effective state space model for the meta-genes of Ω, whose expression dynamics is determined by A˜ and B˜. Because the meta-gene dimension, M1 (M2) is less than T, and much less than N1 (N2), we can estimate A˜ and B˜ as follows. We rewrite Eq 5 as a matrix product on the right side: X˜t+1=A˜X˜t+B˜U˜t=[A˜B˜][X˜tU˜t]. (6) By applying Eq 6 to time points, 2,3, …, T, we then obtain that [X˜2X˜3⋯X˜T]⏟Ζ=[A˜B˜][X˜1X˜2⋯X˜T−1U˜1U˜2⋯U˜T−1]⏟Υ (7), where Ζ∈RM1× (T−1) and Υ∈R (M1+M2) × (T−1). Because of dimension reduction, Υ has more columns than rows so that it has right pseudo-inverse. Thus, the effective internal system matrix A˜ and external control matrix B˜ can be estimated by: [A˜B˜]=ΖΥ* (8), where Υ*∈R (T−1) × (M1+M2) is the right pseudo-inverse of Υ; i. e. , ΥΥ* = I, with M1<N1, M2<N2, M1+M2<T, t = 1,2, …, T. It is worth noting that if we do not reduce the dimensionality, and obtain Eq 7 from Eq 5, then Υ will have much more rows than columns so that it doesn’t have right pseudo-inverse; i. e. , there doesn’t exist a matrix Υ* such that ΥΥ* is a full-rank identify matrix. In addition, the condition of M1+M2<T also makes ΥΥ* be a full-rank identify matrix. The analytic solution to a general first-order linear matrix difference equation [20], Qt+1 = CQt is Qt = CtQ0 = (HEH-1) tQ0 = HEtH-1Q0 = HEtS, where the columns of the matrix H are eigenvectors of C, the diagonal elements of the diagonal matrix E are eigenvalues of C such that CH = HE, and the vector S = H-1Q0. Then, if we rewrite Qt by a linear combination of the time exponential of eigenvalues of C, we have that Qt=HEtS=∑i=1mcαitsiHi=∑i=1mcαitKi, where mc is the total number of eigenvalues of C, αi is the ith eigenvalue of C, si is the ith element of S, Hi is the ith eigenvector of C (i. e. , the ith column of H), and Ki = siHi is the coefficient vector of Qt over the tth time exponential of αi. By Eq 5, the matrix Ã determines the meta-gene states components whose expression dynamics are internally controlled by the meta-genes of Ω. As Eq 2, we define the expression of the meta-gene components driven only internally by themselves in Ω at time t as X˜tINT, an M1-dimensional vector; i. e. , their expression at two adjacent time points have X˜t+1INT=A˜X˜tINT∈RM1×1. According to the above analytic solution, it can be a linear combination of M1 dynamic patterns determined by the eigenvalues of the effective system matrix A˜ as follows: X˜tINT=∑p=1M1λptK˜p; i. e. , the internally driven component of ith meta-gene’s expression across all time points, [X˜1INT (i) X˜2INT (i) …X˜TINT (i) ]=∑p=1M1K˜p (i) [λp1λp2…λpT]⏟pthiPDP (9), where λp and K˜p∈CM1×1 are the pth eigenvalue of A˜ and its coefficient vector from the analytic solution, which determines the pth dynamic pattern driven by effective internal regulations, defined as the pth internal principal dynamic pattern (iPDP) = [λp1λp2…λpT], in which λpt represents the tth power of λp, and Ξ (i) represents ith element of the vector Ξ. C represents the complex number domain. If an eigenvalue λ is complex when A˜ is asymmetric, then its conjugate λ¯ is also an eigenvalue, so we sum its iPDP and its conjugate eigenvalue, λ¯’s iPDP, as a unified iPDP with real elements equal to [λp1+λ¯p1λp2+λ¯p2…λpT+λ¯pT]. The internal principal dynamic patterns (iPDPs) represent canonical temporal expression trajectories, which can be either increasing, or damped oscillation and so on depending on iPDP’s eigenvalues (Fig 3). The iPDPs can be ordered by sorting their eigenvalues. Also by Eqs 2 and 5, the expression of the meta-gene states components driven purely by the external group Ψ at time t is defined as X˜tEXT, an M1-dimensional vector, and its expression dynamics is determined by the equation X˜t+1EXT=B˜U˜t∈RM1×1; i. e. , the externally driven components of meta-gene states at two adjacent time points. In particular, the externally driven component of ith internal meta-gene’s expression across time points: [X˜2EXT (i) X˜3EXT (i) …X˜TEXT (i) ]=∑q=1M2B˜i, q[U˜1 (q) U˜2 (q) …U˜T−1 (q) ]⏟qthePDP (10), where X˜tEXT (i) and U˜t (q) are ith and qth elements of X˜tEXT and U˜t, respectively with t = 1,2, …, T, the vector [U˜1 (q) U˜2 (q) …U˜T−1 (q) ]⏟ is defined as qth external principal dynamic pattern (ePDP), and B˜i, q is the element of B˜ at ith row and qth column, which is also the coefficient of the externally driven component of ith internal meta-gene’s expression over qth ePDP. Based on Eq 2, the expression of the meta-gene components driven by the interactions between internal and external meta-genes is given by X˜tINTER=∑k=1t−2A˜kB˜U˜t−1−k. In this paper, we focus on the purely driven internal patterns (i. e. , iPDPs) and compare them across different biological systems. Because genes and meta-genes have linear relationships in terms of their expression levels as described in Eq 2, the components of gene expression levels in Ω driven by internal regulations, XtINT∈RN1×1 can be also expressed as linear combinations of M1 iPDPs: XtINT=WXX˜tINT=∑p=1M1λptWXK˜p⏟Cp=∑p=1M1λptCp; i. e. , the internally driven component of ith gene’s expression across all time points, [X1INT (i) X2INT (i) …XTINT (i) ]=∑p=1M1Cp (i) [λp1λp2…λpT]⏟pthiPDP (11), where Cp=WXK˜p∈CM1×1 represents the gene coefficient vector for pth iPDP. Similarly, the gene expression components driven by external genes in Ψ, XtEXT∈RN1×1 can be also expressed as linear combinations of M2 ePDPs: XtEXT=WXX˜tEXT=WXB˜⏟DU˜t=DU˜t; i. e. , the externally driven component of ith gene’s expression across all time points, [X2EXT (i) X3EXT (i) …XTEXT (i) ]=∑q=1M2Di, q[U˜1 (q) U˜2 (q) …U˜T−1 (q) ]⏟qthePDP (12), where XtEXT (i) is ith element of XtEXT with t = 1,2, …, T, and Di, q is the element of D=WXB˜ at ith row and qth column, which is also the coefficient of the externally driven component of ith gene’s expression over qth ePDP.
DREISS enables us to compare expression dynamic patterns between two or more temporal gene expression datasets even though they have different numbers of samples, as well as differences in the times at which those samples were collected. For example, we can apply DREISS to two different datasets of the same group of genes, and identify both the common (similar) and the specific (different) dynamic patterns driven by internal regulations captured by the eigenvalues of the effective system matrices between the two datasets. In this paper, we apply DREISS to 3,153 one-to-one orthologous genes between worm (Caenorhabditis elegans) and fly (Drosophila melanogaster) as internal group, Ω to study their expression dynamics during embryonic development [10]. We refer to species-specific TFs as external regulations; i. e. , external group Ψ. We found that worm-fly orthologs have similar internal dynamic patterns, which may be mainly driven by conserved TFs, but have very different external dynamic patterns driven by species-specific TFs between worm and fly embryonic developmental stages. The data is summarized as follows. We define internal group Ω as 3,153 one-to-one orthologous genes between worm and fly during embryonic development, and external group Ψ as all the species-specific TFs (509 worm-specific TFs, 442 fly-specific TFs) [21,22]. We used their temporal gene expression levels (as measured by the RPKM values in RNA-seq) during embryonic development from the modENCODE project [10]. The worm embryonic development dataset includes T = 25 time stages at 0,0. 5,1, 1. 5, …, 12 hours, and the fly dataset includes T = 12 time stages at 0,2, 4, …, 22 hours, but t = 1,2, . . , 25 for worm and t = 1,2, …, 12 for fly are used in this paper, representing the relative time points for the entire embryonic development processes. Because M1+ M2<T in Eq 8, we choose M1 = M2 = 5 meta-genes for fly (T = 12), and find that five meta-genes of Ω and five meta-genes of Ψ capture ~98% of the co-variation of orthologous gene expressions and fly-specific TF gene expressions, respectively. In order to compare worm and fly, we also choose M1 = M2 = 5 meta-genes for worm, which capture ~98% of the co-variation of orthologous gene expressions and worm-specific TF gene expressions. We find that the meta-gene canonical temporal expression trajectories driven by conserved regulatory networks (i. e. , internal principal dynamic patterns, iPDPs) include four major patterns in both the worm and fly embryonic developmental process by order of eigenvalues: 1) a late highly varied pattern; 2) an early fast decaying pattern; 3) a slowly increasing pattern; and 4) an oscillating pattern (Fig 4A); i. e. , the pattern of canonical trajectories (VL, D, I, O) in Fig 3. In contrast to the observed iPDP similarities, we find that worm and fly have very different external principal dynamic patterns (ePDPs) (Fig 4B); i. e. , the expression dynamic patterns driven by species-specific TFs. The principal dynamic patterns driven by the worm-specific regulatory network; i. e. , worm ePDPs, include a varied pattern that decreases until the middle stage and then increases, an increasing pattern, a varied pattern with a peak entering middle stage, a pattern that varies early and then increases during the embryonic development, and a cosine-like oscillating pattern with roughly two periods during the embryonic development. The fly ePDPs, however, have a varied pattern with low expression at the early stage, a sine-like oscillating pattern with roughly one period during the embryonic development, an increasing pattern, another sine-lie oscillating pattern with roughly two periods during the embryonic development, and a varied pattern that is like damped oscillation. In addition, we checked the sensitivity of iPDPs to small perturbations to internal/external regulatory networks by the leave-one-out method; i. e. , we removed one gene in the internal/external group, ran DREISS, and obtained the ordered iPDP eigenvalues for the remaining genes. We repeated the leave-one-out method for all genes, and finally found the ranges in which iPDP eigenvalues vary shown as error bars in S1 Fig. We can see that the iPDP eigenvalues almost stay at the same values (small error bars) for both worm and fly, which implies that the principal dynamic patterns of worm-fly orthologous genes driven by their conserved regulatory network are robust to small changes. The above results suggest that the conserved regulatory networks from orthologous meta-genes between worm and fly have similar effects to orthologous meta-genes, given their similar iPDPs (i. e. , both have four patterns, as described above). The species-specific regulatory networks from species-specific meta-genes (i. e. , worm-specific or fly specific TFs) have effects that differ from the orthologous meta-genes for their different ePDPs. In addition, the expression dynamic patterns driven by the interactions between internal orthologous genes and external species-specific TFs are also different between worm and fly (S2 Fig). In both worm and fly, we observe the similar four types of internally driven canonical temporal expression trajectories; i. e. , four matched internal principal dynamic patterns (iPDPs) (Fig 4A). Thus, we are interested in seeing how individual orthologous genes relate to those dynamic patterns. We find that the worm-fly orthologous genes have correlated coefficients over each of the four iPDPs. Based on Eq 10, we can obtain the coefficients of orthologous genes for each iPDP. We find that their coefficients are significantly correlated between worm and fly iPDPs with a similar pattern (Fig 5): r = 0. 33 (p<2. 2e-16) for the highly varied pattern at late embryonic development stages (first iPDP), r = 0. 66 (p<2. 2e-16) for the fast decaying pattern at early embryonic development stages (second iPDP), r = 0. 67 (p<2. 2e-16) for the slowly increasing pattern during embryonic development (third iPDP), and r = 0. 73 (p<2. 2e-16) for the oscillation pattern during embryonic development (forth iPDP), where r represents Spearman correlation of iPDP coefficients of 3,153 orthologous genes between worm and fly. This implies that, not only do the orthologous meta-genes have similar internal (conserved) regulatory effects (i. e. , similar iPDPs), but the worm-fly orthologous genes also have similar internally-driven expression dynamics as resulted from their significantly correlated coefficients for iPDPs. The ePDPs between worm and fly generally do not show a high degree of matching similarity, but the worm ePDP No. 2, and the fly ePDPs No. 3 are roughly representing the growing patterns. We find that orthologous gene correlation coefficients between these ePDP patterns are very small (Spearman correlation r = -0. 22 of the orthologous gene coefficients of worm ePDP No. 2 and fly ePDP No. 3). The ribosome produces proteins, which is an ancient process and conserved across worm and fly, organisms separated by almost a billion years of evolution. The ribosomal genes are highly expressed during embryogenesis, since intensive cell division and migration require a large amount of proteins to be synthesized. We collected 195 ribosome-related genes based on the GO annotations. We ranked the coefficients of orthologous genes for each iPDP and ePDP in ascending order, and compared the rank values of iPDP and ePDP coefficients of ribosomal genes. We found that their average ranks of iPDP coefficients are significantly larger than ePDP ones in both worm (t-test p<2. 2e-16) and fly (t-test p<2. 6e-13) as shown in Fig 6. This means that the ribosomal gene expression is significantly more influenced by the conserved regulatory network than by the species-specific regulatory network, which is consistent with ribosomal genes having conserved functions during embryonic development. The orthologous genes related to signal transduction for cell-cell communication (a significantly more recent evolutionary adaptation relative to the ribosome) exhibit the opposite trend. We found that 320 signaling genes from GO annotations have significantly larger average rank values of ePDP coefficients than iPDP ones in both worm (t-test p<5. 6e-11) and fly (t-test p<8. 3e-4), as shown in Fig 6. This result implies that the signaling gene expression is significantly more driven by the species-specific regulatory network than by the conserved regulatory network, which is consistent with the signaling genes being commonly associated with species-specific functions, such as body plan establishment and cell differentiation. We next turn to the biological meaning of individual canonical temporal expression trajectory for iPDPs and ePDPs. For the fast-decaying pattern (2nd iPDP), we find that the DNA replication is significantly enriched in Top 300 (~10%) orthologous genes that have the most negative coefficients for this pattern, in both worm (p<1. 6e-8) and fly (p<4. 5e-6). The GO enrichment analysis was performed using DAVID [23]. The very negative coefficients for the fast decaying pattern mean high positive coefficients for a fast-growing pattern (vertically flipped 2nd iPDPs of worm and fly represent a fast-growing pattern), showing a drastic increase at the beginning of embryogenesis, then remain flat during the late embryogenesis (red curves in Fig 7). Most of the cell division of embryogenesis in both worm and fly happens approximately within the first 300 minutes. Then, the cell elongation and migration start to dominate the development [24,25]. The mRNA abundance of the genes involved in DNA replication may change accordingly. This is well reflected by the second iPDP. Interestingly, the original expression patterns of those top orthologous genes actually do not have fast-growing patterns (black curves in Fig 7), probably because of the combined effects of both conserved and species-specific GRN. Maternal mRNAs, which are pre-loaded before fertilization, may also mask the fast growing pattern of DNA replication genes. This pattern could only be observed after we separated the effect of two types of TFs using DREISS. In addition, we did not find any enrichment of DNA replication in top genes of other iPDPs (p>0. 05). Therefore, the fast-growing iPDP patterns identified by our method reveal conserved regulation on the elementary cellular process of both species (i. e. DNA replication). Besides a fast growing pattern driven by conserved worm-fly orthologous TFs, we also identified a fast growing pattern driven by non-conserved TFs for the two species. The Top 300 orthologous genes (~10%) with the fast-growing worm ePDP (ePDP No. 2) (i. e. , driven by species-specific regulatory networks) are enriched in ‘proteasome’ (p<9. 8e-16). Protein degradation is not only a key process in apoptosis, but also throughout the entire course of development [26,27]. For example, eliminating proteins that are no longer needed is a vital process during embryo development; e. g. , the maternal proteins need to be cleaned as the embryogenesis proceeds). Previous reports also showed that different species usually have different maternal mRNA in the oocyte, which indicates that species-specific strategies might be utilized to regulate the protein degradation process [28]. In this study, after separating the effect of conserved and non-conserved regulatory networks, we observed that the protein degradation is significantly enriched in the genes majorly driven by species-specific TFs in worms. In contrast, the Top 300 orthologous genes with fast growing fly ePDP3 are enriched in ‘mitotic cell cycle’ (p<3. 5e-29), ‘translation’ (p<1e-30) and ‘mitochondrion’ (p<7. 7e-20). Those enriched function related to energy generation is probably indicative of the large energy requirement during fly embryogenesis [29], which did not provide the evolutionary conservation of this energy-related gene regulation. Our result reveals that the fly genes associated with respiration are more up-regulated by fly-specific TFs relative to conserved TFs, and that this up-regulation evolved after the separation of worm and fly. Besides the fast-growing pattern driven by species-specific TFs, we also observed some other interesting patterns. For example, worm ePDP3 displays a dramatic peak about 5 hours after fertilization. Among the Top 300 worm orthologous genes of this pattern, genes involved in synaptic transmission (p<5. 6e-9) and cell-cell signaling (p<1e-7) are over-represented, suggesting that they are transiently activated in this stage of embryogenesis by worm-specific TFs. This observation indicates the gene regulatory network for these genes have evolved after the speciation. We applied DREISS to another example (also see supplement) about cancer. We are also interested to identify the gene expression dynamic patterns driven by conserved and human-specific regulatory networks during breast cancer cell cycle. Thus, we applied DREISS to a time-series gene expression data for human estrogen-responsive breast cancer cell line (ZR-75. 1) before and after hormonal stimulation, which has 12 time points covering a complete mitotic cell cycle (0–32 hours) of hormonal stimulated cells [30]. The internal group, Ω is defined as a set of cross-species conserved human genes (i. e. , 1132 worm-fly-human orthologs including 150 orthologous TFs), and the external group, Ψ consists of 1870 human-specific TFs. As shown in S3 Fig, the internally driven principal dynamic patterns (iPDPs) of conserved human genes include an oscillation trajectory whose period is roughly equal to a full cell cycle (iPDP No. 4), but the externally driven patterns (ePDPs No. 2–4) oscillates more frequently than internal one, which suggests that though the evolutionarily conserved TFs regulate the normal cell cycle, the human specific TFs potentially drive the abnormal cycling behaviors of conserved gene expression responding to the hormonal stimulation.
In this paper, we presented a novel computational method, DREISS, which decomposes time-series expression data of a group of genes into the components driven by the regulatory network inside the group (internal regulatory subsystem), and the components driven by the external regulatory network consisting of regulators outside the group (external regulatory subsystem). DREISS is a general-purpose tool that can be used to study the gene regulatory effects of any interested biological subsystems such as protein-coding transcription factors, micro-RNAs, epigenetic factors and so on. As an illustration, we applied DREISS to the time-series gene expression datasets for worm and fly embryonic developments from the modENCODE project [10], and compared the worm-fly orthologous gene expression dynamic patterns driven by the conserved regulatory network (i. e. , regulation effects from orthologous TFs), with the patterns driven by the species-specific regulatory networks (i. e. , regulation effects from worm or fly specific TFs). We found that the conserved TFs drive similar genomic functions, but non-conserved TFs drive species-specific functions of orthologous genes between worm and fly, implying that, in addition to having ancient conserved functions, orthologous genes have been regulated by evolutionarily younger GRNs to execute species-specific functions during the evolution. This work can be easily extended to study the regulatory effects from orthologous TFs and species-specific TFs to species-specific genes. For example, one can find the expression dynamic patterns of worm/fly specific genes driven by specific TFs, and identify the genes with strong patterns associated with worm/fly specific functions, such as body formations. To the best of our knowledge, DREISS is the first method to reveal how the evolution of GRNs affects gene expression during embryogenesis. We emphasize that DREISS is a general-purpose method (a free downloadable R tool available from github. com/gersteinlab/dreiss). Users can define the internal group (Ω) and external group (Ψ) according to their interests. For example, if users want to identify the protein-coding expression patterns driven by miRNAs, they can define miRNAs as an external group and protein-coding genes as an internal group. Additionally, DREISS can be applied to more than two datasets, such as comparing worm, fly and human embryonic stem cell developmental data, and finding their conserved and specific developmental expression patterns. The expression patterns driven by human-specific regulatory factors will potentially help us understand human-specific developmental processes along with the associated human genes. Due to the limited time samples in gene expression datasets, DREISS uses the simple linear state space model (i. e. the first order linear invariant difference equation) to model the temporal gene expression dynamics, and identify principal temporal dynamic patterns. This model assumes that the gene regulatory networks controlling temporal gene expression dynamics does not change across the entire biological process such as (A, B) in Eq 1. Thus, based on the analytic analysis, the principal dynamic patterns (PDPs) must follow a small set of canonical temporal trajectories (Fig 3). With the rapidly increasing gene expression data, we can extend DREISS to more advanced models such as switched and hybrid system models, non-linear models [31], which will allow us to study the gene regulatory networks are time varying, and potentially find the more temporal gene expression patterns capturing the more complex gene regulatory activities.
The dynamics of a biological system can be controlled by its own internal mechanisms and external perturbations. To gain intuition on this, we may draw a comparison with a mass hanging from a spring. The mass will move naturally by itself but its dynamics is also affected by one’s pulling it. That is, the dynamics of the mass is governed by the effect of the external perturbations superimposed on the internal mechanism of the spring (i. e. Hooke’s law). Similarly, given a group of genes, their temporal gene expression dynamics can be controlled by both transcription factors inside the group and external regulatory factors. Therefore, it is useful to identify the expression dynamics that are exclusively controlled by internal or external factors and compare them across various systems. While state-space models have been widely used to decouple the internal and external effects in physical systems, such as the mass and spring, typical biological systems do not have enough time samples to infer all the model’s parameters, and applications of state-space models were not very effective in these instances. Hence, we developed a general-purpose computational method by integrating state-space models and dimensionality reduction to identify temporal gene expression patterns driven by internal and external regulatory networks. We applied our method to the embryonic developmental datasets in the worm and fly (and also in a human cancer context). We successfully identified the temporal expression dynamics of cross-species conserved genes that were driven by conserved and species-specific regulatory networks.
genetic oscillators
developmental biology
mathematics
linear algebra
algebra
gene regulatory networks
embryogenesis
gene expression
genetics
gene regulation
biology and life sciences
physical sciences
computational biology
evolutionary biology
eigenvalues
evolutionary genetics
evolutionary developmental biology
DREISS: Using State-Space Models to Infer the Dynamics of Gene Expression Driven by External and Internal Regulatory Networks