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

  1. Michael Kleyman
  2. Emre Sefer
  3. Teodora Nicola
  4. Celia Espinoza
  5. Divya Chhabra
  6. James S Hagood
  7. Naftali Kaminski
  8. Namasivayam Ambalavanan
  9. Ziv Bar-Joseph  Is a corresponding author
  1. School of Computer Science, Carnegie Mellon University, United States
  2. University of Alabama at Birmingham, United States
  3. University of California, United States
  4. CARady Children’s Hospital San Diego, United States
  5. School of Medicine, Yale University, United States


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:


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 ad-hoc protocols. The first utilized uniform sampling across the duration of the study (Li et al., 2013) 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.


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.

Figure 1 with 3 supplements see all
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.

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.

Figure 2 with 4 supplements see all
Performance of TPS using different sizes for the selected points.

Error comparisons of TPS variants to uniform selection of points and noise. Absolute difference - Greedy iterative addition with absolute difference initialization (Algorithm 1, Appendix Methods). Simulated annealing - Iterating using simulated annealing with absolute difference initialization. Weighted error - Selection based on cluster rather than individual gene errors. See Appendix Methods for details.

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 Figure 3—figure supplement 1 and Figure 3—figure supplement 2 for figures of several other genes and for figures reconstructed by using the best 8 time points as determined by TPS , respectively).

Figure 3 with 2 supplements see all
Reconstructed expression profiles for selected genes.

(a). Pdgfra. , (b). Eln. , (c). Inmt.

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 miRNAs as well. Unlike the mRNA dataset, which utilized prior knowledge to profile less than 1% of all genes, the miRNA dataset contained a much larger number of miRNAs (6^00). 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 with 3 supplements see all
Performance of TPS by on the miRNA data.

(a) TPS reconstruction error when using the mRNA data to select time points for the miRNA experiments. Results of random and uniform selection as well as repeat noise error are also presented for comparison. TPS variants shown are the same two presented in Figure 2. (b) Error of splines with points selected by training TPS on the actual miRNA data itself, using the maximum absolute difference initialization.

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 Figure 5—figure supplement 2 for correlation of all genes.

Figure 5 with 2 supplements see all
Comparison of gene expression and methylation data for selected genes.

(a). Akt1. , (b). Cdh11. , (c). Tnc.


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.

Table 1

Summary of prior high throughput lung development studies.
ReferenceData typesSelected time points (Days)

[Bonner et al., 2003]

mRNA expression

E9, E4, E17, 0, 7, 14, 28

[Melén et al., 2011]

mRNA expression

E16, E18, 0, 7, 14, 28

[Bhaskaran et al., 2009]

microRNA expression

E16, E19, E21, 0, 6, 14, 60

[Dong et al., 2011]

 mRNA and microRNA expression

E12, E14, E16, 0, 2, 10

[Cox et al., 2007]

Protein expression levels

E12, E14, E18, 2, 14, 56

[Schulz et al., 2013]

 mRNA and miRNA expression

0, 4, 7, 14, 42

[Cormack et al., 2010]

 mRNA expression

0, 7, 14, adult

[Mager et al., 2007]

 mRNA expression

 E15, E17, E19, E21, 1, 14, 84

[Mariani et al., 2002]

 mRNA expression

E18, 1, 4, 7, 10, 14, 21, adult

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 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.

Figure 6 with 1 supplement see all
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. Subfigures a,b, and c are a piecewise linear fit over points 0.5, 7.0, 14.0, 28.0 . Subfigures d,e, and f are a piecewise linear fit over points 0.5, 2.0, 14.0, 28.0. Subfigures g,h, and i are a piecewise linear fit over points 0.5, 4.0, 7.0, 14.0, 28.0.

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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={t1,t2,,tT} 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 egtr be the expression value for gene gG at time tT in the r’th repeat for that time point. We define Dg={egtr,tT,rR 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, t1 and tT 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 Dg for genes gG, the number of desired time points k, identify a subset of k-2 time points in T{t1,tT} which minimizes the prediction error for the expression values of all genes in the remaining time points.

Spline assignments

Request a detailed protocol

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 leave-one-out cross-validation (LOOCV).

TPS : Iterative process to select points

Request a detailed protocol

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{t1,tT} and Cn be set of points for iteration n. We are interested in finding a point pair (taCn,tbT-Cn) which minimizes the following error ratio for the next iteration Cn+=Cn{ta}{tb}:

    (1) error ratio=error(Cn+)error(Cn)=gGrRtTCn+(e^gtCn+egtr)2gGrRtTCn(e^gtCnegtr)2

    where e^gtCn is our spline based estimate of the expression of gene g at time t by fitting smoothing spline over points Cn. 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 Cn+1 and continue the iterative process. Otherwise we terminate the process and output Cn 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

Request a detailed protocol

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 RT 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).

