Author response: Selecting the most appropriate time points to profile in high-throughput studies (2016)

Article Figures and data Abstract Introduction Results Discussion Materials and methods Appendix 1 Methods Appendix 2 Results Data availability References Decision letter Author response Article and author information Metrics Abstract 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 https://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 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. 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. Figure 1 with 3 supplements see all Download asset Open asset 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. https://doi.org/10.7554/eLife.18541.002 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 Download asset Open asset 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. https://doi.org/10.7554/eLife.18541.006 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 Download asset Open asset Reconstructed expression profiles for selected genes. (a). Pdgfra. , (b). Eln. , (c). Inmt. https://doi.org/10.7554/eLife.18541.011 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 Download asset Open asset 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. https://doi.org/10.7554/eLife.18541.014 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 Download asset Open asset Comparison of gene expression and methylation data for selected genes. (a). Akt1. , (b). Cdh11. , (c). Tnc. https://doi.org/10.7554/eLife.18541.018 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. Table 1 Summary of prior high throughput lung development studies. https://doi.org/10.7554/eLife.18541.021 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 levelsE12, 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, et al., mRNA expression E18, 1, 4, 7, 14, adult To these 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 et are 0, 4, 7, and 28 28 is last day in our analysis we used it of In TPS selects 6, and As can be in Figure 6, important expression changes in key genes are missed by using the points while the TPS points are able to reconstruct these profiles even though the number of points is the same More the error for the set of selected points is much higher on average results are obtained for the other sampling rates used in the 6, Appendix and when TPS to iterative methods for selecting the set of points to profile supplement 1). This that accurate selection of time points can have a large on the ability of the study to identify key genes and events. See also Appendix Results for a the importance of the differences between the TPS and prior work results for selected genes. Figure with 1 supplement see all Download asset Open asset Comparison of TPS with sampling rates used in studies. curves are the reconstructed profiles based on the points profiled by prior studies. and curves are based on the points selected by TPS As can be seen, even when results from using the same number of points, TPS can identify key events for some of the genes that are missed when using the based sampling and are a fit over points and are a fit over points and are a fit over points Our method on a very small subset of genes that are known to be involved in the process studied for the initial set of experiments. While such set is known for several may be cases very little is known the biological process and so it may be hard to obtain such TPS can still be applied to determine sampling rates for such processes using a small random set of genes. To this we repeated the analysis presented in Results using only the measured values of of genes in our set and the values for the other genes with random As we show in Figure 2—figure supplement 2, even when using such the time points selected by TPS greatly improve upon an set of the same number of time points. in most time series experiments at of the genes are differentially expressed in several cases a much larger et al., et al., 2015) a random selection of genes is to similar results even for processes. the analysis of a specific type of data, several studies have now been profiling multiple types of genomic data over time. studies to on a set of time points which would be common to all experiments so that these types can be integrated to a model (Chang et al., 2013; et al., 2010). To date, the selection of such points relied on ad-hoc The processes being studied were either sampled or based on prior However, known of such systems were often been based on phenotypic observations which may not with the of molecular events. In addition, in many case studies of the or similar processes with to the time points that have been profiled. For example, work on the analysis of cell data in utilized uniform and nonuniform sampling et al., and studies of have a similar et al., et al., more analysis of responses to widely in the sampling rates that were used et al., et al., TPS these by using a principled method for determining sampling 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 gene method provides a very good for expression profiling as well as other types of data including miRNAs and DNA Thus, a experiment can be using our method. While we TPS on several types of high throughput data, we have only tested it so far on data for a specific biological process development in mice). While we that such data is and representative and thus provides a good test case for the analysis of additional datasets may identify new that we have not and we it to work to address including all initialization methods is in and is available on the We that as sequencing technology to more and more studies would integrate types of time series data and TPS in the design of their studies. Materials and methods mRNA and miRNA used in the study a To select the list of 126 genes used in the NanoString profiling we the for genes that have been to the (a) type genes alveolar type alveolar type any (b) genes known to be or during (c) genes known to be in DNA methylation during development, (d) genes known to be involved in genes known to be by miRNA involved in and (f) genes known to be by DNA methylation during Appendix 1 a list of the selected genes and the process for which they were For the miRNA set we used a unbiased, miRNA Expression mRNA and miRNA profiling and analysis a of samples were by from lung at multiple time points to and to The samples were used to was performed by the and were measured by using and NanoString set set and for 126 genes was and the Gene Expression was performed using The data by the were as a and data normalization was performed using the the analysis by DNA methylation analysis a alveolar lung to were at until DNA was using the