PEAK Predicts Gene Regulatory Network Linkages during Sea Urchin Development with High Sensitivity from Gene Expression Data Alone

Gene regulatory network (GRN) inference can now take advantage of powerful machine learning algorithms to predict the entire landscape of gene-to-gene interactions with the potential to complement traditional experimental methods in building gene networks. However, the dynamical nature of embryonic development -- representing the time-dependent interactions between thousands of transcription factors, signaling molecules, and effector genes -- is one of the most challenging arenas for GRN prediction. In this work, we show that successful GRN predictions for developmental systems from gene expression data alone can be obtained with the Priors Enriched Absent Knowledge (PEAK) network inference algorithm. PEAK is a noise-robust method that models gene expression dynamics via ordinary differential equations and selects the best network based on information-theoretic criteria coupled with the machine learning algorithm Elastic net. We test our GRN prediction methodology using two gene expression data sets for the purple sea urchin (S. purpuratus) and cross-check our results against existing GRN models that have been constructed and validated by over 30 years of experimental results. Our results found a remarkably high degree of sensitivity in identifying known gene interactions in the network (maximum 76.32%). We also generated 838 novel predictions for interactions that have not yet been described, which provide a resource for researchers to use to further complete the sea urchin GRN.


Abstract Background
Gene regulatory network (GRN) inference can now take advantage of powerful machine learning algorithms to predict the entire landscape of gene-to-gene interactions with the potential to complement traditional experimental methods in building gene networks. However, the dynamical nature of embryonic development --representing the time-dependent interactions between thousands of transcription factors, signaling molecules, and effector genes --is one of the most challenging arenas for GRN prediction.

Results
In this work, we show that successful GRN predictions for developmental systems from gene expression data alone can be obtained with the Priors Enriched Absent Knowledge (PEAK) network inference algorithm. PEAK is a noise-robust method that models gene expression dynamics via ordinary differential equations and selects the best network based on information-theoretic criteria coupled with the machine learning algorithm Elastic net. We test our GRN prediction methodology using two gene expression data sets for the purple sea urchin (S. purpuratus) and cross-check our results against existing GRN models that have been constructed and validated by over 30 years of experimental results. Our results found a remarkably high degree of sensitivity in identifying known gene interactions in the network (maximum 76.32%). We also generated 838 novel predictions for interactions that have not yet been described, which provide a resource for researchers to use to further complete the sea urchin GRN.

Conclusions
GRN predictions that match known gene interactions can be produced using gene expression data alone from developmental time series experiments.

