Inference of the Xenopus tropicalis embryonic regulatory network and spatial gene expression patterns
 Zhenzhen Zheng†^{1, 4, 5, 7},
 Scott Christley†^{6},
 William T Chiu^{2},
 Ira L Blitz^{2},
 Xiaohui Xie^{3},
 Ken WY Cho^{2} and
 Qing Nie^{1, 4, 5}Email author
https://doi.org/10.1186/1752050983
© Zheng et al.; licensee BioMed Central Ltd. 2014
Received: 16 April 2013
Accepted: 19 December 2013
Published: 8 January 2014
Abstract
Background
During embryogenesis, signaling molecules produced by one cell population direct gene regulatory changes in neighboring cells and influence their developmental fates and spatial organization. One of the earliest events in the development of the vertebrate embryo is the establishment of three germ layers, consisting of the ectoderm, mesoderm and endoderm. Attempts to measure gene expression in vivo in different germ layers and cell types are typically complicated by the heterogeneity of cell types within biological samples (i.e., embryos), as the responses of individual cell types are intermingled into an aggregate observation of heterogeneous cell types. Here, we propose a novel method to elucidate gene regulatory circuits from these aggregate measurements in embryos of the frog Xenopus tropicalis using gene network inference algorithms and then test the ability of the inferred networks to predict spatial gene expression patterns.
Results
We use two inference models with different underlying assumptions that incorporate existing network information, an ODE model for steadystate data and a Markov model for time series data, and contrast the performance of the two models. We apply our method to both control and knockdown embryos at multiple time points to reconstruct the core mesoderm and endoderm regulatory circuits. Those inferred networks are then used in combination with known dorsalventral spatial expression patterns of a subset of genes to predict spatial expression patterns for other genes. Both models are able to predict spatial expression patterns for some of the core mesoderm and endoderm genes, but interestingly of different gene subsets, suggesting that neither model is sufficient to recapitulate all of the spatial patterns, yet they are complementary for the patterns that they do capture.
Conclusion
The presented methodology of gene network inference combined with spatial pattern prediction provides an additional layer of validation to elucidate the regulatory circuits controlling the spatialtemporal dynamics in embryonic development.
Background
Detailed gene regulatory networks (GRNs) in a number of invertebrate species have provided an unprecedented global overview of the genetic program controlling development in sea urchin, Drosophila, and C. elegans[1–4] and have revealed a number of important and conserved regulatory cassettes employed in a diversity of developmental contexts [5]. While generation of such networks will be also extremely valuable in understanding the mechanisms governing cell fate specification in vertebrate systems, similar work in vertebrates is challenging as the number of cell types, genome organization and genes involved in regulating the biological processes are significantly more complex.
In all triploblastic metazoans, establishment of the primary germ layers (endoderm, mesoderm and ectoderm) occurs early, during blastula and gastrula stages. In the Xenopus blastula the presumptive germ layers are arranged along the vegetalanimal axis with endoderm arising from the vegetal cells, mesoderm in an equatorial ring and the ectoderm on the top overlying the blastocoel cavity. This simple spatial arrangement in developing embryos, taken together with a low complexity in terms of numbers of different cell types and the ease in manipulating gene expression, makes the amphibian Xenopus ideally suited to study GRNs in early vertebrate development.
Xenopus developmental biologists have spent nearly 20 years in generating a prototype GRN describing the mesendoderm [6, 7]. Despite this effort, these GRN diagrams are very incomplete and provide only a limited preview of the in vivo condition. New alternative approaches are urgently needed to generate GRNs that incorporate more genes and have predictive features. In this paper, we present a novel method to elucidate gene regulatory circuits from aggregate gene expression measurements in embryos of the frog Xenopus tropicalis using gene network inference algorithms and then test the ability of the inferred networks to predict spatial gene expression patterns.
The primary methodologies for gene network inference include probabilistic graphical models [8–11], informationtheoretic approaches [12, 13], ordinary differential equations (ODEs) (among which include linear ODEs for steadystate data [14–17], linear ODEs for time series data [15, 16, 18–21] and nonlinear ODEs for time series data that adopt heuristic search strategies [22–26]) and linear regression models [11, 20, 27, 28]. There are numerous reviews of these methods and other approaches [29–34].
In this work, we examine gene expression profile changes of hundreds of genes at several developmental stages after lossoffunction analyses. We then employ two inference models with different underlying assumptions, a linear ODE model for steadystate data and a linear Markov model for time series data, to elucidate the core dorsal mesoderm and endoderm regulatory circuits. Both models incorporate sparseness control on the network connections and prior network information, and they can be solved with the same optimization framework. Using one inferred network in combination with known dorsalventral expression pattern images of a subset of genes, we define an optimization problem to predict spatial patterns for all genes in the network. The spatial pattern prediction provides an additional layer of validation for the regulatory circuits controlling the spatialtemporal dynamics in embryonic development.
where x _{ i }(t) is the concentration of mRNA for gene i measured at time t, dx _{ i }(t)/dt is the rate of change for the mRNA concentration of gene i, and p is the number of genes. Each function F _{ i } represents all of the various processes and factors that affect the amount of mRNA for gene i. Previously, we presented a linear steadystate ODE model for gene network inference that incorporates regularization terms for sparseness and prior network information [17]. We showed that inclusion of prior knowledge about the network structure in the inference process increased performance, that incorrect connections in network structure knowledge did not hurt performance, and that a mixture of correct and incorrect connections given as prior knowledge performed better than giving no prior network information.
We employ our steadystate ODE model to gene expression data from the Xenopus embryo. Since the linear steadystate ODE model assumes that observations are made when the experimental system is at a steadystate equilibrium, the model cannot directly incorporate temporal dynamics for the multiple developmental stages present in our data. One technique to account for such dynamics is to approximate the derivatives for the variables (i. e., dx _{ i }(t)/dt), but this approximation can be inaccurate for the long time intervals typical in biological data. An alternative approach proposed by Linde et al. [21, 35] is to consider a firstorder Markov model where the gene expression at time k is a linear function of its regulators at the previous time k–1, i.e., ${x}_{i}^{k}={\displaystyle \sum _{j=1}^{p}{W}_{\mathit{ij}}{x}_{j}^{k1}}$, and W is the linear gene interaction matrix. However, this model suffers from an issue typical of gene network inference models, which is that the number of genes is greater than the number of experimental observations. Therefore, the system is underdetermined and the model tends to produce a dense gene network that overfits the data. Linde et al. utilize a heuristic search strategy to produce sparse networks, however it is not integrated into the optimization problem and thus it is hard to gauge the effectiveness of the heuristic [21, 35]. Various regularization techniques, which are integrated into the optimization problem, have been introduced to prevent overfitting and to perform variable selection including ridge regression [36], LASSO [37–39], and elastic net [40]. Ridge regression tends to achieve better prediction performance through a biasvariance tradeoff among all the variables, while LASSO specifically enforces sparseness by excluding poor predictor variables, and elastic net combines the two techniques. In our prior work, we applied LASSO in our linear steadystate ODE model to produce a parsimonious regulatory network that is optimal as tested by crossvalidation, and we showed how the LASSO regularization operator could be extended to incorporate prior network information [17]. In this paper, we extend the Markov model to include regularization terms that enforce sparseness of the inferred gene network and allow incorporation of prior network information. We apply the model to simulated data from test networks and present results on the model’s ability to recover the network from differing number of observations and mixtures of correct and incorrect connections provided as prior network information. We apply both models to the aggregate gene expression data of the heterogeneous cell types in the Xenopus embryo, and then compare the ability of each model to recover the core regulatory circuits.
Advances in bioimaging and image analysis are allowing gene expression data to be mapped and studied within a spatial context for organisms and tissue [41–45]. This has led to the recognition and challenge of using spatial gene expression data to reconstruct the regulatory circuits responsible for those spatial patterns, such as in a recent case study of reverse engineering the wellstudied gap gene network responsible for segmentation in the embryo of D. melanogaster[46–49]. Our research is the first attempt to our knowledge to apply similar techniques for Xenopus. One of the challenges is quantifying gene expression from spatial pattern images [50], however we take a simpler approach by categorizing the spatial pattern based upon the assessment of a biological expert. Given a set of spatial gene expression image obtained from Xenbase [51], we transform the expression along the dorsalventral axis of the embryo into a onedimensional representation. We then define an optimization problem that takes an inferred gene network, either from the steadystate ODE or Markov model, and a subset of spatial data to predict the spatial patterns for the remaining genes. We characterize the performance for each model in their ability to predict the spatial expression for genes with known patterns, and we discuss hypothesized spatial patterns for genes where no such data exists. Our approach suggests that a single modeling method is not sufficient to capture all aspects of spatial gene expressions, and the differences in the underlying assumptions for each model may provide insights about the spatialtemporal dynamics in embryonic development.
Results and discussion
Simulation results
We generated a set of time series simulation data to test the Markov model. Five random networks containing p = 10 nodes with 2–3 uniform randomly selected incoming edges were generated for a total of exactly 25 edges in the network; each edge had a weight drawn from the normal distribution N (0, 100). A large variance was used to avoid generating simulation data with big values, described in more detail below. Each network was verified to be nonsingular.
where k = 3, 2, 1, x ^{ k } represents the k th timepoint data in one observation and y ^{ k } represents the k th timepoint ideal data without noise in one observation. Here, two kinds of noise were added: the intrinsic noise (e.g., stochastic fluctuations in the underlying biological process) was drawn from N (0,0.1) and the extrinsic noise (e.g., measurement errors) was drawn from N (0,0.3). The weight of each edge in the randomly generated networks was drawn from N (0, 100). The large variance can help avoid generating data with large values that can skew the inference process and produce numerical errors. Using a small variance tends to generate large values in W ^{–1}, thus each time point will produce increasingly larger values for y ^{ k } and x ^{ k }. The generated random networks and time series data used in producing the simulation results are provided in Additional file 1.
For our experiments that utilize existing network information, we provide a Boolean matrix W ^{0}, where an entry ${W}_{\mathit{ij}}^{0}=0$ indicates a directed interaction from gene j to gene i, while ${W}_{\mathit{ij}}^{0}=1$ for all other edges.
Using leaveoneout crossvalidation, we find the values for the regularization parameters, α (sparsity) and β (prior network), for each gene that minimizes the crossvalidation error. A proportional error is calculated to measure the algorithm’s performance throughout this section. Since we introduce noise into the simulated data, the crossvalidation error will vary with the number of observations. Therefore, in each simulation run we divide the minimal crossvalidation error by the minimal leastsquares error obtained using linear regression without any regularization terms. This normalizes the error relative to the minimal possible error achievable through linear regression. Then we take an average across all the random networks to produce the final proportional error.
Inference error decreases as the number of observations increases
Providing valid edges as prior network information increases performance
Providing invalid edges as prior network information does not affect performance
Consistent performance is maintained with a mixture of valid and invalid edges in prior network information
Comparison of the ODE model and the Markov model
Based on the above observations, three common conclusions can be obtained from the ODE model (Eq. 1.2) and the Markov model (Eq. 1.3): (1) the proportional errors decrease as the number of observations increases; (2) providing invalid edges alone does not affect the prediction performance; (3) providing valid edges is generally helpful to improve the performance especially when only a few observations are available. The difference is that the ODE simulation data [17] is separated into two groups and their noise is drawn from N (0,0.3) and N (0,0.1) respectively, while all the Markov simulation data contains both the above noise simultaneously which generates larger noise in the above simulations and weakens the effect of prior valid edges.
Inference of Xenopus tropicalis embryonic regulatory network
The 46 prior connections and the connections inferred from the linear ODE model and Markov model
Prior information  Linear ODE model  Linear Markov model 