Appendix 1 Methods

Selecting the set of 126 genes

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—Table 3. A total of 42 target PCRs were designed by PyroMark Assay Design Software (Qiagen).

Targeted NGBS was done on samples collected at the following time points: E16.5, E18.5, P0.5, P1.5, P2.5, P5, P10, P15, P19 and P26. Mouse genomic DNA (200–500 ng) was bisulfite treated using the EZ DNA Methylation Kit (Zymo Research). Multiplex PCR was performed using 0.5 units of TaKaRa EpiTaq HS (Takara Bio) in 2x master mix.

FASTQ files were aligned using open source Bismark Bisulfite Read Mapper using Bowtie2. Methylation levels were calculated in Bismark by dividing the number of methylated reads by the number of total reads, considering all CpG sites covered by a minimum of 30 total reads. 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.

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:   C0= select initial k time points by absolute difference sorting

3:   e0= error of remaining points by fitting splines to C0i=0

4:   i=0

5:   do

6:    for each pair (ta,tb)(T-Ci)×CiC*=Ci{ta}{tb}e*= do

7:      C=Ci{ta}{tb}

8:      e= estimate error by fitting smoothing spline to C* where regularization parameter is          estimated by LOOCV

9:      if e*<eiCi+1=C*ei+1=e*i=i+1ei+1<ei then

10:       Ci+1=C

11:       ei+1=e

12:      end if

13:      i=i+1

14:    end for

15:   While ei+1<ei

16:   Output Ci and ei

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 Tk-1,2Tk-1,,(k-2)Tk-1 respectively. Then, it uses k interval boundaries including t1 and tT 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 t1 and tT by average absolute difference with respect to its predecessor and successor time points by computing:

(2) mti=gG|Md(egti-1)-Md(egti)|+|Md(egti+1)-Md(egti)|2|G|

where Md(egti) is the median expression for gene g at time ti. We then select the k-2 points with maximum mti as the initial solution.

Finally, we developed an alternative initialization method, based on dyanmic recalculation of a metric on each time point. Metric A is same equal to the equation shown above. Metric B of a time point is the difference absolute difference with respect to its predecessor and successor time points. Metric C of a time point is absolute difference with respect to only its predecessor. The alternative initialization algorithm calculates the given metric on each time point other than the first and last and then places those points in a min heap based on the metric. The top(minimum) point in the heap is removed. The metric is recalculated for the point’s predeccesor and succesor based on thier neighboring points, using only the points remaining in the heap. This process is repated until only k-2 time points remain in the heap. Then the first time point, last time point and the points remaind in the heap are chosen.

(3) MetricAe,ti=gG|(Md(egpreviousti)-Md(egti))+(Md(egnextti)-Md(egti))|2|G|
(4) MetricBe,ti=gG|(Md(egpreviousti)-Md(egti))-(Md(egnextti)-Md(egti))|2|G|
(5) MetricCe,ti=gG|(Md(egpreviousti)-Md(egti))2|G|

Algorithm 2: Init TPS: Iterative initial k point selection

1: Procedure Iterative–Initial Point–Selection

2:  H= Empty min heap

3:  e= matrix where rows are genes and columns are time points, values are expression measurements

4:  for each time point t (other than the first and last) do

5:    valuet=Metrice,tprevioust=t-1nextt=t+1

6:    previoust=t1

7:    nextt=t+1

8:    Add valuet to Hsize(H)>k-2

9:  end for

10:  While size(H) > k2 do

11:   Remove minimum valuem time point m from Hpreviousnextm=previousmnextpreviousm=nextm

12:   previousnextm=previousm

13:   nextpreviousm=nextm

14:   Remove valuepreviousm from H

15:   Remove valuenextm from H

16:   Remove m from evaluepreviousm=Metrice,previousmvaluenextm=Metrice,nextm