Background
In the 1950s and 1960s, it began to be clear to researchers that gene products representing the sum of the entire genome were not present in every cell and, therefore, gene expression was being regulated (Mirsky 1951;Britten and Davidson 1969). The exciting discovery of the role of transcription factors in the process of gene regulation set off a new era of exploration of the roles of genes in phenotypes, disease, evolution, and embryonic development (Jacob and Monod 1961). Gene regulation can be organized and modeled as a hierarchical network, a Gene Regulatory Network (GRN), as rst proposed by Eric Davidson (Arnone and Davidson 1997;Davidson et al. 2002a). GRN models are now routinely used to follow the causal links between regulatory genes that lead to cell fate decisions or cell activities. GRN models are also used to create hypotheses about the function of players in a regulatory program.
Beginning with the construction of the sea urchin endomesoderm GRN (Davidson et al. 2002a(Davidson et al. , 2002b, many GRNs in different model systems have been compiled through extensive experimental perturbations often involving a combination of knockdown techniques combined with visualization of changes in gene expression. Accuracy of the GRN is improved when further experiments con rm cisregulatory interactions at the level of transcription factor binding sites. However, there is now a pressing need to facilitate GRN modeling using computational prediction tools to help ll in the gaps in existing GRNs and to help create new GRN models with testable interaction predictions, especially in systems with limited resources. Various GRN prediction algorithms from gene expression data have been proposed (Delgado and Gómez-Vela 2019), and several have been compared and tested through the DREAM consortium (Marbach et al. 2012a). The overall accuracy of these methods using gene expression data alone is usually below 50% even with a consensus of several computational methods. The use of additional prior information, such as transcription factor binding data, Gene Ontology (GO) annotation, functional association, proteinprotein interactions (PPIs), and experimentally veri ed relations, was shown to improve prediction accuracy (Altarawy et al. 2017). However, this information is not always available for every new model system that could bene t from GRN inference.
For our analysis, we use the Priors Enriched Absent Knowledge (PEAK) network inference algorithm to reconstruct GRN interactions from gene expression data (Altarawy et al. 2017). PEAK uses differential equations, CLR (Context Likelihood of Relatedness), and the machine learning method Elastic net to predict the likely interactions between target genes and transcription factors. PEAK consists of two phases, a coarse-grained phase and a ne-grained phase, to predict a GRN. In the coarse-grained phase, potential regulators for each gene are extracted using mixed context likelihood of relatedness (mixed CLR). In the ne-grained phase, two modi ed versions of Elastic net are used to re ne the predictions and to integrate curated or noisy prior knowledge, when available. Prior knowledge about the network can also be added if available; however, none has been used in this study.
We chose two sea urchin embryonic GRNs governing endomesoderm and ectoderm speci cation as test networks, because they are widely regarded as some of the most well-supported developmental GRN models with many cis-regulatory interactions veri ed at the base-pair level. We also obtained two sea urchin gene expression data sets to use as input. The California purple sea urchin, Strongylocentrotus purpuratus, is a marine invertebrate. The sea urchin is a member of the phylum Echinodermata, which, along with the Hemichordata, are the closest known sister groups to the Chordates, the phylum to which humans belong. The genome of S. purpuratus was sequenced in 2006, which produced an estimated gene count of 23,000 (Sea Urchin Genome Sequencing Consortium et al. 2006). The sea urchin develops rapidly from a single-cell (the fertilized egg) to a late gastrula embryo in 48 hours at 14 o C and then into a prism-shaped larvae by 72 hours. As an indirect developing animal, the sea urchin larvae must undergo metamorphosis to achieve its nal adult structure. Over the last 20 years, GRN models describing the regulation of embryonic development of S. purpuratus have been built by experimentation and collaboration of sea urchin researchers around the world. The GRN models describing development are divided into the network describing the ectodermal tissue layer, which will give rise to the nervous system and outer tissues of the larvae, and the network controlling the endomesodermal tissue layer, which will give rise to the gut and the larval skeleton. The most recent versions of the two GRN models are hosted at BioTapestry, accessible at http://grns.biotapestry.org/SpEndomes/ for the endomesoderm network and at http://grns.biotapestry.org/SpEcto/ for the ectoderm network (Longabaugh et al. 2005).
The goal of our approach is to employ PEAK to predict gene regulatory interactions from gene expression data alone. For many emerging model systems crucial for understanding evolution, gene expression data from developmental time series is readily available but, due to limited resources, these systems lack extensive prior information regarding transcription factor binding sites, ChIPseq data, and proteomics. Therefore, these emerging model systems have a demonstrated need for GRN inference to guide future experiments on regulatory interactions during embryonic development (Tulin et al. 2013;Fernandez-Valverde et al. 2018).
Previous large scale assessment of network inference methods aimed at predicting gene networks using gene expression data alone have a maximum of 50% precision using 800 microarray experiments (Marbach et al. 2012a) and much less accuracy using 300 experiments. Our method of using the PEAK algorithm on gene expression data from developmental time series data alone was able to achieve a maximum of 76.32% sensitivity using 98 experiments.

Results
Differential gene expression analysis A complete GRN model would need to include the entire set of expressed transcripts whose expression is controlled by the regulatory apparatus. To begin to de ne the regulatory gene set to input into the PEAK machine learning algorithm, we started with the sea urchin RNAseq transcript data set covering 0-72hpf (Tu et al. 2014). The RNAseq transcript data set was then ltered to identify the set of transcripts that are differentially expressed during embryonic development and are regulative in nature.
We used three programs (NOISeq, EdgeR, and GFold) to determine the set of differentially expressed genes (DEGs) with the parameters described in the methods section. There was variation in the number of DEGs as determined by NOISeq, EdgeR, and GFold on the transcriptome RNAseq data. We compared the overlap of genes above the threshold identi ed by each method to obtain a core set of DEGs (Fig. 1). There were 20,422 total unique differentially expressed transcripts determined by 5 methods (2 different thresholds for both GFold and NOISeq, and one threshold for EdgeR). Only 10,627 genes were consistently speci ed as differentially expressed no matter which method was used. We found that the result from NOISeq (λ 1 = 0.9) contained the most overlap and the least difference with results from the other methods while maintaining a more selective number of genes determined to be differentially expressed (Table 1). Only .01% of the genes determined by NOISeq (λ 1 = 0.9) are unique to that method; in contrast, more than 5% of the genes from the result of EdgeR are unique to EdgeR. All the genes determined by GFold (λ 4 = ∓ 1.5) overlapped with the results from other methods, but the GFold DEG set was missing 2,755 genes that were identi ed as differentially expressed by the other 4 methods.

Gene Ontology Filter
The NOISeq lter q_value > 0.9 reduced the 21,092 total transcripts to a set of 15,496 differentially expressed transcripts. However, the set of DEGs identi ed by NOISeq (λ 1 = 0.9) is still too large to be effectively used as input into the PEAK prediction algorithm. Therefore, we applied a second lter to the gene set to achieve a more appropriate number of regulatory genes. The second lter was a Gene Ontology (GO) lter for genes related to transcription and signaling. We used the custom GO annotation generated by the authors of the transcriptome (Tu et al. 2012). The GO lter identi ed 1,038 transcripts that are regulatory in nature, of which 544 were also identi ed as differentially expressed (Additional le 1). The 544 transcripts represent 504 unique gene models annotated by a single SPU_ID.

PEAK analysis
We next applied the PEAK GRN prediction algorithm on the ltered gene set of 504 DEGs as determined by NOISeq that are also identi ed by the GO lter. We speci ed to PEAK the set of 258 Transcription Factors (TFs) used during embryonic development according to a compilation of genes speci ed as TFs To evaluate the performance of PEAK, we repeated the PEAK analysis with data relating only to genes present in the ground truth sea urchin GRNs using the same procedure and parameters as with the set of all DGE genes that t the GO lter. In this analysis, we used three measures to perform the evaluation. Sensitivity represents the proportion of our predicted gene interactions that hit the corresponding ground truth GRNs, indicating the percentage of gene interactions that are correctly identi ed by the PEAK algorithm among the total known ones. Because there are still numerous undiscovered gene interactions and other complexities of GRNs such as transient binding, multiple upstream regulators, and indirect interactions (Van den Broeck et al. 2020), new predictions designates the proportion of new gene interactions predicted by the algorithm which is not yet known in the ground truth GRN. And the miss rate represents the proportion of known gene interactions not discovered by the algorithm.
We evaluated the gene set from the ectoderm GRN and the endomesoderm GRN separately and compared the results as summarized in Table 2. For the known ectoderm GRN, there are 36 genes out of 39 that are identi ed as DEGs in the transcriptome RNAseq data. When the gene expression data from the transcriptome RNAseq data set is used as input for PEAK, the algorithm successfully predicted 45 of 78 gene-to-gene interactions present in the ground truth GRN, yielding 57.69% sensitivity ( Table 2). PEAK failed to predict 33 known connections but provided 555 new possible ones. The sea urchin endomesoderm GRN is currently a larger network model, with 58 genes and 146 edges. We found 57 of those 59 genes were identi ed as being differentially expressed by NOISEq. When TFs were speci ed, 83 of the 142 connections were correctly predicted for a sensitivity of 58.45% (Table 2). To achieve a higher sensitivity, we considered what limitations might be present in the approach. A limitation of the transcriptome data when applied to the ground truth GRN is that only 10 timepoints were sampled in total and only 5 of those timepoints covered the period of time during development which the ground truth GRN models describe. Therefore, we sought out a second data set with more timepoints at closer sampling intervals. For the second data set, we chose a high-density embryonic data set where 161 genes critical to early development were sampled once an hour for 48 hours in duplicate and gene expression was quanti ed with Nanostring technology (Materna et al. 2010). Because there are genes in the ground truth GRNs that do not have a match in this second data set, we only retained the genes that appear in both the Nanostring data set and the ground truth GRNs. For the ectoderm gene set, we had a match for 28 of 39 genes, so the number of relevant gene interactions in the ground truth ectoderm GRN was reduced from 85 to 57. With the small number of input genes, we found that PEAK gave similar prediction results. The result was 29 known gene-to-gene interactions were successfully predicted out of 57, yielding a sensitivity of 50.88% (Table 2).
For the endomesoderm data, 41 genes were present in both the ground truth GRN and the Nanostring data set, which reduced the 146 total known interactions to 76 known interactions. When TFs were speci ed, a total of 58 gene-to-gene interactions were successfully predicted out of 76, for a sensitivity of 76.32% (Table 2). Because more known connections in the endomesoderm ground truth GRN are supported by solid evidence at base-pair resolution and the Nanostring data set had far more experimental time points, it was expected that PEAK predictions using the Nanostring data set would yield the highest sensitivity for the endomesoderm GRN.
We further tested the impact of adjusting the parameter describing transcript turnover due to maternal and zygotic degradation mechanisms in terms of transcript half-life on the prediction results. Recent estimates of transcript turnover in sea urchin are in the range of 6 to 9 hours (Gildor et al. 2016). We explored the performance of PEAK for a range of median half-life times from 3 hours to 15 hours (Fig. 3).
In general, the sensitivity of the predicted results did not have any obvious increase or decrease trend with the increase of half-life. We obtained the highest sensitivity on the transcriptome data when the half-life was set to 7 hours for ectoderm genes. But when the half-life was set to 3 hours, we obtained the highest sensitivity on the Nanostring data set. Particularly when we used the set of 45 endomesoderm genes present in the Nanostring data, the sensitivity of the prediction was 76.32%.

Discussion
The computational prediction of a GRN has been explored in many organisms. In bacteria, for example, researchers used an integrative method to predict a GRN in B. subtilis with a large amount of input data (Arrieta-Ortiz et al. 2015). The B. subtilis study used more than 600 gene expression experiments and incorporated prior knowledge from the ground truth network to improve accuracy. The sensitivity of their GRN prediction is 74%, and they predicted 2,258 new regulatory interactions. The scale of experiments used in the Arrieta-Ortiz study is di cult to achieve in multicellular organisms. As is true with most other machine learning applications, the more prior knowledge available to train the algorithm, the better predictions one can expect the algorithm to produce. However, our approach using PEAK was able to yield 76% sensitivity using only gene expression data and no other prior knowledge.
Some organisms bene t from extremely large research communities producing extensive prior knowledge data sets in the form of ChIP assays, protein-protein interaction databases, functional gene annotation, extensive tissue expression information, TF binding sites, and known regulatory interactions. In these cases, the richness of prior knowledge has been harnessed to produce high quality novel network interaction predictions (Marbach et al. 2012b). However, the future of many disciplines will require knowledge from a wider variety of model species to uncover universal truths and mechanisms. The rst "omics'' step in many research programs establishing new model species is often the creation of a transcriptome, which inherently creates a quality source of gene expression data (Henry et al. 2010(Henry et al. , 2017Du et al. 2012;Helm et al. 2013;Tulin et al. 2013;Chen et al. 2014). New model species are especially important in evodevo as these emerging species are all poised to increase our understanding of how metazoan body plans evolved. The investigation of how changes in the gene regulatory program controlling embryonic development have shaped evolution is a particularly active research question and one that can bene t most directly from using GRN inference approaches to create model networks.
Therefore, our approach of only using gene expression data as input to generate regulatory interaction predictions is vital to the success of future research.
One goal of GRN prediction is to narrow the search space for edges to a subset of promising interactions to be further studied and validated experimentally. With 23,300 gene models in the sea urchin, there are ~ 540 million possible gene-to-gene interactions. By using time series gene expression data as input into the PEAK prediction algorithm for the 547 transcripts most likely to be a part of the regulatory program, we generated 14,802 predicted interactions that can be ordered by con dence or searched for speci c regulators or target genes. One caveat to our approach is the limitation that PEAK cannot predict selfinteractions where genes turn on or off their own expression. These self-interactions are known to be important to the precision and robustness of regulatory programs. Experimental design should take this limitation into account and make sure to check for self-regulation when investigating the targets of regulatory genes.
We also tested the ability of PEAK to make predictions using our data sets without specifying the set of TFs. Surprisingly, the sensitivity was not dramatically different than when the list of TFs was speci ed. Since the performance of PEAK was not overtly affected when the list of TF was not speci ed to the algorithm, those model systems where functional gene annotations are incomplete can still make use of our method. It is also possible to specify the transcript half-life to PEAK; however, we found that using the half-life value that most closely corresponds to previously published values for the organism did not always result in the best prediction result. Therefore, our recommendation is to use a range of half-life values to nd the best t to the data experimentally.

Conclusion
We found that using gene expression data as the sole input for machine learning GRN predictions was su cient to generate predictions that match known regulatory interactions when using the PEAK program. Our approach also generated new possible gene-to-gene interactions that are not currently described. The new predictions from our pioneer dataset will serve as a resource for the sea urchin community; a relatively small, but in uential group, in the areas of cis-regulatory biology and developmental regulatory networks. Our future research will include applying a similar approach to an emerging model species that transcriptomic gene expression data but does not yet have a developmental GRN model in order to generate predictions for new regulatory interactions. Our method is broadly applicable to any organism and is more accessible due to the ability to start with gene expression data that is often readily available. Although no prior knowledge is required, PEAK can accept many forms of prior knowledge to improve the quality of predictions (Altarawy et al. 2017). Especially in emerging evodevo model systems, where gene expression datasets are often the rst to be produced, PEAK shows promise to be able to assist network building efforts.

Data sets
Two gene expression data sets were obtained to use as input into the PEAK algorithm. The rst data set comes from the sea urchin transcriptome project (Tu et al. 2012) where 10 embryonic timepoints (labeled with time in units of hours-post-fertilization (hpf)) were assayed for transcript expression by RNAsEq. The sequenced transcripts in the transcriptome data set derived from cDNA collected from: (1) the unfertilized egg, 0 hpf; (2) cleavage stage, 10 hpf; (3) hatched blastula stage, 18 hpf; (4) mesenchyme blastula, 24 hpf; (5) the early gastrula, 30 hpf; (6) mid-gastrula stage, 40 hpf; (7) late-gastrula stage, 48 hpf; (8) prism stage, 56 hpf; (9) late prism stage, 64 hpf (10) the pluteus stage, 72 hpf. All embryonic samples were obtained from a single male and female mating pair, except the 24hr sample which was done separately as a pilot experiment. Only a single replicate was sequenced for each time point, presumably due to limitations of the amount of material that can be obtained from a single spawning event and the decision not to introduce biological variation due to individual differences if multiple urchins were used. Each sample generated 36.5M reads, of which 79% mapped to the S.purpuratus genome v3. Gene models were built by Cu inks, and, after quality ltering, 21,092 transcript models were de ned and assigned an 8digit WHL ID number beginning with "22." These models have been incorporated into the annotated sea urchin gene database and assigned to previously established "SPU ID" numbers. The existing GRN models for sea urchin development cover the time period from 0 to 30 hpf. Of the 10 time points sampled in the transcriptome data set from Tu, et al. only 5 are represented in the range of 0-30hpf. Therefore, we decided to also use data from a high density time series containing 49 time points, one time point per hour over the rst 48 hours of development and the unfertilized egg (Materna et al. 2010). The data set from Materna et al. uses Nanostring technology, which is often used as a gold standard for absolute quantitation due to the fact that it measures RNA directly as opposed to using an enzymatic reaction (Geiss et al. 2008). Even the Tu, et al. transcriptome data set was validated using independent Nanostring quantitation (Tu et al. 2014). A key difference between RNAseq and Nanostring is that Nanostring will only produce data for speci c known gene models for which probes were designed, whereas RNAseq surveys the entire transcriptome. The sea urchin Nanostring data set queried 161 regulatory gene products, including 130 transcription factors and 31 molecules involved in signaling. This gene set overlaps nicely with the gene set present in the ground truth sea urchin ectoderm and endomesoderm GRN models. Speci cally, 68 genes in total overlapped between the Nanostring probe set and the ground truth GRN genes. The overlap included 31 of 39 genes overlapping with the ectoderm GRN and 44 of 58 genes overlapping with the endomesoderm GRN. There were 7 genes in the Nanostring probe set that appear in both the ectoderm and endomesoderm GRNs (namely, not, otxa, foxA, eve, myc, bra, and hnf6). The normalized RNA counts produced by the Nanostring's Ncounter were used directly in our analysis.

Sea urchin GRN models
The current versions of the S. purpuratus GRNs for endomesoderm and ectoderm development are hosted by the Institute for Systems Biology and can be accessed online using the web application BioTapestry Interactive Network Viewer. The endomesoderm GRN can be found at http://grns.biotapestry.org/SpEndomes/, and the ectoderm GRN can be found at http://grns.biotapestry.org/SpEcto/. These two GRN models were built by a collaboration of sea urchin labs over the last thirty-some years. Each regulatory interaction is depicted as a directional line connecting two gene nodes, and each interaction is supported by experimental evidence, which can be accessed in the BioTapestry viewer directly. We obtained lists of the genes present in each network and a list of every gene-to-gene interaction present in the current version of the models from the BioTapestry director William Longabaugh. There are 39 genes represented in the ectoderm GRN and 58 genes represented in the endomesoderm GRN. The interaction list we obtained includes direct interactions and indirect interactions that are driven by signaling molecule intermediates. Interactions derived from signaling intermediates were not used in our comparison list. We also removed interactions where a gene regulates its own expression because the PEAK algorithm is not designed to be able to predict this type of interaction. The nal list used in our analysis contained 85 unique gene-to-gene interactions present in the ectoderm GRN (Additional le 5) and 146 unique gene-to-gene interactions present in the endoderm GRN (Additional le 6). These unique interactions made up our ground truth GRN models, which were used for comparison to the interactions predicted by the PEAK algorithm. The number of genes and connections included in our analysis when requiring a match between gene expression data set and corresponding ground truth network are described in Table 3. Preprocessing and differential gene expression determination The RNAseq data set was constructed with only one biological replica. Multiple methods have been developed to perform differential gene expression analysis on RNAseq data when only a single biological replica is available. These methods include: NOISeq (Tarazona et al. 2011), based on the multinomial distribution; GFold (Feng et al. 2012), based on the posterior distribution of log fold change; and EdgeR (Robinson et al. 2010), based on the negative binomial (NB) distribution. To discover quantitative changes in expression levels between experimental time points, we rst applied NOISeq, GFold, and EdgeR to determine the set of differentially expressed genes for further analysis.
For NOISeq, we normalized the data by RPKM (Reads Per Kilobase of transcript per Million mapped reads), which takes into account that more sequencing reads are generated from longer transcripts. The length of each transcript was obtained from the sea urchin database, Echinobase (https://www.echinobase.org). We omitted genes that had no record in the gene database. We set the simulation parameters as recommended in the NOISeq handbook, where the percentage (pnr) of the sequencing depth is pnr = 0.2, the number of samples to be simulated (nss) for each condition is nss = 5 and a small variability (v) is v = 0.02. We selected the differentially expressed genes with the higher NOISeq probabilities based on our chosen thresholds λ 1 = 0.9, λ 2 = 0.85. For GFold, we set the thresholds for the GFold value to λ 3 = ∓ 1, λ 4 = ∓ 1.5, since the GFold value is similar to the log2 fold change that is reliably used to select differentially expressed genes. For EdgeR, we set the log2 fold change (log2fc) cutoff as 2 and the edgeR dispersion as 0.01.
We used the gene database at Echinobase to map all WHL IDs to SPU_IDs to ensure that the IDs we analyzed before and after are consistent and unique.
Computational GRN prediction  Figure 1 Intersection of unique differentially expressed genes determined by NOISeq, EdgeR and GFold. The intersection of all 5 methods contains 10,627 unique genes speci ed as differentially expressed. NOISeq ( λ1= 0.9) and GFold2 ( λ4= ∓ 1.5) identify the fewest unique genes at 2 and 0, respectively. The intersection of the other 4 programs was increased by 2,755 when GFold2 (( λ4= ∓ 1.5)) was removed from the overlap. NOISeq ( λ1= 0.9) had the best balance of containing the most agreed upon gene set while adding the fewest unique genes.

Figure 2
A subset of PEAK predicted interactions for genes in the ground truth network with the top-5 predictions with the highest con dence scores for each gene. Both known interactions (red arrows) and new predicted interactions (black arrows) are included among these predictions.

Figure 3
Sensitivity as a measure of accuracy for the prediction of ectoderm and endomesoderm gene regulatory relations calculated with 5 different median mRNA half-life settings (3hrs, 5hrs, 7hrs, 10hrs, 15hrs) on the transcriptome data and the nanostring data.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.