sox17a regulates hnf1b  Inferred  Inferred 
sox17a regulates foxa4a  Inferred  Inferred 
sox17a regulates foxa1  Inferred  Inferred 
sox17a regulates foxa2  Inferred  Inferred 
sox17a regulates gata4  Inferred  
sox17a regulates gata5  Inferred  
sox17a regulates gata6  Inferred  Inferred 
sox17a regulates bix1.2  Inferred  
sox17b.1 regulates hnf1b  Inferred  Inferred 
sox17b.1 regulates foxa4a  Inferred  
sox17b.1 regulates foxa1  Inferred  Inferred 
sox17b.1 regulates foxa2  Inferred  Inferred 
sox17b.1 regulates gata4  Inferred  
sox17b.1 regulates gata5  Inferred  Inferred 
sox17b.1 regulates gata6  Inferred  Inferred 
sox17b.1 regulates bix1.2  Inferred  Inferred 
sox7 regulates sox17a  
sox7 regulates sox17b.1  
gata4 regulates sox17a  
gata4 regulates sox17b.1  Inferred  
gata5 regulates sox17a  Inferred  
gata5 regulates sox17b.1  Inferred  Inferred 
gata6 regulates sox17a  Inferred  
gata6 regulates sox17b.1  
bix1.2 regulates sox17a  Inferred  Inferred 
bix1.2 regulates sox17b.1  Inferred  
vegt regulates mix1  Inferred  
vegt regulates mixer  Inferred  
vegt regulates sox17a  Inferred  
vegt regulates sox17b.1  Inferred  Inferred 
foxh1 regulates otx2  Inferred  Inferred 
foxh1 regulates lhx1  Inferred  Inferred 
foxh1 regulates mix1  Inferred  Inferred 
foxh1 regulates mixer  Inferred  Inferred 
foxh1 regulates bix1.2  Inferred  Inferred 
foxh1 regulates t  
foxh1 regulates ventx2.2  Inferred  
foxh1 regulates sox17a  Inferred  Inferred 
foxh1 regulates sox17b.1  Inferred  
foxh1 regulates frzb  Inferred  Inferred 
foxh1 regulates gsc  Inferred  
foxh1 regulates hhex  Inferred  Inferred 
foxh1 regulates msx1  Inferred  Inferred 
ventx2.2 regulates ventx1.2  Inferred  Inferred 
ventx1.2 regulates myf5  Inferred  Inferred 
t regulates myf5  Inferred 
The inferred network from the ODE model contains 694 edges and 34 out of 46 (pvalue = 0.006) prior connections are correctly inferred. The full list of inferred connections from the ODE model are in Additional file 1: Table S1. The inferred network from the Markov model contains 410 edges and 32 (pvalue = 0) edges are consistent with the prior information. The full list of inferred connections from the Markov model are in Additional file 1: Table S2. Details for pvalue calculations are provided in Additional file 1. 25 connections are shared among the connections in the inferred ODE network and Markov network (details are in Table 1). The Markov network is more sparse than the ODE network, e.g., the average degree of all nodes is 11.39 for the Markov network and 19.28 for the ODE network. There are 172 common connections among all the connections in both networks.
The crossvalidation procedure calculates values for the regularization parameters, α (sparsity) and β (prior network) such that the crossvalidation error is minimized. However, the resultant network does not necessarily contain all of the connections provided in the prior network as with our study where 34 (ODE) and 32 (Markov) of the 46 prior networks connections are in the inferred network. This can occur for a number of reasons, for example: 1) the prior network information may be incorrect and thus excluded, 2) the experimental data may lack a discriminatory signal that the algorithm can use to infer the connection, 3) the crossvalidation error may be too stringent by excluding connections with minimal support, or 4) the nonlinear dynamics of the prior network connection may not be sufficiently captured by the linear model. If crossvalidation for the sparsity control parameter α in (Eqs. 1.2 and 1.3) is not used, α can be varied to produce more or less prior network connections. In the ODE model, by setting α = 0.006 then 42 prior network connections are obtained, while α = 0.0 provides all 46 prior connections. Likewise for the Markov model, setting α = 0.0117 provides 35 prior connections, α = 0.008 provides 40 prior connections, and α = 0.0 provides all 46 prior connections. However, for such cases the crossvalidation errors are not as good as the one obtained through the learning algorithm on the sparsity parameter. For example, the crossvalidation errors are 0.5179, 0.5162 and 0.5129 respectively in the above three settings of α for the Markov model, while the optimal crossvalidation error we obtained was 0.5014. Furthermore, increasing the value of the sparsity control parameter may decrease the number of prior connections by enforcing more sparsity and eliminating connections that are least consistent with the experimental data. Therefore, the number of prior network connections within the inferred network should not be considered as a strict measure of the accuracy of the algorithm, instead it is a relative indication of the information provided within the experimental data that is consistent with the prior network, while taking into account the tradeoff of generalization versus overfitting by the inference algorithm.
Inference of the core dorsal endoderm circuit
Both inference models recovered the core circuitry controlling dorsal endoderm specification including direct regulation of hnf1β, foxa1, foxa2, foxa4a, gata5, gata6 and bix1 by sox17; direct or indirect regulation of gata4 by sox17; as well as the direct or indirect regulatory feedback of gata46 and bix1 onto sox17. Both models predicted vegt regulation of sox17. However, the two models predict the regulation of two sox17 genes, sox17a and sox17b, which are paralogs, to be different. The ODE model usually does not differentiate regulatory action of sox17a and sox17b: both are regulators of hnf1β, foxa1, foxa2, gata4, gata5, gata6 and bix1. The Markov model infers that both sox17a and sox17b are regulators of hnf1b, foxa1, foxa2, foxa4a and gata6. However, the Markov model sometimes splits the regulatory action of sox17a and sox17b: sox17b is a direct regulator of gata5 and bix1 while sox17a is not. Given the differing model assumptions with the ODE model assuming steadystate and the Markov model assuming temporal change, the predictions could suggest that sox17a and sox17b have different temporal actions in the context of the feedback loop with gata46 and bix1[52], even though sox17a and sox17b are similar in their expression and activity [53]. One hypothesis is that sox17b is the primary driver of temporal change for the feedback loop, while sox17a stabilizes those changes. Better understanding of the other factors involved in the feedback loop could help resolve this difference.
Inference of the core dorsal mesoderm circuit
Both inference models recovered the core circuitry controlling dorsal mesoderm specification with foxh1 being a direct regulator of mix1, mixer, lhx1, bix1, otx2, sox17, frzb, msx1 and hhex. Both models predicted ventx2 regulation of ventx1 and ventx1 regulation of myf5 in the ventrolateral mesoderm. The Markov model also predicted foxh1 regulation of gsc in the dorsal mesoderm, and ventx2.2 in the ventrolateral mesoderm. Only the ODE model was able to predict Vegt regulation of mix1 and mixer. Both inference models recovered some core dorsal mesoderm circuit with slightly different gene sets.
Inference of Xenopus tropicalis embryonic spatial gene expression
Predicted gene spatial patterns based on the Markov spatial prediction model for 36 genes from averages with 1000 random initial patterns
Gene  Predicted pattern  Prior pattern 