17:   valuepreviousm=Metrice,previousm

18:   valuenextm=Metrice,nextm

19:   Add valuenextmtoH

20:   Add valuepreviousmtoH

21:  end while

22:   Ouput all t left as valuet in H + first time point + last time point

23: end procedure

We found that for our particular dataset, the dynamic initialization with MetricAe,ti performed best for selections of time points smaller than one third of the the initial dense time series, while the non dynamic mti 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(kGT2Q) to O(kGT2bQ) 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 er is lower than error of current solution ei whereas we move to a solution with probability e-C(er-ei) if erei. 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.

Fitting smoothing spline

TPS uses splines for fitting expression curves. Regularized smoothing spline satisfies the piecewise cubic polynomial μ(t)=ai+bi(t-ti)+ci(t-ti)2+di(t-ti)3 for t[ti,ti+1),i1,,T-1 as shown in Wahba (1990). Then, according to (Reinsch, 1967; De Boor, 1978), regularized smoothing spline objective can also be expressed as in:

(6) min(y-a)(y-a)+λcRc

where a=(a1,a2,,aT), c=(c2,c3,,cT-1), and R is a (n-2)2 tridiagonal symmetric matrix with entries ri,i=2(hi+hi+1)3, ri,i+1=hi+13 where hi=ti+1-ti. The continuity restrictions imply that:

(7) Rc=Qa

where Q is an n×(n-2) tridiagonal matrix with entries qi,i+1=1hi+1, qi+1,i=1hi+1 and qi,i=-(1hi+1hi+1). Thus, we may write Equation 6 as:

(8) min(y-a)(y-a)+λaQR-1Qa

where a can be derived as in:

(9) a=(I+λQR-1Q)-1y

Once a is estimated, b, c, d are estimated by corresponding Equations in Reinsch (1967).

For our specific setting, we also introduce a regularization parameter to enable us to determine the number of control points. Let Ig={(t,Md(egt)),tC}, 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:

(10) min(t,yt)Ig(yt-μ(t))2+λt1tTμ′′(x)2dx

where λ is the regularization parameter which prevents overfitting by affecting the number of knots selected. We estimated λ 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.

Appendix 2 Results

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 2—figure 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 single-organism cellular process, positive regulation of biological process, regulation of metabolic process, etc. See Supporting Website for complete results.

