Testing DTW-S using heterochrony simulations
To test the efficiency of time shift estimation by DTW-S, we first applied the algorithm to a number of simulated datasets with known heterochrony.
Sine function
We tested whether our method could be applied to periodic expression changes, such as expression changes associated with the cell cycle or circadian rhythm. Across the trials, we used the same error distribution parameters and time shift, while we varied the three parameters n, N, and M (these parameters reflect numbers of actual and interpolated time points, see Methods for Details). For simplicity, we assumed numbers of actual and interpolated time points to be equal (n = N). The expression changes with time were modeled as y = sin(π* 25/t)+ε, with error ε~N(0,0.3) and fixed time shift Δ = 5 for all time points. Using these error distribution parameters, residual variance constituted approximately 15% of the total variance in the simulated datasets. This percentage of error-related variance is within the range of technical and biological variance, relative to factors such as age and species identity, observed in actual microarray experiments (e.g. [15]). These results demonstrate that DTW-S estimates are robust with respect to the numbers of original data points and interpolated time points. Furthermore, DTW-S could estimate the true time shift value (Figure 2). As expected, variation among time shift estimates decreased with increasing number of time points.
Linear and polynomial functions
We tested DTW-S on non-periodic gene expression changes, such as changes described by linear and polynomial curves. Most organismal processes, e.g. development, aging, or response to stress, involve such non-periodic patterns. Specifically, we tested the algorithm's ability to determine constant or variable time shifts, as well as it's performance at different levels of error variance (ε) and with a different sample size. We tested the algorithm using the following parameters: n = 20, N = 20 and M = 40.
First, we applied DTW-S on simulated time series based on linear and polynomial curves with a constant shift. We used the function y = a+bt+dt2+ε, with a = (0,1), b = (-1,0,1), d = (-1,0,1), to simulate expression time series. In addition, we used a fixed time shift Δ = 2 for all time points. We set the error distribution parameters at ~10% or ~20% of the total variance, i.e. 10% or 20% of the total variance is attributable to error. This proportion is within the error variance range observed in actual microarray experiments. Under both error rates and for both linear and polynomial curves, our method yielded accurate and robust time shift estimates (Additional file 1, Figures S2, S3).
Second, we tested the ability of our method to determine variable time shifts between two time series. In a biological system, it can be reasonably assumed that the difference in timing of expression changes between two processes might increase or decrease gradually. We therefore used four sets of variant time shift (the linear shifts: C1 and C2, and the polynomial shifts: C3 and C4) to simulate expression time series, two linear and two quadratic (y = 3+2t-t2+ε and y = 3-2t+t2+ε), with error ε~Ν(0,4), such that residual variance constituted ~20% of the total variance. Our results showed that under these conditions DTW-S can effectively identify the modeled variable time shift between two time series (Figure 3). Specifically, in simulations for non-linear shift patterns (C3 and C4 in Figure 3), the predicted time shifts were significantly better described by the original polynomial shift models, than they were by linear models (ANOVA p < 0.01).
Next, we assigned the same four sets of variant time shift (C1, C2, C3, and C4) to time series, modeled using combinations of the expression profile parameters a, b, d (i.e. linear and polynomial curves), with ~20% residual variance. We found that, irrespective of the expression curve type and the time shift type, our method robustly identified the original shift profile between time series (Additional file 1, Figure S4).
Lastly, we tested our method's performance under 30% residual variance. Overall, DTW-S performed well at this relatively high error rate. Compared to the original time shift, the estimated time shifts showed higher variance at the middle of the age distribution, and biases towards positive and negative time shift estimates at the beginning and end of the age distribution respectively (Additional file 1, Figure S5). We note that, under even higher levels of noise, DTW-S might not be able to distinguish a non-linear time shift trajectory from a linear one. In these cases, it might still be possible to accurately estimate the magnitude and type of shift by averaging across multiple genes to minimize noise.
Application to primate brain development
Human and chimpanzee cortical development
We first applied the DTW-S algorithm to investigate the heterochrony patterns between the brain developmental profiles of humans and chimpanzees. Prevalence of a specific type of heterochrony, delayed development in humans (or neoteny), was previously reported in human and chimpanzee brains [15]. This result was based on expression time series measured in the dorsolateral prefrontal cortex in 39 humans (0-47 years), 14 chimpanzees (0-44 years) and one outgroup species: 9 rhesus macaques (1-18 years). For all three species, the dataset covered most of the postnatal ontogenesis and maturation, from birth until adulthood (Additional file 1, Figure S1A). In this study, we re-examined this work using DTW-S. Our hypothesis was that expression heterochrony patterns might vary among genes, even among those that show the same type of heterochrony, with regard to the timing of the shift, as well as shift magnitude. Following the original analysis, we first selected 1183 genes satisfying the following criteria: (1) significant expression changes with age (using polynomial regression models, at F-test p < 0.05), (2) significant expression difference between humans and chimpanzees (using analysis of covariance (ANCOVA), at F-test p < 0.05) and, (3) significant positive correlation between human and chimpanzee expression profiles (Pearson correlation p < 0.05 and r > 0). We then applied the DTW-S algorithm to these 1183 genes to identify expression trajectories with significant time shift between the two species. Here, we aligned the chimpanzee time series to the human age-scale. For each gene, we set the p-value cutoff for the significant time shift points at 5%, and the gene significance cutoff as more than one third of the gene's time points having significant time shift (Methods). At these cutoffs, 482 out of the 1183 genes tested had significant time shift, and the FPR was estimated at 10.7% (Figure 4A).
We assigned these genes to one of four phylo-ontogenetic categories, using our previously reported classifications [15]: (a) human-specific neotenic expression (here neoteny denotes a slow-down or delay in development), (b) human-specific accelerated expression, (c) chimpanzee-specific neotenic expression and, (d) chimpanzee-specific accelerated expression (Additional file 1, Figure S6). To do so, we first determined the direction of the time shift (acceleration, neoteny, or none) between the human and chimpanzee time series, for all 482 genes with significant time shift. If at least 70% of time points showed consistent shift direction, we classified the related gene as accelerated or neotenic. We then determined whether expression change has happened to the human or to the chimpanzee evolutionary lineage, using the rhesus macaque time series as the outgroup (Methods). For all 482 genes, we could assign 337 (at lineage assignment test p < 0.05), or 267 (at lineage assignment test p < 0.01), to one of the four phylo-ontogenetic categories. This analysis yielded a clear excess of genes in the human-specific neotenic expression category (Figure 4B), which adds support to our previously published results [15].
To test the robustness of our method, we reversed the entire time shift identification procedure by aligning the human time series to the chimpanzee trajectory. This yielded 642 genes with significant shift (FPR = 8.7%, p-value < 0.05, number of significant time shift points per gene = 7). The difference in the numbers of genes with significant time shift compared to the previous procedure, is most likely related to sample size differences between the chimpanzee and the human time series (Methods). Despite this variation, the two sets of genes identified in the two procedures largely overlap: 369 between 482 and 642 genes (Fisher exact test: p < 0.0001). The proportions of genes assigned to the four categories were also highly consistent (Figure 4B).
In the previous study, any differences in time shift magnitude among the heterochronic genes, were ignored. Furthermore, the time shift was assumed to remain constant across all time points [15]. To test these assumptions, we analyzed the time shift patterns estimated by DTW-S for 118 human-specific neotenic genes. Using hierarchical clustering we found that, instead of forming a single pattern, the time shift profiles of these genes fell into three distinct time shift clusters (Additional file 1, Figure S7). Clustering time shift profiles using the k-means algorithm also produced consistent segregation into three clusters. In 828 out of 1,000 iterations, the k-means algorithm yielded the same distribution of time shift patterns across the three clusters. The three clusters also showed marked differences in both the timing and the amplitude of the time shift (Figure 4C). Specifically, all three clusters showed time shift peaks at 20-30 years of age in humans, but at different ages in the chimpanzees: 1 year (Cluster2), 2-4 years (Cluster3), and 4-6 years (Cluster1) (Figure 4C).
To test the robustness of the time shift estimates for the three clusters, we combined the estimates obtained in the chimpanzee-to-human alignment, and in the human-to-chimpanzee alignment, and repeated the clustering using the k-means algorithm. This approach also led to robust segregation into three clusters. In addition, for the 170 genes with significant time shift estimates in at least one of the two alignments, both direct and reverse alignments produced consistent time shift estimates between human and chimpanzee expression profiles within each of the three clusters (Figure 4D).
Finally, we investigated whether the various patterns of heterochrony identified by our approach may be associated with specific biological functions. We first tested the characteristics of genes in the three time shift clusters using published data for gene expression specific to human brain gray and white matter [21], as well as on cell type-specific expression in the mouse central nervous system [22]. Genes in Cluster1 were significantly enriched among genes preferentially expressed in human gray matter (hypergeometric test [HT] p = 0.001), while genes in the other two clusters did not show any cell type specificity ([HT] p > 0.5) (Figure 4E). Consistently, genes in Cluster1, but not in the other two clusters, tended to be enriched among neuron-specific genes ([HT] p = 0.058) (Figure 4E). These results indicate that differences in time shift patterns, identified by DTW-S, might reflect differences in ontogenetic timing of gene expression changes among different human brain cell types, and between histological locations.
Development of the human and macaque prefrontal cortex and cerebellum
To further test the utility of the DTW-S in uncovering novel biological features, we studied time shift patterns in human and rhesus macaque brain development in two distinct regions: prefrontal cortex and cerebellum. For this, we used a published expression dataset containing time series measured in the dorsolateral prefrontal cortexof 23 humans and 26 rhesus macaques, and in the cerebellum of 22 humans and 24 rhesus macaques [23]. In both species, the dataset covered most of the lifespan, humans: 0-98 years, macaques: 0-28 years (Additional file 1, Figure S1B).
Based on these data, we selected 1084 genes in the cortex and 950 genes in the cerebellum that satisfied the following criteria: (1) significant expression changes with age (using polynomial regression models, at F-test p < 0.05), (2) significant expression differences between humans and rhesus macaques (using analysis of covariance (ANCOVA), at F-test p < 0.05), (3) significant positive correlation between human and rhesus macaque expression profiles (Pearson correlation p < 0.05 and r > 0) and, (4) significant time shift between human and rhesus macaque expression profiles (FPR = 11.7% for cortex and 11.5% for cerebellum). Taking the union of the genes that passed these criteria in the cortex or cerebellum yielded 1735 genes.
For these genes, we determined the direction of the time shift (acceleration, neoteny, or none) between human and rhesus macaque time series, in both cortex and cerebellum, by aligning rhesus macaque data to the human age scale. If at least 70% of time points showed consistent time shift direction, we classified this gene as accelerated or neotenic. Following this procedure, the vast majority of genes were classified as neotenic, with 1444 (83%) and 1469 (85%) in the cortex and cerebellum respectively. This result is consistent with the faster rate of macaque brain development and maturation that has been previously reported [24].
Overall, time shift measurements correlated positively between the prefrontal cortex and cerebellum (Additional file 1, Figure S8), with 560 out of 1750 genes tested showed strong positive correlation (Pearson p < 0.05, r > 0.5). The time shift profiles of these genes could be assigned to four clusters using the k-means clustering method (Figure 5A). Notably, the vast majority of these genes (452 genes in clusters 1 and 2) showed nearly identical time shift profiles in the prefrontal cortex and cerebellum. Does this similarity of time shift profiles reflect similarity in gene expression profiles in the two brain regions? Indeed, for some genes in both clusters 1 and 2, gene expression profiles in cortex and cerebellum followed the same trajectories (Figure 5B). Other genes, however, showed clearly distinct expression profiles in the two brain regions, while still sharing the same time shift profile. In a search for a biological explanation of this phenomenon, we tested cell type specificity and enrichment in biological functions specified by Gene Ontology (GO) annotation [25]. We found that genes grouped into Cluster1 based on their time shift profiles, and showing identical expression patterns in cortex and cerebellum (Cl1-group1), were enriched in white matter ([HT] p = 0.0002), and were annotated in GO terms: "lipid metabolic process" ([HT] p = 0.007) and "cellular lipid metabolic process" ([HT] p = 0.007). By contrast, Cluster1 genes showing different expression patterns in cortex and cerebellum (Cl1-group2) were enriched in gray matter ([HT] p = 0.0004) and annotated by GO as "nervous system development" ([HT] p = 0.05) (Figure 5C). Similarly, genes grouped into Cluster2 based on their time shift profiles, and showing identical expression patterns in cortex and cerebellum (Cl2-group2), were enriched in white matter ([HT] p = 0.02) and mature oligodendrocytes ([HT] p = 0.03), while genes showing different expression profiles (Cl2-group1) were enriched in gray matter ([HT] p = 0.006). Thus, our results revealed an interesting biological phenomenon: within one species, ontogenetic profiles are shared between the prefrontal cortex and cerebellum for genes expressed in white matter, but distinct for genes expressed in gray matter. Importantly, however, the time shift between the human and rhesus macaque ontogenetic profiles is perfectly synchronized for both white and gray matter genes. On an organismal level this observation might not be surprising. Changes in the rate of ontogenesis might be expected to operate on the brain as a whole, leading to synchronized delay in white and gray matter development in humans, compared to rhesus macaques. Our results confirm that, on the gene expression level, such synchronization can indeed be observed.