‘bix1.2’  ‘d’  ‘b’ 
‘bmp4’  ‘v’  ‘v’ 
‘ctnnb1’  ‘v’  ‘d’ 
‘foxa1’  ‘u’  ‘m’ 
‘foxa2’  ‘v’  ‘m’ 
‘foxa4a’  ‘d’  ‘b’ 
‘foxh1’  ‘u’  ‘u’ 
‘foxh1.2’  ‘v’  [ ] 
‘frzb’  ‘d’  ‘d’ 
‘gata4’  ‘d’  ‘m’ 
‘gata5’  ‘u’  ‘m’ 
‘gata6’  ‘v’  ‘m’ 
‘gsc’  ‘d’  ‘d’ 
‘hhex’  ‘u’  ‘d’ 
‘hnf1b’  ‘u’  [ ] 
‘lhx1’  ‘u’  ‘d’ 
‘mespb’  ‘v’  [ ] 
‘mix1’  ‘u’  ‘u’ 
‘mixer’  ‘u’  ‘m’ 
‘msx1’  ‘v’  ‘v’ 
‘myc’  ‘v’  [ ] 
‘myf5’  ‘u’  [ ] 
‘myod1’  ‘u’  [ ] 
‘nodal3’  ‘u’  ‘d’ 
‘nodal6’  ‘u’  ‘m’ 
‘otx2’  ‘d’  ‘d’ 
‘sox17a’  ‘u’  ‘m’ 
‘sox17b.1’  ‘u’  ‘m’ 
‘sox21’  ‘v’  [ ] 
‘Sox7’  ‘v’  [ ] 
‘t’  ‘v’  ‘b’ 
‘vegt’  ‘u’  ‘u’ 
‘ventx1.2’  ‘v’  ‘v’ 
‘ventx2.2’  ‘v’  ‘v’ 
‘wnt11’  ‘b’  ‘b’ 
‘xbp1’  ‘u’  ‘b’ 
Predicted spatial patterns for genes with unknown patterns based on the ODE spatial prediction model and 28 prior gene patterns
Gene  Predicted pattern 