miRNA reconstruction

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. Several of these miRNAs are known to be involved in regulation of lung development. For example, mmu-miR-100 is known to regulate Fgfr3 and Igf1r, 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 ( 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 ( 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 non-optimized 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 ( 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.

Appendix 2—table 1

List of genes used for the Nanostring analysis and the rational for their inclusion.
Ensembl gene IDAccession numberGene nameRationale
ENSMUSG00000024130NM_001039581.2Abca3Alveolar Type II cell marker
ENSMUSG00000031378NM_007435.1Abcd1important in other processes (IPF, COPD etc)
ENSMUSG00000029802NM_011920.3Abcg2Mesenchymal cell marker
ENSMUSG00000035783NM_007392.3Acta2Fibroblast cell marker
ENSMUSG00000029580NM_007393.1ActbCommon house-keeping gene
ENSMUSG00000036040NM_029981.1Adamtsl2Altered DNA methylation during septation
ENSMUSG00000015452NM_007425.2AgerAlveolar Type I cell marker
ENSMUSG00000001729NM_001165894.1Akt1Altered DNA methylation during septation
ENSMUSG00000053279NM_013467.3Aldh1a1Important for septation
ENSMUSG00000013584NM_009022.3Aldh1a2Potentially important for septation
ENSMUSG00000022244NM_008537.4Amacrimportant in other processes (IPF , COPD etc)
ENSMUSG00000044217NM_009701.4Aqp5Alveolar Type I cell marker
ENSMUSG00000026576NM_009721.5Atp1b1Lung fluid clearance
ENSMUSG00000060802NM_009735.3B2mCommon house-keeping gene
ENSMUSG00000102037NM_009742.3Bcl2a1aApoptosis regulator
ENSMUSG00000056216NM_009884.3CebpgImportant for lung development
ENSMUSG00000029084NM_007646.4Cd38Airway smooth muscle cell functional responses
ENSMUSG00000018774NM_009853.1Cd68Monocyte cell marker
ENSMUSG00000031673NM_009866.4Cdh1Epithelial cell marker
ENSMUSG00000064246NM_007695.2Chil1Monocyte cell marker
ENSMUSG00000040809NM_009892.1Chil3Increased during septation
ENSMUSG00000022512NM_016674.3Cldn1Tight junction protein
ENSMUSG00000070473NM_009902.4Cldn3Tight junction protein (mostly epithelial)
ENSMUSG00000041378NM_013805.4Cldn5Tight junction protein
ENSMUSG00000018569NM_016887.6Cldn7Tight junction protein (mostly epithelial)
ENSMUSG00000001506NM_007742.3Col1a1Fibroblast cell marker
ENSMUSG00000063063NM_009819.2Ctnna2Altered DNA methylation during septation
ENSMUSG00000031360NM_001168571.1Ctps2important in other processes (IPF , COPD etc)
ENSMUSG00000040856NM_010052.4Dlk1Decreased during septation
ENSMUSG00000020661NM_007872.4Dnmt3aAltered DNA methylation during septation
ENSMUSG00000046179NM_001013368.5E2f8Altered DNA methylation during septation
ENSMUSG00000000303NM_009864.2Cdh1Epithelial cell marker
ENSMUSG00000020122NM_207655.2EgfrImportant for lung development
ENSMUSG00000029675NM_007925.3ElnAltered DNA methylation during septation
ENSMUSG00000045394NM_008532.2EpcamEpithelial cell marker
ENSMUSG00000052504NM_010140.3Epha3Involved in lung development
ENSMUSG00000028289NM_001122889.1Epha7Involved in lung cancer, potential role in development
ENSMUSG00000021055NM_010157.3Esr2Important regulator of multiple processes
ENSMUSG00000061731NM_010162.2Ext1Altered DNA methylation during septation
ENSMUSG00000039109NM_001166391.1F13a1Involved in lung injury , cancer
ENSMUSG00000057967NM_008005.1Fgf18Important for septation
ENSMUSG00000030849NM_010207.2Fgfr2Important regulator of multiple processes
ENSMUSG00000078302NM_008242.2Foxd1Pericyte cell marker
ENSMUSG00000042812NM_010426.1Foxf1Involved in lung development
ENSMUSG00000038402NM_010225.1Foxf2Altered DNA methylation during fibrosis
ENSMUSG00000001020NM_011311.1S100a4Fibroblast cell marker
ENSMUSG00000057666NM_001001303.1GapdhCommon house-keeping gene
ENSMUSG00000005836NM_010258.3Gata6Important regulator of multiple processes
ENSMUSG00000029992NM_013528.3Gfpt1important in other processes (IPF, COPD etc)
ENSMUSG00000041624NM_001033322.2Gucy1a2Important for septation
ENSMUSG00000025534NM_010368.1GusbCommon house-keeping gene
ENSMUSG00000021109NM_010431.2Hif1aHypoxia signaling
ENSMUSG00000058773NM_020034.1Hist1h1bDecreased during septation
ENSMUSG00000061615NM_175660.3Hist1h2abDecreased during septation
ENSMUSG00000032126NM_013551.2HmbsCommon house-keeping gene
ENSMUSG00000029919NM_019455.4Hpgdsimportant in other processes (IPF, COPD etc)
ENSMUSG00000025630NM_013556.2HprtCommon house-keeping gene
ENSMUSG00000020053NM_001111274.1Igf1Regulating miRNA altered during septation
ENSMUSG00000020427NM_008343.2Igfbp3Altered DNA methylation during septation, fibrosis
ENSMUSG00000003477NM_009349.3InmtIncreased during septation
ENSMUSG00000026768NM_001001309.2Itga8Involved in lung development
ENSMUSG00000040029NM_001081113.1Ipo8important in other processes (IPF, COPD etc)
ENSMUSG00000030786NM_001082960.1ItgamMonocyte cell marker
ENSMUSG00000030789NM_021334.2ItgaxMonocyte cell marker
ENSMUSG00000090122NM_021487.1Kcne1limportant in other processes (IPF, COPD etc)
ENSMUSG00000063142.10XM_006518608.1Kcnma1Altered DNA methylation during septation
ENSMUSG00000079852NM_010649.3Klra4Increased during septation
ENSMUSG00000023043NM_010664.2Krt18Epithelial cell marker
ENSMUSG00000061527NM_027011.2Krt5Basal cell marker
ENSMUSG00000029570NM_008494.3LfngImportant for septation
ENSMUSG00000024529NM_010728.2LoxAltered DNA methylation during fibrosis
ENSMUSG00000028003NM_023624.4LratIncreased during septation
ENSMUSG00000027070NM_001081088.1Lrp2Altered DNA methylation during septation
ENSMUSG00000061068NM_010779.2Mcpt4Decreased during septation
ENSMUSG00000026110NM_173870.2Mgat4aInvolved in acute lung injury
ENSMUSG00000043613NM_010809.1Mmp3Increased during septation
ENSMUSG00000018623NM_010810.4Mmp7Important in lung fibrosis
ENSMUSG00000066108XM_006508653.1Muc5bImportant in lung fibrosis
ENSMUSG00000037974NM_010844.1Muc5acEpithelial cell marker
ENSMUSG00000024304NM_007664.4Cdh2Tight Junction/Adhesion
ENSMUSG00000054008NM_008306.4Ndst1Involved in pathologic airway remodeling
ENSMUSG00000031902NM_010901.2Nfatc3Important for lung development
ENSMUSG00000073435NM_019730.2Nme3Apoptosis-related gene
ENSMUSG00000026575NM_138314.3Nme7Important for stem cell renewal
ENSMUSG00000014776NM_030152.4Nol3Regulating miRNA altered during septation
ENSMUSG00000051048NM_177161.4P4ha3Important in lung fibrosis
ENSMUSG00000068039NM_013686.3Tcp1Basal cell marker
ENSMUSG00000029998NM_025823.4Pcyox1important in other processes (IPF , COPD etc)
ENSMUSG00000029231NM_011058.2PdgfraImportant for septation
ENSMUSG00000024620NM_008809.1PdgfrbPericyte cell marker
ENSMUSG00000028583NM_010329.2PdpnAlveolar Type I cell marker
ENSMUSG00000062070NM_008828.2Pgk1important in other processes (IPF , COPD etc)
ENSMUSG00000053398NM_016966.3Phgdhimportant in other processes (IPF, COPD etc)
ENSMUSG00000005198NM_009089.2Polr2aimportant in other processes (IPF, COPD etc)
ENSMUSG00000071866NM_008907.1PpiaCommon house-keeping gene
ENSMUSG00000024997NM_007452.2Prdx3Mitochondrial oxidative stress regulator
ENSMUSG00000026134NM_008922.2Prim2Expressed in placenta and crucial for mammalian growth.
ENSMUSG00000033491NM_178738.3Prss35Decreased during septation
ENSMUSG00000032487NM_011198.3Ptgs2Regulating miRNA altered during septation
ENSMUSG00000056458NM_011973.2MokAlveolar Type I cell marker
ENSMUSG00000037992NM_001177302.1RaraImportant for septation
ENSMUSG00000022883NM_019413.2Robo1Altered DNA methylation during septation
ENSMUSG00000066361NM_008458.2Serpina3cIncreased during septation
ENSMUSG00000022097NM_011359.1SftpcAlveolar Type II cell marker
ENSMUSG00000021795NM_009160.2SftpdAlveolar Type II cell marker
ENSMUSG00000050010NM_001033415.3Shisa3Altered DNA methylation during septation
ENSMUSG00000032402NM_016769.3Smad3Important for septation
ENSMUSG00000042821NM_011427.2Snai1Important for lung development and injury
ENSMUSG00000000567NM_011448.4Sox9Altered DNA methylation during septation
ENSMUSG00000027646NM_001025395.2SrcAltered DNA methylation during septation
ENSMUSG00000014767NM_013684.3TbpCommon house-keeping gene , involved in multiple processes
ENSMUSG00000000094NM_172798.1Tbx4Altered DNA methylation during septation
ENSMUSG00000032228NM_011544.3Tcf12Involved in multiple developmental processes
ENSMUSG00000022797NM_011638.3TfrcCommon house-keeping gene
ENSMUSG00000002603NM_011577.1Tgfb1Important for septation
ENSMUSG00000045691NM_153083.5Thtpaimportant in other processes (IPF, COPD etc)
ENSMUSG00000032011NM_009382.3Thy1Fibroblast cell marker
ENSMUSG00000028364NM_011607.1TncAltered DNA methylation during septation
ENSMUSG00000044986NM_009437.4Tstimportant in other processes (IPF, COPD etc)
ENSMUSG00000026803NM_009442.2Ttf1Important for lung development
ENSMUSG00000008348NM_019639.4UbcCommon house-keeping gene
ENSMUSG00000023951NM_001025250.3VegfaAngiogenesis; Altered DNA methylation during septation
ENSMUSG00000026728NM_011701.4VimMesenchymal cell marker
ENSMUSG00000020218NM_011915.1Wif1Altered DNA methylation during septation
ENSMUSG00000022285NM_011740.2YwhazCommon house-keeping gene
Appendix 2—table 2

Summary of methylation dataset
GeneNumber of lociGeneNumber of loci
Appendix 2—table 3

Target regions for each gene for methylation analysis
GeneEnsembl gene IDEnsembl transcript IDAssay IDTarget locationFwd TmRev Tm% GCCoordinates (GRCm38/mm10)
Akt1ENSMUSG00000001729ENSMUST00000001780ADS33333’ UTR6865.531.5chr12:112654548–112654709
Akt1ENSMUSG00000001729ENSMUST00000001780ADS3332Intron 9/Exon 1068.369.838.3chr12:112657120–112657273
Cdh11ENSMUSG00000031673ENSMUST00000075190ADS3308Intron 366.868.336.9chr8:102677609–102677766
Cdh11ENSMUSG00000031673ENSMUST00000075190ADS3318Intron 164.169.737chr8:102784569–102784722
Dnmt3aENSMUSG00000020661ENSMUST00000020991ADS632Intron 16464.732.2chr12:3834382–3834592
Dnmt3aENSMUSG00000020661ENSMUST00000020991ADS3328Exon 6/Intron 664.76431.8chr12:3901545–3901764
Dnmt3aENSMUSG00000020661ENSMUST00000020991ADS3329Intron 666.866.125.4chr12:3907514–3907765
ElnENSMUSG00000029675ENSMUST00000015138ADS3319Intron 1667.167.947.8chr5:134721191–134721447
ElnENSMUSG00000029675ENSMUST00000015138ADS3309Intron 7/Exon 8/Intron 864.167.437.8chr5:134729221–134729526
Foxf2ENSMUSG00000038402ENSMUST00000042054ADS4505Promoter63.16542.5chr13: 31625470–31625556
Igfbp3ENSMUSG00000020427ENSMUST00000020702ADS3301Exon 4/Intron 470.57033chr11:7208306–7208481
Igfbp3ENSMUSG00000020427ENSMUST00000020702ADS5133Intron 168.368.526.1chr11:7212803–7213043
LoxENSMUSG00000024529ENSMUST00000171470ADS4512Exon 26970.931.3chr18: 52529184–52529315
LoxENSMUSG00000024529ENSMUST00000171470ADS4513Exon 465.764.828.5chr18:52526887–52527023
Sox9ENSMUSG00000000567ENSMUST00000000579ADS3311Intron 169.768.534.7chr11:112783358–112783605
Sox9ENSMUSG00000000567ENSMUST00000000579ADS3310Exon 366.463.126.2chr11:112784760–112784885
SrcENSMUSG00000027646ENSMUST00000109533ADS4514Intron 164.865.935.9chr2:157423925–157424027
SrcENSMUSG00000027646ENSMUST00000109533ADS4515Intron 466.568.837.6chr2:157457351–157457520
SrcENSMUSG00000027646ENSMUST00000109533ADS4516Exon 1465.565.633.7chr2:157469741–157469912
TncENSMUSG00000028364ENSMUST00000107377ADS3324Intron 1463.362.223chr4:63982645–63982818
TncENSMUSG00000028364ENSMUST00000107377ADS3325Intron 1462.561.620.2chr4:63982799–63982986
TncENSMUSG00000028364ENSMUST00000107377ADS3323Exon 36567.535.2chr4:64017478–64017721
VegfaENSMUSG00000023951ENSMUST00000071648ADS3335Intron 2/Exon 364.663.829.5chr17:46025336–46025620
Wif1ENSMUSG00000020218ENSMUST00000020439ADS3303Intron 4/Exon 5/Intron 560.960.131.8chr10:121083800–121083997
Wif1ENSMUSG00000020218ENSMUST00000020439ADS3304Exon 10/3-UTR66.667.524.8chr10:121099752–121099973
Zfp536ENSMUSG00000043456ENSMUST00000056338ADS4509Exon 468.469.935.4chr7:37567973–37568130
Appendix 2—table 4

Mean and standard deviation of mean squared error over all 126 genes by TPS selecting 5 points and piecewise linear fits over 3 sets of points identified heuristically in the literature.
MethodMeanStd dev
TPS (0.5, 6, 9.5, 19 and 28)0.403063359620.2206665163
Piecewise linear over 0.5, 7, 14, 280.5940727194940.399642079492
Piecewise linear over 0.5, 2, 14, 280.7109670613490.721681860787
Piecewise linear over 0.5, 4, 7, 14, 280.5609902305010.364739525724

Data availability

The following data sets were generated
    1. Kleyman M
    2. Sefer E
    3. Nicola T
    4. Espinoza C
    5. Chhabra D
    6. Hagood JS
    7. Kaminski N
    8. Ambalavanan N
    9. Joseph ZB
    (2016) miRNA Data of Mouse Lung Developement
    Publicly available at the Systems Biology Group, School of Computer Science, Carnegie Mellon University website.


  1. Book
    1. Bishop CM
    Pattern recognition and machine learning
    1. Cormack M
    2. Lin M
    3. Fedulov A
    4. Mentzer SJ
    5. Tsuda A
    Age-dependent changes in gene expression profiles of postnatally developing rat lungs exposed to nano-size and micro-size cuo particles
    The FASEB Journal 24:612–618.
    1. De Boor C
    Mathematics of Computation
    A practical guide to splines, Mathematics of Computation, New York, Springer, NY.
  2. Book
    1. Hartigan JA
    Clustering Algorithms
    Hoboken, NJ: John Wiley & Sons, Inc.
  3. Book
    1. McLachlan G
    2. Peel D
    Finite Mixture Models
    Hoboken, NJ: John Wiley & Sons.
  4. Conference
    1. Singh R
    2. Palmer N
    3. Gifford D
    4. Berger B
    5. Bar-Joseph Z
    Active learning for sampling in time-series experiments with application to gene expression analysis
    Proceedings of the 22nd International Conference on Machine Learning (ICML-05). pp. 832–839.
  5. Book
    1. Wahba G
    Spline Models for Observational Data, Vol. 59
    Philadelphia, PA: SIAM: Society for Industrial and Applied Mathematics.
    1. Weinstein M
    2. Xu X
    3. Ohyama K
    4. Deng CX
    FGFR-3 and FGFR-4 function cooperatively to direct alveogenesis in the murine lung
    Development 125:3615–3623.

Article and author information

Author details

  1. Michael Kleyman

    Machine Learning and Computational Biology, School of Computer Science, Carnegie Mellon University, Pittsburgh, United States
    MK, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Emre Sefer
    Competing interests
    The authors declare that no competing interests exist.
  2. Emre Sefer

    Machine Learning and Computational Biology, School of Computer Science, Carnegie Mellon University, Pittsburgh, United States
    ES, Conceptualization, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft
    Contributed equally with
    Michael Kleyman
    Competing interests
    The authors declare that no competing interests exist.
  3. Teodora Nicola

    Division of Neonatology, Department of Pediatrics, University of Alabama at Birmingham, Birmingham, United States
    TN, Resources, Data curation, Investigation
    Competing interests
    The authors declare that no competing interests exist.
  4. Celia Espinoza

    1. Division of Respiratory Medicine, Department of Pediatrics, University of California, San Diego, United States
    2. CARady Children’s Hospital San Diego, San Diego, United States
    CE, Resources, Data curation, Investigation
    Competing interests
    The authors declare that no competing interests exist.
  5. Divya Chhabra

    1. Division of Respiratory Medicine, Department of Pediatrics, University of California, San Diego, United States
    2. CARady Children’s Hospital San Diego, San Diego, United States
    DC, Resources, Data curation, Investigation
    Competing interests
    The authors declare that no competing interests exist.
  6. James S Hagood

    1. Division of Respiratory Medicine, Department of Pediatrics, University of California, San Diego, United States
    2. CARady Children’s Hospital San Diego, San Diego, United States
    JSH, Resources, Data curation, Investigation
    Competing interests
    The authors declare that no competing interests exist.
  7. Naftali Kaminski

    Section of Pulmonary, Critical Care and Sleep Medicine, School of Medicine, Yale University, New Haven, United States
    NK, Resources, Data curation, Investigation
    Competing interests
    The authors declare that no competing interests exist.
  8. Namasivayam Ambalavanan

    Division of Neonatology, Department of Pediatrics, University of Alabama at Birmingham, Birmingham, United States
    NA, Resources, Data curation, Investigation
    Competing interests
    The authors declare that no competing interests exist.
  9. Ziv Bar-Joseph

    Machine Learning and Computational Biology, School of Computer Science, Carnegie Mellon University, Pittsburgh, United States
    ZB-J, Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3430-6051


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.


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.


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.

Version history

  1. Received: June 6, 2016
  2. Accepted: January 23, 2017
  3. Accepted Manuscript published: January 26, 2017 (version 1)
  4. Version of Record published: February 21, 2017 (version 2)


© 2017, Kleyman et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 2,508
  • 494
  • 28

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Michael Kleyman
  2. Emre Sefer
  3. Teodora Nicola
  4. Celia Espinoza
  5. Divya Chhabra
  6. James S Hagood
  7. Naftali Kaminski
  8. Namasivayam Ambalavanan
  9. Ziv Bar-Joseph
Selecting the most appropriate time points to profile in high-throughput studies
eLife 6:e18541.

Share this article

Further reading

    1. Cell Biology
    2. Computational and Systems Biology
    N Suhas Jagannathan, Javier Yu Peng Koh ... Lisa Tucker-Kellogg
    Research Article

    Bats have unique characteristics compared to other mammals, including increased longevity and higher resistance to cancer and infectious disease. While previous studies have analyzed the metabolic requirements for flight, it is still unclear how bat metabolism supports these unique features, and no study has integrated metabolomics, transcriptomics, and proteomics to characterize bat metabolism. In this work, we performed a multi-omics data analysis using a computational model of metabolic fluxes to identify fundamental differences in central metabolism between primary lung fibroblast cell lines from the black flying fox fruit bat (Pteropus alecto) and human. Bat cells showed higher expression levels of Complex I components of electron transport chain (ETC), but, remarkably, a lower rate of oxygen consumption. Computational modeling interpreted these results as indicating that Complex II activity may be low or reversed, similar to an ischemic state. An ischemic-like state of bats was also supported by decreased levels of central metabolites and increased ratios of succinate to fumarate in bat cells. Ischemic states tend to produce reactive oxygen species (ROS), which would be incompatible with the longevity of bats. However, bat cells had higher antioxidant reservoirs (higher total glutathione and higher ratio of NADPH to NADP) despite higher mitochondrial ROS levels. In addition, bat cells were more resistant to glucose deprivation and had increased resistance to ferroptosis, one of the characteristics of which is oxidative stress. Thus, our studies revealed distinct differences in the ETC regulation and metabolic stress responses between human and bat cells.

    1. Computational and Systems Biology
    2. Neuroscience
    Sara Ibañez, Nilapratim Sengupta ... Christina M Weaver
    Research Article

    Normal aging leads to myelin alterations in the rhesus monkey dorsolateral prefrontal cortex (dlPFC), which are positively correlated with degree of cognitive impairment. It is hypothesized that remyelination with shorter and thinner myelin sheaths partially compensates for myelin degradation, but computational modeling has not yet explored these two phenomena together systematically. Here, we used a two-pronged modeling approach to determine how age-related myelin changes affect a core cognitive function: spatial working memory. First, we built a multicompartment pyramidal neuron model fit to monkey dlPFC empirical data, with an axon including myelinated segments having paranodes, juxtaparanodes, internodes, and tight junctions. This model was used to quantify conduction velocity (CV) changes and action potential (AP) failures after demyelination and subsequent remyelination. Next, we incorporated the single neuron results into a spiking neural network model of working memory. While complete remyelination nearly recovered axonal transmission and network function to unperturbed levels, our models predict that biologically plausible levels of myelin dystrophy, if uncompensated by other factors, can account for substantial working memory impairment with aging. The present computational study unites empirical data from ultrastructure up to behavior during normal aging, and has broader implications for many demyelinating conditions, such as multiple sclerosis or schizophrenia.