Selecting the most appropriate time points to profile in high-throughput studies

Biological systems are increasingly being studied by high throughput profiling of molecular data over time. Determining the set of time points to sample in studies that profile several different types of molecular data is still challenging. Here we present the Time Point Selection (TPS) method that solves this combinatorial problem in a principled and practical way. TPS utilizes expression data from a small set of genes sampled at a high rate. As we show by applying TPS to study mouse lung development, the points selected by TPS can be used to reconstruct an accurate representation for the expression values of the non selected points. Further, even though the selection is only based on gene expression, these points are also appropriate for representing a much larger set of protein, miRNA and DNA methylation changes over time. TPS can thus serve as a key design strategy for high throughput time series experiments. Supporting Website: www.sb.cs.cmu.edu/TPS DOI: http://dx.doi.org/10.7554/eLife.18541.001


Introduction
Time series experiments are very commonly used to study a wide range of biological processes. Examples include various developmental processes (Roy et al., 2010), stem cell differentiation (Sperger et al., 2003), immune responses (Yosef and Regev, 2011), stress responses (Gitter et al., 2013) and several others. Indeed, analysis of the largest repository of gene expression experiments, the Gene Expression Omnibus (GEO), determined that roughly a third of these datasets come from experiments profiling dynamic processes over time (Zinman et al., 2013).
While mRNA gene expression data have been the primary source of high-throughput time series data, more recently several other genomic regulatory features are profiled over time. These include miRNA expression data (Schulz et al., 2013), ChIP-Seq studied to determine TF targets (Chang et al., 2013) and several types of epigenetic markers including DNA methylation (Singer et al., 2014), histone modifications (Paige et al., 2012) and more. In fact, with the rise in our ability to perform such high-throughput time series analysis, many researchers are now combining a few or several of these time series profiling experiments in a single experiment (Chang et al., 2013;Buenrostro et al., 2015) and then integrate these datasets to obtain a better understanding of cellular activity.
While integrated analysis of high-throughput genomic datasets can greatly improve our ability to model biological processes and systems, it comes at a cost. From the monetary point of view, these costs include the increased number of Seq experiments required to profile all types of genomic features. While such costs are common to all types of studies utilizing high-throughput data, they can be prohibitively high for time series based studies since they are multiplied by the number of time points required, the number of repeats performed for each time point and the number of different types of data being profiled. Importantly, even if the budget is not an issue, the ability to obtain enough samples for profiling all genomic features at all time points may be challenging, if not completely prohibitive.
One of the key determinants of the experimental and sample acquisition costs associated with time series studies is the number of time points that are being profiled. In most studies, the first and last time point can usually be determined by the researcher (for example, the time from birth to full lung structural development and maturation in mice). However, the number of samples required between these two points and the sampling frequency (given a fixed budget) are often hard to determine based on phenotypic observations since the molecular events of interest may precede such phenotypic events. To date, sampling rates have largely been determined using one of two adhoc protocols. The first utilized uniform sampling across the duration of the study  with the number of samples constrained by the available budget and samples. The second relied on some (conceived or real) knowledge of the process, often based on phenotypic observations. These studies, especially for responses though also for development, have often used nonuniform sampling (Schulz et al., 2013;Bar-Joseph et al., 2003a) though it is hard to determine if such sampling misses important molecular events between the sampled points.
Relatively, little work has focused so far on the selection of time points to sample in high throughput time series studies. Singh et al (Singh et al., 2005) and Rosa et al (Rosa et al., 2012) presented an iterative process which starts with profiling a small number of time points and then selects the next time point either based on an Active Learning method (Singh et al., 2005) or based on using prior related experiments (Rosa et al., 2012). Next the selected point is profiled and the process is repeated until a stopping criteria has been reached. Both of these methods require several iterations until the final time series is profiled, which can drastically lengthen the experiment time and can introduce additional biases making them less useful in practice. In addition, these methods employ a stopping criteria that does not take into account the full profile and also require that related time series expression experiments be used to select the point, which may be problematic when studying new processes or treatments.
Here, we propose the first non iterative method to address the issue of sampling rates across all different genomic data types. Our method starts by selecting a small set of genes that are known to be associated with the process being studied (while the full set is often unknown, for most processes a small set is usually known in advance). Next, we use a cheap array-based technology to sample these genes at a high, uniform rate across the duration of the study. Note that unlike standard curve fitting algorithms, a method for selecting time points for these experiments is required to accommodate over a hundred curves (for all genes) simultaneously, and we discuss various ways to formulate this as an optimization problem. To solve this optimization problems, we developed the Time Points Selection method (TPS), an algorithm that uses spline based analysis and combinatorial search to select a subset of the points that, when combined, provide enough information for reconstructing the values for all genes across all time points. The number of points selected can either be set in advance by the user (for example, based on budget constraints) or can be defined as a function of the reconstruction error. The selected time points are then used for the larger, genome-wide experiments across the different types of data being profiled.
To test and evaluate the method we applied it to study lung development in mice. Normal development of lung alveoli through the process of alveolar septation is a dynamic, coordinated process that requires the accurate spatial and temporal integration of signals. We currently lack a comprehensive understanding of the dynamic networks that govern normal alveolar septation. Thus, lung development can serve as an ideal test case for TPS since a variety of different time series genomic datasets are needed to enable accurate reconstruction of networks regulating this process. As we show, TPS was able to successfully identify time points for reconstructing the mRNA profiles of selected genes and these points improved upon uniform based sampling for such points. Further, we show that the set of points selected based on the analysis of this limited set of highly sampled mRNAs is also appropriate for sampling a much larger, unbiased, set of miRNA profiles as well as to determine the temporal protein levels of over 1000 proteins. Finally, we show that the mRNA samples can also be used to determine the optimal sampling points for a DNA methylation study of the same developmental process.