‘foxh1.2’  ‘v’ 
‘hnf1b’  ‘m’ 
‘mespb’  ‘b’ 
‘myc’  ‘m’ 
‘myf5’  ‘d’ 
‘myod1’  ‘b’ 
‘sox21’  ‘b’ 
‘sox7’  ‘d’ 
In the Markov spatial prediction model, we supposed the known spatial gene expression patterns were observed at a particular time point t 2 and they were evolved from a few gene expression patterns at a time point t 1. It was found that two gene expression patterns (e.g., bmp4 and ctnnb1) at t 1 could result in 11 (i.e., 39.3%) correctly predicted gene expression patterns at t 2 using the inferred network from the Markov model for the 28 genes. The correctly predicted genes were: bmp4, foxh1.1, frzb, gsc, hhex, mix1, msx1, otx2, vegt, ventx1.2, ventx2.2. Expression patterns of 11 genes (i.e., 39.3%) at t 2 were correctly predicted when using the same two genes at t1 and the 36gene network inferred from the Markov model for the 36 genes. The prediction results are listed in Table 2.
The ODE and Markov spatial prediction models were applied to steadystate patterns and time series patterns, respectively. The models were able to predict the spatial patterns for some of the key genes involved in mesendoderm specification including vegt, sox17, gata4, ventx and gsc when only provided a subnetwork of genes. However, the majority of the predictions between the two models do not overlap and each model predicts a slightly different subset of the core circuitry. Despite only having a single time point of experimental images and using a predefined spatial pattern for two genes at an early time point, the Markov model has analogous performance to the ODE model in terms of predicting spatial expression patterns. The correct prediction of different gene subsets of the core mesoderm and endoderm regulatory circuitry by the two models may be suggestive of different underlying spatial dynamics for those genes.
Conclusions
There is an increasing need to integrate the approaches that unravel the complicated networks of gene regulatory processes and the works that focus on the spatialtemporal multicellular phenomena of pattern formation and morphogenesis. Currently, the networkcentric studies produce volumes of regulatory interactions typically with little regard to how these networks specify cellular fate in the context of spatial patterns of gene expression. While developmentcentric studies focus on small sets of genes, they require laborintensive approaches, and do not fully embed those genes within the larger regulatory network. Our study represents an initial attempt to integrate these disparate approaches into a single methodology based on biological gene perturbations combined with constraints from spatial modeling. With such an approach, one can make more meaningful predictions for spatial patterns and developmental programs, constrained by the observed complex regulatory networks and in response to changes in gene expression that can be tested experimentally.
We have applied two gene regulatory network inference models with different underlying assumptions to Nanostring experimental data from heterogeneous cell populations from the Xenopus embryo. One inference model is an ODE model that assumes steadystate data, and we have previously developed an optimization framework for this model that incorporates prior network information. The other inference model is a Markov model for time series data. We have shown that the Markov model fits within our optimization framework, and extended the model so that prior network information and sparseness constraints can be incorporated directly into the optimization task. We have tested the extended Markov model on simulated network data and showed that existing network information improves performance and can perform well even when some of the existing network information is partially incomplete. In this regard, we have recently obtained ChIPseq data for sox17 and RNAseq after sox17 MO knockdown [Cho, Blitz and Zorn, unpublished results]. Based on this observation we were able to confirm 4 out of 7 sox17 connections predicted by the steadystate ODE model, and 3 out of 7 connections predicated by the Markov model. Some of these confirmed interactions include newly predicted connections.
Both inference models were able to recover the core circuitry for controlling dorsal endoderm specification and dorsal mesoderm specification. Differences in the model predictions suggest different dynamics that may be related to the underlying assumptions for each model. For the dorsal endoderm circuitry, the ODE model usually does not differentiate regulatory action of sox17a and sox17b, while the Markov model sometimes splits their regulatory action and places more connections for sox17b. This suggests that even though sox17a and sox17b are similar in their expression and activity, they may play different roles in their temporal dynamics in their feedback loop with gata46 and bix1. A putative hypothesis is that sox17b is the primary driver of temporal change for the feedback loop, while sox17a stabilizes those changes. For the dorsal mesoderm circuitry, both inference models recovered some core dorsal mesoderm circuit with slightly different gene sets.
Recent experimental results have provided the opportunity to compare our predictions. A T/T2 double knockdown was performed by microinjection of sequencespecific morpholino antisense oligonucleotides, and RNAseq data of perturbed embryos was obtained at stage 32 [54]. The differential expressions of the selected targets were retrieved from the analyzed datasets for the 36 genes used in this study (see Additional file 7). For comparison purposes, we consider genes with >1.5 fold change to be directly or indirectly regulated by t where a negative log fold change indicates positive regulation, while a positive log fold change indicates negative regulation. The results indicate a total of 18 genes regulated by gene t (Additional file 7). Among them, foxa4a, gsc, mespb, myf5, mix1, bix1.2, myod1 (i.e., totally 7 genes) are positively regulated by t, and bmp4, sox21, vegt, ventx2.2, ventx1.2, msx1, wnt11, foxh1.1, nodal6, mixer, nodal3 (i.e., totally 11 genes) are negatively regulated by t. Using the forward ODE and Markov models (described in Methods), the forward ODE model predicts the positively regulated genes are foxa4a, myf5 and bix1.2; and predicts the negatively regulated genes are sox21, vegt, msx1, wnt11, foxh1.1, nodal6 and mixer. The forward Markov model predicts the positively regulated genes are foxa4a, gsc, myf5, mix1 and myod1; and predicts the negatively regulated genes are bmp4, vegt, foxh1.1, nodal6 and mixer.
Given the two inferred regulatory networks from the ODE model and the Markov model, we additionally constrained these networks by using them to predict spatial gene expression in the Xenopus embryo. Both models were able to predict the spatial patterns for some of the key genes involved in mesoendoderm specification. Interestingly, each model tended to correctly predict a different subset of genes suggesting that those genes are playing different roles in the spatialtemporal dynamics.
The spatial prediction model is dependent upon the provided inferred network for how well it can predict spatial patterns. Similar to the inferred network, the number of prior predicted spatial patterns should not be directly interpreted as the accuracy of the algorithm, because there are numerous reasons why not all of the prior patterns were predicted. A number of (notinclusive) reasons include 1) incomplete or incorrect connections in the inferred network, 2) incorrect or coarse pattern classification for the biological spatial images, or 3) the nonlinear spatiotemporal dynamics are not accurately captured in our 1dimensional abstract model. However, the fact that a statistically significant portion of the prior patterns is predicted suggests that the spatial prediction algorithm is effectively utilizing the information it is given. Given a more accurate inferred network, the spatial predictions should improve.
One limitation of the spatial prediction with regards to the Markov model is that we did not have access to temporal spatial gene expression images. Instead, we took the approach of assuming initial spatial patterns for a small set of wellstudied genes for an earlier time point, then tested the model ability to predict spatial patterns at a later time point. Despite being given only limited initial data, the Markov model was able to correctly predict 39.3% of the gene expression patterns and thus suggesting that the model is accurately capturing some aspects of the temporal dynamics involved in early Xenopus development. In the future, in situ images of gene expression patterns at multiple embryo stages could improve the prediction capability of the Markov model.
Methods
Morpholino knockdown and NanoString analysis of gene expression
Synchronously developing Xenopus tropicalis embryos were obtained by in vitro fertilization, dejellied in pH 7.6 3% cysteine prepared in 1/9^{th}X MMR, and cultured on agarosecoated plates until they reached the 24cell stage [55]. Embryos were microinjected in 1X MMR with translationblocking morpholino antisense oligonucleotides (MO) that targeted foxh1 (22.5 ng), vegt (22.5 ng), sox17α/β (20 ng), or ctnnb1 (βcatenin) (10 ng) and cultured in 1/9^{th}X MMR at 25°C. Sequences of the MOs used in this study were:

foxh1 (TCATCCTGAGGCTCCGCCCTCTCTA)

vegt (TGTGTTCCTGACAGCAGTTTCTCAT)

ctnnb1 (TTTCAACAGTTTCCAAAGAACCAGG)

sox17a (AGCCACCATCAGGGCTGCTCATGGT)

sox17b1/b2 (AGCCACCATCTGGGCTGCTCATGGT)
Total RNA from three biologically independent samples were prepared at the stages of interest using the acid guanidinium thiocyanate phenol chloroform method [56]. Measurements of transcript abundances were performed using the NanoString platform [57]. In brief, RNA sample (total 100 ng) was hybridized to probe sets at 65°C for a minimum of 18 hours. The hybridized probes were recovered using the NanoString Prep Station, and immediately evaluated using the NanoString nCounter. For each reaction 1155 fields of view were counted. Detailed protocols for NanoString transcript counting followed the manufacturer’s instruction manual (http://www.nanostring.com/).
A steadystate ODE model
where ${x}_{i}^{k}$ is the concentration of mRNA for gene i measured at the k th observation (e.g., time point, sample, wild type expression, knockout expression etc.), $\frac{d{x}_{i}^{k}}{\mathit{dt}}$ is the rate of change for the mRNA concentration of gene i at the k th observation, p is the number of genes, W is the gene interaction matrix which is to be inferred.
where n is the number of observations, α is a positive parameter that enforces the sparsity of the interaction network, ${\u2225W\u2225}_{1}={\displaystyle \sum _{i=1}^{p}{\displaystyle \sum _{j=1}^{p}\left{W}_{\mathit{ij}}\right}},{W}^{0}$ is a Boolean network containing existing network information where ${W}_{\mathit{ij}}^{0}=0$ indicates a directed interaction from gene j to gene i and thus is not penalized while ${W}_{\mathit{ij}}^{0}=1$ for all the other edges and β ≥ 0 indicates the strength of the penalty. ${\overline{x}}^{k}={x}^{k}$ if x ^{ k } is a wildtype data point. If x ^{ k } is a knockout data, we previously put forward one method from a perturbation aspect [17]. Here we suggest another but simple way as follows. Since a knocked out gene does not contribute to the change of other genes’ concentrations, we can let ${\overline{x}}_{j}^{k}={x}_{j}^{k}$ for all j = 1,…, p except ${\overline{x}}_{i}^{k}=0$ if gene i is knocked out at a single observation k. All these notations are used throughout the paper with the same meanings.
A Markov model
where obs _{max} represents the number of observations. Here one observation means one sample (e.g., wildtype data or knockdown data) with a complete time series, k _{max} represents the number of time points in each observation, and all the other notations have the same meanings as before. In particular, ${\overline{x}}^{k}={x}^{k}$ if x ^{ k } is a wildtype data point. If gene i is knocked down at time k in one observation obs we can let ${\left({\overline{x}}_{j}^{k}\right)}^{\mathit{obs}}={\left({x}_{j}^{k}\right)}^{\mathit{obs}}$ for all j = 1,…, p except i and ${\left({\overline{x}}_{i}^{k}\right)}^{\mathit{obs}}=0$. Notice, here we do not require j ≠ i as we did in (Eq. 1.2) since for gene i the concentration at time $k\left(\mathrm{i}.\mathrm{e}.,{\overline{x}}_{i}^{k}\right)$ could contribute to that at time $k+1\left(\mathrm{i}.\mathrm{e}.,{x}_{i}^{k+1}\right)$.
W _{ i } denotes the ith row of the matrix W, D ^{ i } is a p × p zero matrix for i = 1,…, p.
where the matrices D ^{ i } and U are defined in (Eq. 1.8) or (Eq. 1.9).
The optimization problem (Eq. 1.10) can be solved by combining an iterative coordinate descent algorithm for a given pair of parameters (α, β) and a leaveoneout crossvalidation to find the optimal values of (α, β) which provide the minimal crossvalidation error [17, 59]. Here ‘leaveoneout’ means to leave one ‘observation’ out, which implies all the timeseries data in one observation are left out for the Markov model. We perform an exponential search starting from max U _{ ij }, where U is defined in (Eq. 1.8) or (Eq. 1.9), and going down. Since the incoming edges for each gene in the gene network can be considered independent from the incoming edges of other genes, we used a separate α and β for each gene. For each gene and (α, β), zero initial and an accuracy control of 10^{5} are used in the optimization procedure.
Knockout simulations
 (1)
A Forward ODE Model
where W ∈ R ^{ p × p } is the inferred network derived from the ODE model (Eq. 1.2) with zero diagonal elements, x ∈ R ^{ p } represents a vector of p gene expressions.
 (2)
A Forward Markov Model
where W ∈ R ^{ p × p } is the inferred network derived from the Markov model (Eq. 1.3), x ^{ k } ∈ R ^{ p } (k = 1,2,…) is the kth iteration vector of p gene expressions.
We define the initial vector x ^{1} as a random vector drawn from the uniform distribution between 0 and 1. We consider the similar knockout as before. The only difference is to add a superscript k, i.e., keep ${x}_{i}^{k}=0$ (i ∈ K) for all k.
Spatial prediction model
Given a regulatory network, known spatial expression patterns for genes in the network can be used to predict the unknown spatial patterns for the remaining genes in the network. Examples of typical spatial gene expression patterns in Xenopus gastrula stage embryos are shown in Figure 7. For simplicity, we regard the dorsalventral axis as a onedimensional interval, which is partitioned into three regions, i.e., right (dorsal), middle (vegetal) and left (ventral). We took a set of 28 in situ images of the 36 total genes being considered in this study, and categorized the gene expression level for the three regions as either low, medium or high within a given embryo. For computational purposes, we assigned arbitrary values to these levels of 0.1 (low), 0.4 (medium) and 1.0 (high) (the effect of varying these values is provided in Additional file 1). The overall spatial expression pattern could then be assigned a category of dorsal (‘d’), ventral (‘v’), both dorsal and ventral (‘b’), vegetal (‘m’) or uniform (‘u’). The expression values across the three regions are (0.1, 0.4, 1.0) for “dorsal” and (1.0, 0.4, 0.1) for “ventral”. The other categories have multiple possible values with “both” having higher dorsal and ventral values than the middle region (e.g. 1.0, 0.1, 1.0), “vegetal” having a higher value in the middle region than the dorsal and ventral regions (e.g. 0.1, 1.0, 0.1), and “uniform” having the same values across all regions (e.g. 0.4, 0.4, 0.4). The categorized spatial gene expression patterns are shown in Table 2.
where pos represents the three onedimensional regions of the embryo, Index is the set of genes with known patterns, ${y}_{k}^{\mathit{pos}}$ represents the concentration of gene k with known patterns at position pos (i.e. one of 1, 0.4 and 0.1), and W is the inferred network from the steadystate ODE model. This optimization problem (Eq. 1.13) considers the inferred network W as fixed, thus directly using the coefficients and genetogene interaction of W to define the structure of the spatial prediction model. Some of the ${X}_{i}^{\mathit{pos}}$ values are provided as known gene spatial expression patterns, while the remaining ${X}_{i}^{\mathit{pos}}$ values are free variables. The resultant values of the free variables are interpreted as spatial predictions for their associated genes. The optimization problem (Eq. 1.13) satisfies the known gene spatial patterns while simultaneously constraining the model with the genetogene interaction topology from the steadystate equations in the ODE model (Eq. 1.1), so we call (Eq. 1.13) ODE spatial prediction model, which is a quadratic programming and can be solved by the MATLAB function ‘lsqlin’.
where pos represents the three onedimensional regions of the embryo, N is the number of genes, ${y}_{k}^{\mathit{pos},1}$ represents the concentration of gene k with known patterns at position pos (i.e. one of 1, 0.4 and 0.1) at time t 1, ${x}_{k}^{\mathit{pos},\phantom{\rule{0.1em}{0ex}}1}$ and ${x}_{k}^{\mathit{pos},2}$ represent the concentrations of gene k at position pos at time t 1 and t 2, respectively, W is the inferred network from the Markov model, Index1 is the index set of genes observed at time t1. This optimization problem (Eq. 1.14) considers the inferred network W as fixed, thus directly using the coefficients and genetogene interaction of W to define the structure of the spatial prediction model. Since no information is available for the initial patterns of the other genes (i.e., ${x}_{k}^{\mathit{pos},1},k\in \left\{1,\dots ,N\right\}\backslash \mathit{Index}1$), when we solved the optimization problem (1.14) we used random initial patterns (for genes k, k ∈ {1, …, N}\Index 1) for 1000 optimization runs and averaged the values of ${x}_{k}^{\mathit{pos},2}$ as the final prediction. Since $\left\{{x}_{k}^{\mathit{pos},1}\right\}$ are all defined, this model is actually a direct computation for $\left\{{x}_{k}^{\mathit{pos},2}\right\}$ with a nonnegative constraint. We attempted to leave the expression values for the other genes at time t1 to be variables and have the optimization model predict their values. However this results in an underdetermined model with an infinite set of solutions.
 1.
Check if $\frac{\mathit{min}}{\mathit{max}}\ge \mathit{TH}$ or max < 0.1,
 2.
Check if m < l and m < r and $\frac{\mathit{abs}\left(lr\right)}{\mathit{abs}\left(mr\right)}\le 0.5$ and $\frac{\mathit{abs}\left(lr\right)}{\mathit{abs}\left(lm\right)}\le 0.5$,
 3.
Check if m > l and m > r,
 4.
Check if l > m and l > r,
 5.
Check if r > m and r > l,
If yes, output ‘Dorsal (‘d’)’; otherwise, output ‘Not one of the five patterns’.
Notes
Declarations
Acknowledgements
We thank E. Davidson (Caltech) for the use of NanoString. This work was supported by NIH grants 5T32HD06055502 (WTC), RO1HD056219 (KWYC), and P50GM76516 (QN), and NSF grant DMS1161621 (QN).
Authors’ Affiliations
References
 Davidson EH, Levine MS: Properties of developmental gene regulatory networks. Proc Natl Acad Sci USA. 2008, 105 (51): 2006320066. 10.1073/pnas.0806007105.PubMed CentralView ArticlePubMedGoogle Scholar
 Maduro MF: Endomesoderm specification in caenorhabditis elegans and other nematodes. Bioessays. 2006, 28 (10): 10101022. 10.1002/bies.20480.View ArticlePubMedGoogle Scholar
 Christiaen L: The transcription/migration interface in heart precursors of ciona intestinalis. Science. 2008, 320 (5881): 13491352. 10.1126/science.1158170.View ArticlePubMedGoogle Scholar
 Biemar F: Comprehensive identification of drosophila dorsalventral patterning genes using a wholegenome tiling array. Proc Natl Acad Sci USA. 2006, 103 (34): 1276312768. 10.1073/pnas.0604484103.PubMed CentralView ArticlePubMedGoogle Scholar
 Erwin DH, Davidson EH: The evolution of hierarchical gene regulatory networks. Nat Rev Genet. 2009, 10 (2): 141148.View ArticlePubMedGoogle Scholar
 Loose M, Patient R: A genetic regulatory network for xenopus mesendoderm formation. Dev Biol. 2004, 271 (2): 467478. 10.1016/j.ydbio.2004.04.014.View ArticlePubMedGoogle Scholar
 Koide T, Hayata T, Cho KW: Xenopus as a model system to study transcriptional regulatory networks. Proc Natl Acad Sci USA. 2005, 102 (14): 49434948. 10.1073/pnas.0408125102.PubMed CentralView ArticlePubMedGoogle Scholar
 Friedman N: Using bayesian networks to analyze expression data. J Comput Biol. 2000, 7: 601620. 10.1089/106652700750050961.View ArticlePubMedGoogle Scholar
 Friedman N: Inferring cellular networks using probabilistic graphical models. Science. 2004, 303: 799805. 10.1126/science.1094068.View ArticlePubMedGoogle Scholar
 Perrin BE: Gene networks inference using dynamic bayesian networks. Bioinformatics. 2003, 19 (Suppl 2): ii138ii148.View ArticlePubMedGoogle Scholar
 HaibeKains B: Predictive networks: a flexible, open source, web application for integration and analysis of human gene networks. Nucleic Acids Res. 2012, 40 (Database issue): D866D875.PubMed CentralView ArticlePubMedGoogle Scholar
 Faith JJ: Largescale mapping and validation of escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biol. 2007, 5 (1): e810.1371/journal.pbio.0050008.PubMed CentralView ArticlePubMedGoogle Scholar
 Margolin AA: Reverse engineering cellular networks. Nat Protoc. 2006, 1 (2): 662671. 10.1038/nprot.2006.106.View ArticlePubMedGoogle Scholar
 Gardner TS: Inferring genetic networks and identifying compound mode of action via expression profiling. Science. 2003, 301: 102105. 10.1126/science.1081900.View ArticlePubMedGoogle Scholar
 Bonneau R: The inferelator: an algorithm for learning parsimonious regulatory networks from systemsbiology data sets de novo. Genome Biol. 2006, 7 (5): R3610.1186/gb200675r36.PubMed CentralView ArticlePubMedGoogle Scholar
 Gupta R: A computational framework for gene regulatory network inference that combines multiple methods and datasets. BMC Syst Biol. 2011, 5: 5210.1186/17520509552.PubMed CentralView ArticlePubMedGoogle Scholar
 Christley S, Nie Q, Xie X: Incorporating existing network information into gene network inference. PLoS ONE. 2009, 4 (8): e679910.1371/journal.pone.0006799.PubMed CentralView ArticlePubMedGoogle Scholar
 Gustafsson M, Hornquist M, Lombardi A: Constructing and analyzing a largescale genetogene regulatory network–lassoconstrained inference and biological validation. IEEE/ACM Trans Comput Biol Bioinform. 2005, 2 (3): 254261. 10.1109/TCBB.2005.35.View ArticlePubMedGoogle Scholar
 Guthke R: Dynamic network reconstruction from gene expression data applied to immune response during bacterial infection. Bioinformatics. 2005, 21 (8): 16261634. 10.1093/bioinformatics/bti226.View ArticlePubMedGoogle Scholar
 Li Z: Largescale dynamic gene regulatory network inference combining differential equation models with local dynamic Bayesian network analysis. Bioinformatics. 2011, 27 (19): 26862691. 10.1093/bioinformatics/btr454.View ArticlePubMedGoogle Scholar
 Linde J: Regulatory interactions for iron homeostasis in aspergillus fumigatus inferred by a systems biology approach. BMC Syst Biol. 2012, 6: 610.1186/1752050966.PubMed CentralView ArticlePubMedGoogle Scholar
 Yip KY: Improved reconstruction of in silico gene regulatory networks by integrating knockout and perturbation data. PLoS One. 2010, 5 (1): e812110.1371/journal.pone.0008121.PubMed CentralView ArticlePubMedGoogle Scholar
 Vu TT, Vohradsky J: Nonlinear differential equation model for quantification of transcriptional regulation applied to microarray data of saccharomyces cerevisiae. Nucleic Acids Res. 2007, 35 (1): 279287.PubMed CentralView ArticlePubMedGoogle Scholar
 Kimura S: Inference of Ssystem models of genetic networks using a cooperative coevolutionary algorithm. Bioinformatics. 2005, 21 (7): 11541163. 10.1093/bioinformatics/bti071.View ArticlePubMedGoogle Scholar
 Vohradsky J: Neural network model of gene expression. FASEB J. 2001, 15 (3): 846854. 10.1096/fj.000361com.View ArticlePubMedGoogle Scholar
 To CC, Vohradsky J: Measurement variation determines the gene network topology reconstructed from experimental data: a case study of the yeast cyclin network. FASEB J. 2010, 24 (9): 34683478. 10.1096/fj.10160515.View ArticlePubMedGoogle Scholar
 D’Haeseleer P: Linear modeling of mRNA expression levels during CNS development and injury. Pac Symp Biocomput. 1999, 4: 4152.Google Scholar
 Holter NS: Dynamic modeling of gene expression data. Proc Natl Acad Sci USA. 2001, 98 (4): 16931698. 10.1073/pnas.98.4.1693.PubMed CentralView ArticlePubMedGoogle Scholar
 Jong HD: Modeling and simulation of genetic regulatory systems: a literature review. J Comput Biol. 2002, 9: 67103. 10.1089/10665270252833208.View ArticlePubMedGoogle Scholar
 Bansal M: How to infer gene networks from expression profiles. Mol Syst Biol. 2007, 3: 78PubMed CentralView ArticlePubMedGoogle Scholar
 Bonneau R: Learning biological networks: from modules to dynamics. Nat Chem Biol. 2008, 4 (11): 658664. 10.1038/nchembio.122.View ArticlePubMedGoogle Scholar
 De Smet R, Marchal K: Advantages and limitations of current network inference methods. Nat Rev Microbiol. 2010, 8 (10): 717729.PubMedGoogle Scholar
 Hurley D: Gene network inference and visualization tools for biologists: application to new human transcriptome datasets. Nucleic Acids Res. 2012, 40 (6): 23772398. 10.1093/nar/gkr902.PubMed CentralView ArticlePubMedGoogle Scholar
 Marbach D: Wisdom of crowds for robust gene network inference. Nat Methods. 2012, 9 (8): 796804. 10.1038/nmeth.2016.PubMed CentralView ArticlePubMedGoogle Scholar
 Linde J: Regulatory network modelling of iron acquisition by a fungal pathogen in contact with epithelial cells. BMC Syst Biol. 2010, 4: 14810.1186/175205094148.PubMed CentralView ArticlePubMedGoogle Scholar
 Hoerl AE, Kennard RW: Ridge regression: biased estimation for nonorthogonal problems. Technometrics. 2000, 42: 8086. 10.1080/00401706.2000.10485983. URL: http://www.jstor.org/stable/1271436 10.1080/00401706.2000.10485983View ArticleGoogle Scholar
 Efron B: Least angle regression. Ann Stat. 2004, 32 (2): 407451. 10.1214/009053604000000067.View ArticleGoogle Scholar
 Friedman J, Hastie T, Tibshirani R: Sparse inverse covariance estimation with the graphical lasso. Biostat (Oxford, England). 2008, 9 (3): 432441. 10.1093/biostatistics/kxm045.View ArticleGoogle Scholar
 Tibshirani R: Regression shrinkage and selection via the lasso. Journal of the royal statistical society. Series B Methodol. 1996, 58 (1): 267288.Google Scholar
 Zou H, Hastie T: Regularization and variable selection via the elastic net. J Roy Stat Soc B. 2005, 67: 301320. 10.1111/j.14679868.2005.00503.x.View ArticleGoogle Scholar
 Frise E, Hammonds AS, Celniker SE: Systematic imagedriven analysis of the spatial drosophila embryonic expression landscape. Mol Syst Biol. 2010, 6: 345PubMed CentralView ArticlePubMedGoogle Scholar
 Kerwin J: The HUDSEN atlas: a threedimensional (3D) spatial framework for studying gene expression in the developing human brain. J Anat. 2010, 217 (4): 289299. 10.1111/j.14697580.2010.01290.x.PubMed CentralView ArticlePubMedGoogle Scholar
 Pepperkok R, Ellenberg J: Highthroughput fluorescence microscopy for systems biology. Nat Rev Mol Cell Biol. 2006, 7 (9): 690696. 10.1038/nrm1979.View ArticlePubMedGoogle Scholar
 Spencer WC: A spatial and temporal map of C. elegans gene expression. Genome Res. 2011, 21 (2): 325341. 10.1101/gr.114595.110.PubMed CentralView ArticlePubMedGoogle Scholar
 Vermot J, Fraser SE, Liebling M: Fast fluorescence microscopy for imaging the dynamics of embryonic development. HFSP J. 2008, 2 (3): 143155. 10.2976/1.2907579.PubMed CentralView ArticlePubMedGoogle Scholar
 Crombach A: Efficient reverseengineering of a developmental gene regulatory network. PLoS Comput Biol. 2012, 8 (7): e100258910.1371/journal.pcbi.1002589.PubMed CentralView ArticlePubMedGoogle Scholar
 Jaeger J: Dynamic control of positional information in the early drosophila embryo. Nature. 2004, 430 (6997): 368371. 10.1038/nature02678.View ArticlePubMedGoogle Scholar
 Mjolsness E, Sharp DH, Reinitz J: A connectionist model of development. J Theor Biol. 1991, 152 (4): 429453. 10.1016/S00225193(05)803911.View ArticlePubMedGoogle Scholar
 Perkins TJ: Reverse engineering the gap gene network of drosophila melanogaster. PLoS Comput Biol. 2006, 2 (5): e5110.1371/journal.pcbi.0020051.PubMed CentralView ArticlePubMedGoogle Scholar
 Botman D, Kaandorp JA: Spatial gene expression quantification: a tool for analysis of in situ hybridizations in sea anemone nematostella vectensis. BMC Res Notes. 2012, 5: 55510.1186/175605005555.PubMed CentralView ArticlePubMedGoogle Scholar
 Bowes JB: Xenbase: gene expression and improved integration. Nucleic Acids Res. 2010, 38 (Database issue): D607D612.PubMed CentralView ArticlePubMedGoogle Scholar
 Sinner D: Global analysis of the transcriptional network controlling xenopus endoderm formation. Development. 2006, 133 (10): 19551966. 10.1242/dev.02358.View ArticlePubMedGoogle Scholar
 Clements D, Woodland HR: Changes in embryonic cell fate produced by expression of an endodermal transcription factor, Xsox17. Mech Dev. 2000, 99 (1–2): 6570.View ArticlePubMedGoogle Scholar
 Gentsch GE: In vivo Tbox transcription factor profiling reveals joint regulation of embryonic neuromesodermal bipotency. Cell Rep. 2013, 4 (6): 11851196. 10.1016/j.celrep.2013.08.012.PubMed CentralView ArticlePubMedGoogle Scholar
 Khokha MK: Techniques and probes for the study of Xenopus tropicalis development. Dev Dyn. 2002, 225 (4): 499510. 10.1002/dvdy.10184.View ArticlePubMedGoogle Scholar
 Chomczynski P, Sacchi N: Singlestep method of RNA isolation by acid guanidinium thiocyanatephenolchloroform extraction. Anal Biochem. 1987, 162 (1): 156159.View ArticlePubMedGoogle Scholar
 Geiss GK: Direct multiplexed measurement of gene expression with colorcoded probe pairs. Nat Biotechnol. 2008, 26 (3): 317325. 10.1038/nbt1385.View ArticlePubMedGoogle Scholar
 Hecker M: Integrative modeling of transcriptional regulation in response to antirheumatic therapy. BMC Bioinformatics. 2009, 10: 26210.1186/1471210510262.PubMed CentralView ArticlePubMedGoogle Scholar
 Friedman J, Hastie T, Tibshirani R: Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010, 33 (1): 122.PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.