Results
The time points selection (TPS ) method We developed TPS to select a subset of k time points from an initial larger set of n points such that the selected subset provides an accurate, yet compact, representation of the temporal trajectory. Figure 1 presents an overview of the method. TPS utilizes splines to represent temporal profiles and implements a cross-validation strategy to evaluate potential sets of points. Following initialization which is based on the expression values, we employ a greedy search procedure that adds and removes points until a local minima is reached (Materials and methods). The resulting set is then used for the larger genomic and epigenetic experiments.
To test the usefulness of TPS , we used it to determine time points for a lung development study in mice. We first profiled the expression of 126 genes known or suspected to be involved in lung development using NanoString (See Appendix Methods for a list of the selected genes and the reason each was selected). We then used TPS analysis of these experiments to select a subset of time points for profiling the expression of a larger, unbiased, set of miRNAs. Finally, we have used TPS to design time series experiments to study DNA methylation patterns for a subset of the genes.

TPS identifies subset of important time points across multiple genes
We have tested the performance of TPS by using it to select subsets of points ranging from 3 to 25 and evaluating how well these can be used to determine the values of non-sampled points. To determine the accuracy of the reconstructed profiles using the selected points, we computed the average mean squared error for points that were not used by the method (Materials and methods). The results are presented in Figure 2. The figure includes a comparison of our method with two baseline methods: a random selection of the same number of points and uniform sampling of points within the range being studied, a method that is commonly used for time series expression profiling as discussed above. We have also compared the performance of the different strategies for initializing the set of points as discussed in Appendix Method (sorting by absolute differences or by equal partition) and between different methods for searching for the optimal subset (simulated annealing, weighting genes by cluster size, and adding/removing multiple time points per iteration, Appendix Methods). Finally, Figure 2 also presents the repeat noise values which is the theoretical limit for the performance of any profile reconstruction method.
As expected, we find significant performance improvement when using TPS when compared to randomly selected points. Importantly, we also see a significant and consistent improvement (for all sizes of selected time points) over uniform sampling highlighting the advantage of condition-specific sampling decisions. Sorting initial points by absolute values further improves the performance highlighting the importance of initialization when searching large combinatorial spaces. Simulated annealing, weighting, and multiple point selection improve performance as well. As the number of points used by TPS increases, it leads to results that are very close to the error represented by noise in the data (0.108) ( Figure 2-figure supplement 1). Figure 3 presents the reconstructed and measured expression values when using TPS to select 13 time points (less than a third of the points that were profiled). Note that even though each of these genes has distinct trajectory and inflection points, the selected set of time points enable TPS to fit all quite accurately without overfitting (See Identified time points using mRNA data are appropriate for miRNA profiling To test the usefulness of our method for predicting the correct sampling rates for other genomic datasets, we next profiled mouse miRNAs for the same developmental process. miRNAs have been known to regulate lung development (Sessa and Hata, 2013) and several miRNAs are differentially expressed during this developmental process (Williams et al., 2007). Several of these are also coordinately activated with various TFs to control specific transitions during development (Schulz et al., 2013). Thus, any large scale effort to model lung development would require the profiling of miR-NAs as well. Unlike the mRNA dataset, which utilized prior knowledge to profile less than 1% of all Figure 1. The TPS method. Clockwise from top left. Given a dense sampling of a selected subset of genes (a) we select an initial set of points (b) using the initialization method described in the text. Next, we fit a spline to the selected points for each gene (c) and evaluate the error on all other points. We perform a greedy search process (d) which iteratively removes and adds points to improve the test data fit resulting in the final set of points (e). The reconstructed curves are fitted to all genes (f) and an overall error is computed and compared to the theoretical limit (noise) to determine the ability of the selected number of points to fit the data. DOI: 10.7554/eLife.18541.002 The following figure supplements are available for figure 1: genes, the miRNA dataset contained a much larger number of miRNAs (600). Thus, the miRNA data represent an unbiased sample providing information on whether using one type of genomic data can be helpful for determining rates for other types. In our analysis, we normalized miRNA values by variance mean normalization (Bolstad et al., 2003).
To test TPS on this dataset, we used the mRNA expression data to select time points and then used the miRNA expression values for the selected time points to reconstruct the complete trajectories for each miRNA. The results are presented in Figure 4. As can be seen, when using the points selected based on the mRNA data we achieve a much lower error when compared to the error resulting from using the same number of uniform or random points (p<0:01 for random based on randomization analysis) highlighting the relationship between the two datasets and the ability to use one to determine points for the other. More generally, even though the noise in the miRNA data is      higher than for the mRNA dataset, relative ordering of the performance of each of the methods is similar to the mRNA results in Figure 2. This serves as a strong indication that mRNAs can serve as a general proxy for selecting time points for other genomic datasets. Figure 4b presents the error achieved when using the miRNA data itself to select the set of points (evaluated on the miRNA data). As expected, the performance when using the miRNA data itself is better than when using the mRNA data. However, when taking into account the inherent noise in the data the differences are not large. For example, when using the 13 selected mRNA points, the average mean squared error is 0.4312 whereas when using the optimal points based on the miRNA data itself the error is 0.4042. Figure 4-figure supplement 1 presents the reconstructed and measured expression values for a few miRNAs based on time points identified using the mRNA dataset. As with the mRNA data, the ability to accurately reconstruct different miRNA profiles highlights the importance of selecting a global set of points that can fit all genes and miRNAs in our study.
We have also analyzed the performance of TPS when using the mRNA data to select sampling time points for profiling the levels of more than 1000 proteins. We observed results that are very similar to the results obtained for the miRNA time point selection. Specifically, the points selected by TPS lead to reconstruction errors that are lower than those observed for uniform sampling or for a random set of the same number of points further demonstrating the general applicability of our method. See Appendix Results for details.

Using TPS to select time points for DNA methylation analysis
In addition to mRNA and miRNA expression data, epigenetic data have been increasingly studied in time series experiments (Talens et al., 2010;Schneider et al., 2010). To test the ability of the mRNA data to determine the appropriate points for DNA methylation analysis, we used targeted bisulfite sequencing to profile three CpG-enriched regions for 13 genes at 8 of the 42 time points used for the mRNA and miRNA studies (Materials and methods). We next applied TPS to the mRNA data of these 8 points to select the best subest of 4 points and compared the selected points to those that would have been selected using the methylation data itself. The 4 points identified using the mRNA data (0:5, 5, 15, 26) were exactly the same as the ones selected using the methylation data indicating again that mRNA data is a good proxy for the dynamics of the epigenetic data as well. Figure 5-figure supplement 1 presents the reconstructed splines over the identified points for several genomic methylation loci. Figure 5 presents the methylation and expression curves for 3 genes: Akt,1 Cdh11, and Tnc. These were the genes with the strongest negative correlation between their methylation and expression. As can be seen, in several cases we observed strong negative or positive correlations between the two datasets in the time points we used serving as another indication for the ability to use one dataset to select the sampling points for the other. See

Discussion
Time series gene expression experiments are widely used in several studies. More recently, advances in sequencing and proteomics are enabling the profiling of several other types of genomic data over time. Here we focused on lung development in mice with the goal of identifying an optimal set of time points for profiling various genomic and proteomic data types for this process.
An important question is: Whether a better selection of time points really leads to observations that are missed when using an inferior set of points (even if the number of points is the same)? To answer this question we looked at several prior studies that profiled mouse lung development over time using various high throughput assays. Table 1 presents 9 representative studies and lists the biological data that was profiled and the time points that were used. As can be seen, while certain time points seem to be widely used across studies (for example, 7d) others were profiled in only one or two of the studies (2d, 10d, three weeks). This raises several issues. First, it is very hard to compare or combine these datasets (for example, protein levels were not profiled on day 7 (Cox et al., 2007) whereas all mRNA levels were). It also makes it hard to determine if differences between DE genes or miRNAs between these studies are the result of differences in the underlying conditions studied (for example, when testing for mutants or treatments) or simply the result of different sampling. Finally, each of these studies may have missed key genes, proteins or miRNAs because of the sampling used restricting the ability of downstream analysis to use the data to model causal and regulatory events in lung development.
To illustrate these problems we compared the resulting curves using three of the sampling rates from Table 1 to the reconstructed curves obtained by using TPS to select the optimal 5 and 8 time points. For example, the points selected by Schulz et al. (2013) are 0, 4, 7, 14 and 28 (since 28 is last day in our analysis we used it instead of 42). In contrast, TPS selects 0:5, 6, 9:5, 19 and 28. As can be seen in Figure 6, important expression changes in key genes are missed by using the arbitrary points while the TPS points are able to correctly reconstruct these profiles even though the total number of points is the same (5). More globally, the error for the arbitrary set of selected points is much higher on average (Appendix 2- Table 4). Similar results are obtained for the other sampling rates used in the past ( Figure 6, Appendix 2- Table 4) and when comparing TPS to iterative methods previously suggested for selecting the set of points to profile (Figure 1-figure supplement 1). This indicates that accurate selection of time points can have a large impact on the ability of the study to identify   6. Comparison of TPS with sampling rates used in previous studies. Dark green curves are the reconstructed profiles based on the points profiled by prior studies. Light green and red curves are based on the points selected by TPS . As can be seen, even when comparing results from using the same number of points, TPS can identify key events for some of the genes that are missed when using the phenotype based sampling rates. Figure 6 continued on next page key genes and events. See also Appendix Results for a discussion about the importance of the differences between the TPS and prior work results for selected genes.
Our method relies on a very small subset of genes that are known to be involved in the process studied for the initial (highly sampled) set of experiments. While such set is known for several processes, there may be cases where very little is known about the biological process and so it may be hard to obtain such set. TPS can still be applied to determine sampling rates for such processes using a small random set of genes. To illustrate this we repeated the analysis presented in Results using only the measured values of 25% of genes in our original set and replacing the values for the other genes with random profiles. As we show in Figure 2-figure supplement 2, even when using such set, the time points selected by TPS greatly improve upon an arbitrary set of the same number of time points. Since in most time series experiments at least 25% of the genes are differentially expressed (and in several cases a much larger fraction, (Zhou et al., 2009;Shi et al., 2015) a random selection of genes is likely to exhibit similar results even for poorly understood processes.
Beyond the analysis of a specific type of data, several studies have now been profiling multiple types of genomic data over time. Such studies need to agree on a set of time points which would be common to all experiments so that these diverse types can be integrated to form a unified model (Chang et al., 2013;Roy et al., 2010). To date, the selection of such points relied on ad-hoc methods. The processes being studied were either sampled uniformly or based on prior knowledge. However, known properties of such systems were often been based on phenotypic observations which may not necessarily agree with the timing of molecular events. In addition, in many case studies of the same, or similar processes differed with respect to the time points that have been profiled. For example, early work on the analysis of cell cycle data in yeast utilized both uniform and nonuniform sampling (Spellman et al., 1998) and recent studies of circadian rhythms have followed a similar pattern (Storch et al., 2002;Ueda et al., 2002). Similarly, more recent analysis of responses to flu diverged widely in the (nonuniform) sampling rates that were used (Shapira et al., 2009;Li et al., 2011).
TPS addresses these problems by using a principled method for determining sampling rates. An important goal in the development of TPS was to enable it to be successfully applied to different types of biological datasets. As we show, a relatively inexpensive, gene centric, method provides a very good solution for RNA expression profiling as well as other types of data including miRNAs and DNA methylation. Thus, a combined experiment can be fully designed using our method.
While we evaluated TPS on several types of high throughput data, we have only tested it so far on data for a specific biological process (lung development in mice). While we believe that such data is both challenging and representative and thus provides a good test case for the method, analysis of additional datasets may identify new challenges that we have not addressed and we leave it to future work to address these.
TPS, including all initialization methods discussed, is implemented in Python and is available on the supporting website. We hope that as sequencing technology continues to advance, more and more studies would integrate diverse types of time series data and will utilize TPS in the design pipeline of their studies.

Materials and methods mRNA and miRNA used in the study
To select the list of 126 genes used in the NanoString profiling we searched the literature for genes that have been linked to the following processes: (a) Cell type specification genes (e.g. alveolar type I epithelial, alveolar type II epithelial, any epithelial, basal, endothelial, mesenchymal, pericyte, fibroblast, monocyte), (b) genes known to be up or down regulated during septation, (c) genes known to be altered in DNA methylation during development, (d) genes known to be involved in septation, (e) genes known to be regulated by miRNA involved in septation, and (f) genes known to be regulated by DNA methylation during fibrosis. Appendix 2- Table 1 contains a list of the selected genes and the process for which they were selected.
For the miRNA set we used a commercially available, unbiased, array (the nCounter Mouse miRNA Expression Assay Kit, NanoString).

mRNA and miRNA profiling and analysis
A total of 240 samples were isolated by Laser Capture Microscopy (LCM) from murine lung at multiple time points (E16.5, P.05 to P14 every 12 hr, and P15 to P28 every 24 hr). The samples were used to prepare total RNA. RNA extraction was performed by miRNeasy MicroKit (Qiagen) following the manufacturer's protocol. RNA concentration and integrity were measured by using NanoDrop ND-2000 and 2200 Tape Station. A custom NanoString probe set (Reporter Code set and Capture Probe set) for 126 genes was designed and the nCounter Gene Expression Assay was performed using 50 ng total RNA. The data files produced by the nCounter Digital Analyzer were exported as a Reporter Code Count (RCC) file and data normalization was performed using the nSolver, the analysis software provided by Nanostring.

DNA methylation analysis
Mouse alveolar lung tissues attached to LCM caps were stored at À80˚C until processing. DNA was extracted using the ZR Genomic DNA-Tissue MicroPrep kit (Zymo Research). Incubation with Digestion buffer and proteinase K was done overnight at 55˚C in inverted tubes. 13 genes were chosen for targeted NextGen bisulfite sequencing (NGBS): Igfbp3, Wif1,Cdh11,Eln,Sox9,Tnc,Dnmt3a,Akt,Vegfa,Lox,Foxf2,Zfp536 and Src, based on published data (Cuna et al., 2015). Targeted NGBS was done on samples collected at: E16.5, E18.5, P0.5, P1.5, P2.5, P5, P10, P15, P19 and P26. Multiplex PCR was performed using 0.5 units of TaKaRa EpiTaq HS (Takara Bio, Kusatsu, Japan) in 2x master mix. FASTQ files were aligned using open source Bismark Bisulfite Read Mapper using Bowtie2. Methylation levels were calculated in Bismark. Sites where the difference in methylation was less than 5% over the entire time period, those where there was a difference of >20% at a single time point and those with less than 3 non zero values were removed from the analyses.

Problem statement
Our goal is to identify a (small) subset of time points that can be used to accurately reconstruct the expression trajectory for all genes or other molecules being profiled. We assume that we can efficiently and cheaply obtain a dense sample for the expression of a very small subset of representative genes (here we use nanostring to profile less than 0.5% of all genes) and attempt to use this subset to determine optimal sampling points for the entire set of genes.
Formally, let G be the set of genes we have profiled in our dense sample, T ¼ ft 1 ; t 2 ; . . . ; t T g be the set of all sampled time points. We assume that for each time point we have R repeats for all genes. We denote by e r gt be the expression value for gene g 2 G at time t 2 T in the r'th repeat for that time point. We define D g ¼ fe r gt ; t 2 T; r 2 R as the complete data for gene g over all replicates and time points T.
To constrain the set of points we select, we assume that we have a predefined budget k for the maximum number of time points we can sample in the complete experiment (i.e. for profiling all genes, miRNAs, epigenetic marks etc. using high-throughput seq experiments). We are interested in selecting k time points from T which, when using only the data collected at these k points, minimizes the prediction error for the expression values of the unused points. To evaluate such a selection, we use the selected values to obtain a smoothing spline (De Boor, 1978;Bar-Joseph et al., 2003a;Wahba, 1990) function for each gene and compare the predicted values based on the spline to the measured value for the non-selected points to determine the error. In our problem, t 1 and t T define the first and end points, so they are always selected. The rest of the points are selected to maximize the following objective 1: Problem statement: Given D g for genes g 2 G, the number of desired time points k, identify a subset of k À 2 time points in T n ft 1 ; t T g which minimizes the prediction error for the expression values of all genes in the remaining time points.

Spline assignments
Before discussing the actual procedure we use to select the set of time points, we discuss the method we use to assign splines based on a selected subset of points for each gene. There are two issues that need to be resolved when assigning such smoothing splines: (1) The number of knots (control points) and (2) their spacing. Past approaches for using splines to model time series gene expression data have usually used the same number of control points for all genes regardless of their trajectories (Subhani et al., 2010;Bar-Joseph et al., 2003b), and mostly employed uniform knot placements. However, since our method needs to be able to adapt to any size of k as defined above, we also attempt to select the number of knots and their spacing. We do this by using a regularization parameter for the fitted cubic smoothing spline where number of knots is increased until the smoothing condition is satisfied (Wahba, 1990). The regularization parameter is estimated by leaveone-out cross-validation (LOOCV).

TPS : Iterative process to select points
Because of the highly combinatorial nature of the time points, we rely on a greedy iterative process to select the optimal points as summarized in Figure 1 (See Appendix Methods for pseudocode).
There are three key steps in this algorithm which we discuss in detail below.
. Selecting the initial set of points: When using an iterative algorithm to solve non-convex problems with several local minima, a key issue is the appropriate selection of the initial solution set (Hartigan, 1975;McLachlan and Peel, 2004)]. We have tested a number of methods for performing such initializations and results for some of these are presented in Figure 1-figure supplement 2. Since the goal of the method is to optimize a specific function (error on the left out set of expression values measured at time points not used), all initialization methods can be tested for each dataset and the solution minimizing the left out error can be used. See Appendix Methods for details.
. Iterative improvement step: After selecting the initial set, we begin the iterative process of refining the subset of selected points. In this step we repeat the following analysis in each iteration. We exhaustively remove all points from the existing solution (one at a time) and replace it with all points that were not in the selected set (again, one at a time). For each pair of such point, we compute the error resulting from the change (using the splines computed based on the current set of points evaluated on the left out time points), and determine if the new point reduces the error or not. Formally, let T À ¼ T n ft 1 ; t T g and C n be set of points for iteration n. We are interested in finding a point pair ðt a 2 C n ; t b 2 T À n C n Þ which minimizes the following error ratio for the next iteration C nþ ¼ C n n ft a g [ ft b g: whereê Cn gt is our spline based estimate of the expression of gene g at time t by fitting smoothing spline over points C n . If there are pairs which lead to an error ratio of less than 1 in the above function, we select the best (lowest error), assign it to C nþ1 and continue the iterative process. Otherwise we terminate the process and output C n as the optimal solution. While the process is guaranteed to converge, given the large combinatorial search space convergence can be slow. This makes adequate initialization an important issue which we have focused on. In practice we find that the search usually converges very fast (within 10 -15 iterations).
. Fitting smoothing spline: The third key step of our approach is fitting a smoothing spline to every gene independently for the selected subset of time points. As discussed above, this is done by using a regularized version of approximating splines which allow us to determine a unique number of control points and spacing for each of the genes. See Appendix Methods for more details.

Individual vs. cluster-based evaluation
So far, we assumed that error of each gene has the same contribution to the overall error. However, this assumption ignores the fact that the expression profiles of genes are correlated with the expression of other genes. To take the correlation between gene profiles into account, we also performed cluster based evaluation of genes where we analyzed the error by weighting each gene in terms of inverse of the numbers of genes in the cluster it belongs. This scheme ensures that each cluster contributes equally to the resulting error rather than each gene. We find clusters by k-means algorithm over time series-data by treating each gene as a point in R T space as well as over a vector of randomly sampled T time points on fitted spline (Bishop, 2006). We use Bayesian Information Criterion (BIC) to determine the optimal number of clusters (Schwarz, 1978).

Acknowledgement
We thank the LungMAP consortium for useful comments regarding the methods and analysis presented in this paper. Work supported in part by NIH grant U01 HL122626.

Funder Grant reference number Author
National Institutes of Health U01 HL122626 Ziv Bar-Joseph The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Ethics Animal experimentation: This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols (APN 10042) of the University of Alabama at Birmingham. All lungs were isolated immediately following euthanasia using approved protocols.

Major datasets
The following dataset was generated: Database, license, and accessibility Table 1 provides the list of genes used for the nanostring analysis and the rational for their inclusion.

DNA Methylation analysis
Mouse alveolar lung tissues attached to LCM caps were stored at À80˚C until processing. DNA was extracted using the ZR Genomic DNA-Tissue MicroPrep kit (Zymo Research). Incubation with Digestion buffer and proteinase K was done overnight at 55˚C in inverted tubes. 13 genes were chosen for targeted NextGen bisulfite sequencing (NGBS): Igfbp3, Wif1,Cdh11,Eln,Sox9,Tnc,Dnmt3a,Akt,VEGF,Lox,FoxF2,ZFP536 and Src, based on published data (Cuna et al., 2015). The presence of CpG islands in 5-UTR, gene body and 3-UTR was interrogated using NCBI Epigenomics database, as well as CpG island searcher (Takai and Jones, 2002), and EMBOSS Cpgplot (Rice et al., 2000). Targeted NGBS was done by Epigendx Inc. Gene sequences from selected regions were acquired from the Ensembl database. Gene IDs, transcript IDs, simplex PCR IDs, and target regions for each gene are listed in Appendix 2-

TPS Algorithm
A pseudocode for the TPS algorithm is presented in Algorithm 1.3.
Algorithm 1. TPS : Iterative k-point selection 1: Procedure ITERATIVE-TEMPORAL-SELECTION 2: C 0 ¼ select initial k time points by absolute difference sorting 3: e 0 ¼ error of remaining points by fitting splines to C 0 4: for each pair ðt a ; t b Þ 2 ðT À n C i Þ Â C i do 7: e Ã ¼ estimate error by fitting smoothing spline to C Ã where regularization parameter is estimated by LOOCV 9: if e Ã <e i then 10: end if 13: i ¼ i þ 1 14: end for 15: While e iþ1 <e i 16: Output C i and e i 17: end procedure Selecting the initial set of points When using an iterative algorithm to solve non-convex problems with several local minima, a key issue is the appropriate selection of the initial solution set. We have tested a number of methods for performing such initializations. The simplest method we tried is to uniformly select a subset of the points (so if k ¼ T=4 we use each 4'th point). Another method we tested is to partition the set of all time points T into k À 1 intervals of almost equal size. This method determines these boundaries by estimating the cumulative number of points until each time point and selecting time points with cumulative values T kÀ1 ; 2 T kÀ1 ; . . . ; ðk À 2Þ T kÀ1 respectively. Then, it uses k interval boundaries including t 1 and t T as initial solution. We also tested a method that relies on the changes between consecutive time points to select the most important ones for our initial set. Specifically, we sort all points except t 1 and t T by average absolute difference with respect to its predecessor and successor time points by computing: where Mdðe gti Þ is the median expression for gene g at time t i . We then select the k À 2 points with maximum m ti as the initial solution.  We found that for our particular dataset, the dynamic initialization with MetricA e;ti performed best for selections of time points smaller than one third of the the initial dense time series, while the non dynamic m ti method works best for selections of time points between one third and and one half of the initial time series. The dyanmic metric and non dynamic metrics can be compared in their performance on our data in Figure 1-figure supplement 2. However, all of the metrics performed much better than a selection of random points as shown in Figure 1-figure supplement 3.

Further improvements to the iterative points selection procedures
We tested the following possible search strategies to improve the iterative points removal and addition in TPS.
. We add and remove b time points in each iteration instead of a single point. This increases the complexity of each iteration from OðkGT 2 QÞ to OðkGT 2b QÞ where Q is the complexity of fitting a smoothing spline.
. We use simulated annealing to escape from local minima (Kirkpatrick et al., 1983). In this case, we do not always move to a pair of points with the minimum error in each iteration, but instead move to a solution with random pair of points with probability 1 if its error e r is lower than error of current solution e i whereas we move to a solution with probability e ÀCðe r Àe i Þ if e r ! e i . Here, C is the temperature that increases by increasing number of iterations and the probability of moving to a solution with larger error decreases over time.
In practice, even though both approaches should in theory be better able to escape local minima than the greedy approach described above, for the data we analyzed they do not perform significantly better as Figure 2 in the main text demonstrates.
For our specific setting, we also introduce a regularization parameter to enable us to determine the number of control points. Let I g ¼ fðt; Mdðe gt ÞÞ; t 2 Cg, and be the spline we are interested in fitting, smoothing spline can be found by the following optimization problem which minimizes penalized least-squares error: where l is the regularization parameter which prevents overfitting by affecting the number of knots selected. We estimated l by leave-one-out cross-validation (LOOCV) in our experiments (See Appendix Methods for details of smoothing spline fitting).

Proteomics analysis
Proteins were extracted using tissue protein extraction reagent (T-PER, Thermo) as per manufacturer's instructions, carried out directly on the micro-dissection cap. Protein concentrations was determined with the EZQ protein assay (Life Sciences). The proteins were digested overnight at 37C, followed by acidification to pH 3 À 4 with 10%formic acid (FA), and extracted as per manufacturer's instructions, then concentrated to near completion using a Savant SpeedVac Concentrator (Thermo) and diluted with 0:1% FA to a final concentration of~100 ng \uL for analysis by LCMS. The LCMS data were converted to a universal MzXML file format prior to being searched using SEQEST (Thermo) against a Mouse subset of the UniRef100 database. These data were then uploaded to Scaffold (Proteome Software) in order to filter and group each peptide ID to specific proteins with peptide probability scores set at 80%, and protein probability scores set at 99%. Using only proteins presenting with 2 or more peptides per protein, the confidence interval was set to ~99:9% with and FDR <0:1. Quantification was carried out using Scaffold Q + using normalized spectral counts.

Example of a TPS run
Here we discuss a specific setting for TPS that allows us to discuss the set of points selected and their relevance. Specifically, to test TPS , we fixed three set points in advance (first (0:5'th day) and last (28'th day), which are required for any setting and day 7 which was previously determined to be of importance to lung development. Next, we have asked TPS to further select 10 more points (for a total of 13). For this setting, the method selected the following points: 0:5, 1:0, 1:5, 2:5, 4, 5, 7, 10, 13:5, 15, 19, 23, 28. While we do not know the ground truth, the larger focus on the earlier time points determined by the method (with 7 of the 13 points for the first 7 days) makes sense in this context as several aspects of lung differentiation are determined in the first week (Guilliams et al., 2013). The other 3 weeks were more or less uniformly sampled by our TPS . This highlights the usefulness of an unbiased approach to sampling time points rather than just uniformly sampling through the time window.

TPS identifies subset of important time points across multiple genes
To understand whether gene-expression profiles over time has a simple trend, we also compare the reconstruction performance of TPS with fitting piecewise linear curves between initial and middle time points and between middle and last time points. The reconstruction error by TPS is significantly better than the piecewise linear reconstruction for 102 genes out of 126 genes.
We have plotted the comparison of reconstruction for several of these genes as in Figure 2figure supplement 3. The distribution of error difference between these methods looks significantly different than normal distribution (p<0:0001 by Shapiro-Wilk test).

miRNA clusters are enriched for several biological processes
While the mRNA datasets includes only a handful of genes (less than 0.5% of all genes) the miRNA data includes more profiles and so further analysis of this data can be perfromed. We have performed clustering of the miRNA data using k-means (Hartigan, 1975) where the number of clusters is selected by Bayesian Information Criteria (Schwarz, 1978) leading to 8 stable miRNA clusters Figure 4-figure supplement 2. Next, we mapped miRNA's to predicted targets using TargetScan (Agarwal et al., 2015), and performed gene-enrichment analysis by FuncAssociate (Berriz et al., 2003). We find clusters to be enriched for several Gene Ontology biological processes (Ashburner et al., 2000). For instance, cluster 4 is enriched for singleorganism cellular process, positive regulation of biological process, regulation of metabolic process, etc. See Supporting Website for complete results.  , mmu-miR-136 targets Tgfb2, mmu-miR-152 targets Meox2, Robo1, Fbn1, Nfya (Popova et al., 2014). Additional figures for all miRNAs and mRNAs are avialable on the supporting website.

TPS application to select time points for proteomics analysis
We used mass spectrometry to profile the levels of 1020 proteins over the optimal 13 time points determined by TPS (using the mRNA expression data): ½0:5; 1:0; 1:5; 2:5; 4:0; 5:0; 7:0; 10:0; 13:5; 15:0; 19:0; 23:0; 28:0). To test the ability of TPS to determine the optimal time points for the proeomics data (based only on the mRNA data) we performed a similar analysis to the analysis performed for the miRNA data. Specifically, we used TPS to select subset of 4 to 12 of these points based on the mRNA data and compared the error using these points to random and uniform selection of the same number of points. The results are presented in . In addition to comparing TPS to random and uniform we have also compared different strategies for initializing the set of points as discussed in Method. Finally, the figure also presents the repeat noise values which is the theoretical limit for the performance of any profile reconstruction method.
As for the miRNA data, we see a significant and consistent improvement (for all number of selected time points) over uniform sampling highlighting the advantage of condition specific sampling decisions. Again, as the number of points used by TPS increases, it leads to results that are very close to the error represented by noise in the data (17:47).

Analysis of methylation data
Methylation data included 3 repeats for time points 0:5, 1:5, 2:5, 5, 10, 15, 19, 26 for 266 loci belonging to 13 genes. Among these genes all except Zfp536 were also profiled in our nanostring mRNA analysis. Appendix 2- Table 2 summarizes the number of loci for each gene in the methylation dataset. We used shifted percentage of methylation at each time point in our analysis which is obtained by subtracting the median percentage of methylation at initial time point (baseline) from all data points for each gene. Figure 5-figure supplement 2 presents the best positive or negative correaltion observed between the methylation data and the gene expression data for these genes (note that we do not expect all up stream regions to show a correlated profile since it is likely that only a subset, or even a single, region is responsible for the changes in expression observed which is why we look for the most correlated or anti-correlated region).

Importance of correct determination of expression profiles
As shown in Figure 6 in the main text, TPS results differ from prior methods when reconstructing expression profiles for several genes. Below we discuss the significance of these differences and their impact on the ability to correctly assign function to that gene: . Nol3: Nucleolar protein 3 (apotosis repressor with CAR domain) gene (also called ARC) encodes a protein that inhibits apoptosis, by decreasing activities of Caspases 2 and 8 and tumor protein p53. Evaluation of the TPS profile suggests that the increase in Nol3 correlates with postnatal lung development, with a rapid increase from birth until 2 weeks of age, followed by stabilization, while the prior sampling rates show only an initial peak and then decrease. While the exact role of Nol3 in lung development has not been established, it is known that Nol3 protects pulmonary arterial smooth muscle cells from hypoxia-induced death and facilitates growth factor-induced proliferation and hypertrophy, and is probably involved in human pulmonary hypertension (Turi et al., 1990). Nol3 is a regulator of myogenic differentiation (Hunter et al., 2007) and its pattern of expression suggests that it may be important in regulating pulmonary airway and vascular smooth muscle development and differentiation.
. Esr2: The gene estrogen receptor beta encodes a receptor for estrogen, and is important in regulating lung development and modulating differences in lung development between males and females (Gortner et al., 2013). Evaluation of the TPS profile suggests that the Esr2 decreases briefly after birth, followed by an increase from around day 5 until day 20 whereas non-optimized profile suggests a relatively flat profile. While fetal mouse lungs express both Esr2 alpha and beta, adult mouse lungs express only Esr2 beta consistent with the TPSresults (Carvalho and Goncalves, 2012).
. Igfbp3: Insulin-like growth factor binding protein 3 ( Igfbp3) belongs to the Igfbp family and has a Igfbp domain and a thyroglobulin type-I domain (http://www.ncbi.nlm.nih.gov/gene/ 3486). The TPS profile for Igfbp3 is very different from the non-optimized profile, suggesting that important biological information is lost when not using the TPS profile. Igfbp3 regulates the induction of TNC by TGF-beta (Brissett et al., 2012) and both these molecules are critical in lung alveolar septation.
. Wif1: Wnt inhibitory factor 1 ( WIF1) inhibits Wnt proteins, that are well known to be critical in many stages of lung development. The TPS profile is very different from the non-optimized profile, as the TPS profile indicates a much earlier and higher peak of WIF1 during postnatal lung development that may be critical in alveolar septation. WIF1 is a target gene for Smad1, one of the BMP receptor proteins important in lung development and maturation. A regulatory loop of Bmp4-Smad1-Wif1-Wnt/beta-catenin may coordinate BMP and Wnt pathways to control lung development (Xu et al., 2011), and dysregulation of the Smad1/Wif1 axis is associated with lung hypoplasia (Fujiwara et al., 2012).
. Inmt: Indolethylamine N-methyl transferase (Inmt) gene encodes an enzyme that N-methylates indoles such as tryptamine (http://www.ncbi.nlm.nih.gov/gene/11185). The TPS profile for Inmt is very different from the non-optimized profile, as the TPS profile indicates a much lower and prolonged reduction of Inmt during postnatal lung development. Methyl conjugation is an important pathway in the metabolism of many drugs, neurotransmitters, and xenobiotic compounds (Thompson and Weinshilboum, 1998). While it is known that Inmt expression varies over the course of human lung development (Kopantzev et al., 2008), its exact role in lung development is not known.
. Fgf18: Fibroblast growth factor 18 (Fgf18) is a member of the fibroblast growth factor family, and the Fgfs are well known to be critical in multiple stages of lung development. The nonoptimized profile indicates a smaller and later peak, and is not similar to the TPS profile which suggests a much more improtant role.Fgf18 is a pleiotropic growth factor that stimulates proliferation in a number of tissues (http://www.ncbi.nlm.nih.gov/gene/8817). Fgf18 is highly expressed in the developing lung as the TPS profile indicates (Ohbayashi et al., 1998), and Fgfr3 is important in postnatal alveolar development (Weinstein et al., 1998).
The role of Fgf18 in regulating fibroblast proliferation (Hu et al., 1999) may be important in alveolar septation, as Fgf18 increases after birth with a peak around P10, with reduction after completion of alveolar septation.