Training with genetically homogeneous and heterogeneous data

The broad range of transcriptomes from microarray and RNA-seq data that we use is detailed in Note A in S1 Appendix as are the methods used to prepare the data for use with TimeTeller. The data that is used to prepare TimeTeller’s probability model is referred to as the training data. The data that is then analysed using this probability model is called test data. In this paper we use four different training datasets and more details about these are in Note A in S1 Appendix.

Timecourse and intergene normalisation.

When combining training data from multiple tissues, for each gene in the REP we study the variation across the tissues in that gene’s expression time-series. This analysis (Fig B in S1 Appendix) shows that for RNA-seq data this variation is significantly greater than that found, for example, in the Affymetrix MoGene 1.0 ST and GeneChip Human Genome U133 Plus 2.0 microarray platforms that we have analysed. Therefore, for the RNA-seq training data, it is usually necessary if we are combining data from multiple tissues to carry out what we call timecourse normalisation.

Each of our training data sets is organised into time series for each gene in the REP with times t_k, k = 1, …, K, that are usually independent of the particular gene. We can normalise the data by replacing each of these time series by a normalised version which has mean expression zero and standard deviation 1. We call this approach timecourse normalisation. Following such normalisation of a training dataset, if we wish to test an independent test sample REV from a given tissue and gene we will have to normalise the REV using the offsets and scalings that were used in the timecourse normalisation of the training data for this tissue and gene. Such normalisation of test data is called timecourse-matched.

There is, however, a cost in using timecourse-matched normalisation because the test data from a particular tissue has to be normalised using the adjustments calculated for that tissue in the training data. This means that one can only use test data for tissues where we have a training time-series. Moreover, when using timecourse-matched normalisation on test data it is crucial that the training data are produced by the same transcriptomics platform.

Intergene normalisation avoids this. When timecourse normalisation is unnecessary or impossible because we do not have a training data set for the test data tissue, the data is normalised using intergene normalisation where, if g = (g_i) is a REV, the normalised levels are given by where μ and σ² are the mean and variance of the entries g_i. Essentially, this maps the REV onto its shape as a vector. It is also possible to usefully combine timecourse and intergene normalisation (see Table 1). Though the use of timecourse normalisation typically improves timing performance, we will show that intergene normalisation can also be remarkably effective (e.g., see Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Mean and median absolute timing errors for the training datasets.

Column 1 shows the normalisation used. Columns 2 and 3 show respectively the mean and median absolute timing error. Columns 5 and 6 show the mean and median absolute timing error after a correction is made using the timing displacement for the tissues or individuals as relevant. A-E. The apparent timing errors for the training datasets when a leave-one-out cross-validation approach was used. For the Zhang et al. microarray data [1] we compare using all the data (2h resolution) to only a subset giving 6h resolution. F. Timing results for Zhang et al. RNA-seq test data [1] when Zhang et al. microarray [1] is used as training data. G. As F. but with datasets swapped.

https://doi.org/10.1371/journal.pcbi.1011779.t001

We can also apply such timecourse normalisation to test data when this contains a time series; as several experimental model datasets do. However, any difference in amplitude between the training and test dataseries is then removed. On the other hand analysis using timecourse-matched normalisation for the test data maintains such a change in amplitude.

Similar considerations to the above apply when combining data across individuals instead of tissues as we do with the Bjarnason et al. human data ([30] and Note A in S1 Appendix). Timecourse normalisation can also be very useful when analysing microarray data and we have found it necessary when the training and test data come from different microarray platforms (e.g., as in Figs C and D in S1 Appendix). Table A in S1 Appendix. summarises the normalisations that were used for all the analyses shown in the Figs 1–5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. The color of data points etc (when not black) corresponds to the time when the data was sampled.

This coloring is used in a consistent way across all figures. A-C. Using local PCA projection to visualise data. The identity of each data point can be read from the legends to the right of each example. Only one projection for each example is shown but the differently timed projections have a similar quality. Examples showing all of the projections are in Figs E & G in S1 Appendix. A. The CT18 local PCs of the Zhang et al. microarray data [1] but using timecourse normalisation. B. A detail from a projection of the the Fang et al. test data [32] together with the Zhang et al. training data [1] as in A showing coherence of the Fang et al. WT data [32] and the gap between this and the KO data. Intergene normalisation is used. C. The Kinouchi et al. RNA-seq skeletal muscle test data for FED and FAST mice [34] plotted against the Zhang et al. RNA-seq training data [1]. Timecourse and timecourse matched normalisation is used. The ellipsoids shown are of the form (x − μ)^TΣ⁻¹(x − μ) = ε where μ is the mean of the estimated P(g|t) and Σ is its estimated covariance with ε chosen so that the ellipsoid should contain 97.3% of the training data (i.e. 3 standard deviations). This enables visualisation of the variation and covariation in the data.

https://doi.org/10.1371/journal.pcbi.1011779.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

A. The centred LRFs for the Zhang et al. mouse microarray data [1] using a leave-one-out analysis. They are centred in that the maximum of the curve is moved to noon. This makes the shapes of the curves clearer and more comparable. The black curve is the curve C(t|T) for T = 12 (Methods) that is used in the calculation of Θ. B,C. A Leave-one-out analysis of the Bjarnason et al. oral mucosa data [30]. B. Examples of the likelihood curves. C. Boxplots showing the apparent timing errors found for each individual ordered by their means. This shows the substantial timing displacements of some individuals. D. Showing how Θ is calculated using the LRF and the curve C(t|T). Θ is the proportion of time the LRF spends above C(t|T) i.e. the proportion of the time in the horizontal red curves. This is contributed to by the LRF around the highest peak, and by any secondary peaks or flat regions that go above C(t|T). E-H. Analysis of the Fang et al. Rev-erb-α KO data [32]. Uses l_thresh = −5. E. Plots of the Θ value against the estimated time T. The vertical lines show the true time with colours indicating the sampling time. WT timings are close to the true sample times and the KO times deviate from them. F and G. Boxplots of the maximum likelihood and Θ values showing significant differences between the WT and KO groups with p-values from the Wilcoxon rank sum test calculated using the Matlab ranksum function. Note that the smallest MLs are around e⁻⁴ which is why l_thresh was taken to be -5. Taking l_thresh = −4 gives entirely similar results. H. The centred LRFs for the WT and KO samples. I-L Analysis of the Kinouchi et al. FED/FAST skeletal muscle data [34]. This analysis used a logthresh of −12. The plots J-L are as for F-H but for the Kinouchi et al. data [34]. M-P Analysis of the Koronowski et al. data [39]. comparing WT, Arntl KO and Liver-RE data using l_thresh = −12. M The signed error boxplots show the timing dysfunction in the KO data as well as good recovery in the reconstituted Liver-RE clock but with a clear phase advance. N,O,P Boxplots of ML and Θ values, and centred LRFs for the three genotypes.

https://doi.org/10.1371/journal.pcbi.1011779.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

A-C. Analysis of the Boyle et al. data [43]. A. Boxplots of the Θ values for the smoker and never smoker individuals showing a statistically significant difference in the distributions. There is no statistically significant (Wilcoxon test) difference for the maximum likelihoods (Fig L(B) in S1 Appendix). B. The centred likelihood curves for the smokers and never smokers. C. The black, orange and light blue curves are estimates of the probability p_rand for random choices of the bad clock group of different sizes as in the legend. The red curve is for p_Θ(m). There were 5000 iterations for each curve shown. which gave a similar result to 10,000. The number of DE genes was decided using the BH adjustment method with p < 0.05 without any restriction on the minimum log fold change. The inset shows a blow up of these curves for m ≤ 10. From the blue curve (p_Θ(m) for 30 ≤ m ≤ 40) we see that for this range of m (unlike 8 ≤ m ≤ 20) it is very likely that only a very small number of DEGs are found. For the 66% of cases where a DEG is found there is a 99% chance that PER3 is among them and a 68% chance of NR1D2 being present. D-I. Analysis of the Feng et al. data [50]. D. Boxplots of the Θ values for the samples from individuals in the normal, cancer and dysplasia subgroups. These show a statistically significant (Wilcoxon test) difference in the distributions between the normal and cancer groups and the cancer and combined normal and dysplasia subgroups. E. Boxplots showing the predicted timing of the samples. F. Boxplots showing the predicted timing when all samples timed as before 7am and with a second peak are given the timing of the second peak. Of the 108 such samples 93 have moved. This suggests that the mistimed samples are primarily so because the wrong peak has a higher likelihood. G. Centred LRFs for the three subgroups. H. A study of differential effects between between those n individuals with the worse clocks according to the Θ stratification and those with better clocks. The black, orange and light blue curves are estimates of the probability p_rand as in C above but for the Feng et al. data [50]. The red curve is for p_Θ(m). I. Scatter plot of the projection of each REV in the Feng et al. data [50] with 12 < T < 16 (after using the second peaks if the first gives T < 7) against timing T. The red curve is a kernel smoothed estimate of the mean of . J. Distribution of the deviations in H. For each data point this is the horizontal difference between the data point and the red curve. A simple analysis shows that this is largely independent of and hence its standard deviation can be used as an upper bound for that of P(T|g).

https://doi.org/10.1371/journal.pcbi.1011779.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Examples of PCP plots for the Bjarnason et al. data [30].

These (full set in Fig M in S1 Appendix) show the strong linear relationship between the REP gene phases and the timing deviations in this data. Each point corresponds to an individual. The regression was carried out using Matlab’s fit function and Cosinor [52] was used to estimate the gene phases from the time series of each individual. The p-values test the hypothesis that the slope of the line is non-zero and are given by the F-test using the Matlab functions coefTest and fitlm.

https://doi.org/10.1371/journal.pcbi.1011779.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5.

A-E. Analysis of the Acosta-Rodríguez et al. data [37]. This data is analysed as test data using timecourse-matched normalisation and timecourse normalised Zhang et al. RNA-seq data [1] for training. We use l_thresh = −8. The timing, Θ and ML for each sample point shown is a suitably averaged value for the two replicates with the same feeding condition and age. A. Box plots of the apparent timing error for the different conditions and ages. The mean value of each box plot gives the timing displacement for the condition and age. Only 12 of the 66 comparisions are not significant at the p = 0.05 level using the Wilcoxon test. B. Box plots of the ML for the different conditions and ages. C. Box plots of the Θ value for each condition and age. See Table D in S1 Appendix for statistical analysis of the differences. D. Centred LRFs for each condition and age. E. Legend. F. PDP plots for the genes in the REP, some cell cycle genes and some of the genes highlighted in [37]. See Fig N in S1 Appendix for more information. The gene phases were measured by Cosinor [52]. G-K. Analysis of the Mure et al. data [2]. In each plot the color corresponds to the tissue as shown in I. The data from the central tissues is used for training. G. TimeTeller predicted time T vs the sample time using leave-one-out analysis for each sample from the 33 tissues. H. Box plots of the signed apparent errors for the samples for each of the 33 tissues in order of increasing timing displacement. Each of the 7 leftmost boxplots is significantly different from each of the 7 rightmost boxplots at the p = 0.005 level. I. Legend for G-K. J. Timing displacements for each tissue. K. Some examples showing PDP plots of the gene phases and the TimeTeller timing displacement for all 33 tissues from Mure et al. [2]. The p-value and r² values for the other genes are in Table C in S1 Appendix. Note that the r²s (Fig 4) for core clock genes are much smaller than those for the Bjarnason et al. human data [30] and the Acosta-Rodríguez et al. data [37] in F.

https://doi.org/10.1371/journal.pcbi.1011779.g005

Multidimensional visualisation provides important information about phenotype

When constructing the probability model, the TimeTeller algorithm projects the G-dimension REVs into fewer dimensions using a local version of principal component analysis (Methods, and Note E and Figs E and F in S1 Appendix). This gives a different projection for each time in the dataset and the algorithm extends this to all times around the day. If for the G-dimensional data the distributions P(g|t) are approximately multivariate normal (MVN) then the corresponding distributions of the projected data optimise the capture of the dominant gene-gene correlations after projection (see Methods). We find that for our datasets d = 3 is sufficient for this (e.g., see Fig F in S1 Appendix) and the resulting 3-dimensional model of the clock provides a very informative visualisation.

Fig 1A shows such a visualisation for the mouse multi-organ microarray training data from Zhang et al. [1] when timecourse normalisation has been applied. TimeTeller actually produces such a local projection visualisation for each time in the training dataset as shown in Figs E and G in S1 Appendix but normally inspection of just one of these is adequate and we only show one in Fig 1. With each such visualisation we also show the curve given by the means of the estimated distributions P(g|t) as t varies over the day. Also in such plots we often provide for a sample of times t an ellipsoid showing the covariance structure of the estimated distribution P(g|t) (see caption of Fig 1). We color the training data points and mean curve by time with a color coding as given in the legend of Fig 1 using the sample time for the data points. The same color coding is used throughout the paper.

Fig 1B plots microarray test data from Fang et al. [32] comparing it with the Zhang et al. microarray training data [1] in Fig 1A. This test data compares liver samples of Nr1d1 (Rev-erb α) knock-out (KO) and wild-type (WT) mice entrained to light-dark (LD)12:12 cycles. The gene Nr1d1 is a core clock gene of the mammalian circadian clock important in one of the interlocked feedback loops and a key link to metabolism [33]. Knocking it out leaves a functional but perturbed clock when compared to WT mice [32]. Since Nr1d1 is a member of the default REP it would not be surprising that TimeTeller could distinguish Nr1d1 KO mice from WT mice, and indeed this is the case. Therefore, for this validation, we exclude Nr1d1 from the REP genes. The visualisation shows that while the WT data appears to fit well with the training data, the KO data has a consistent substantial difference. TimeTeller is able to detect this apparent difference in each of the four KO samples and shows a coherent difference from WT. This suggests that the Nr1d1 KO mice have a significantly perturbed clock when compared to WT mice (Fig 1B). However, it is still somewhat functional as it gives approximately correct timing and the level of sample variation between WT and KO is similar. We investigate this further below.

The other test data we visualise (Fig 1C) in this figure is from Kinouchi et al. [34]. This contains samples analysed by RNA-seq from mouse skeletal muscle taken around the clock in LD 12:12 [34], and compares mice that had been fed ad libitum (FED) with mice that had been starved for exactly 24 hrs prior to point of sampling (FAST). On the one hand, while FED samples align with the RNA-seq training data, the FAST samples are substantially perturbed (Fig 1E). On the other hand, the FAST samples show consistency in that for a given sample time they tend to cluster together. A similar visualisation for the liver samples from [34] is given in Fig H in S1 Appendix. It should be noted that the test samples from Kinouchi et al. [34] have been collected in LD whereas the training dataset was collected on the first three days in constant conditions. Interestingly, there is little difference between FED (control) and training dataset mice, which might be due to the fact that the free-running period of these WT C57Bl/6 mice is around 23.8 hours.

Other examples demonstrating the utility of such visualisation are discussed below.

Analysis of single test samples

TimeTeller’s estimate of P(t|g) from the training data is used to analyse test data. For a normalised test data REV g our estimate L_g(t) of P(t|g), which we regard as a function of t, is referred to as the likelihood curve (LC) for the corresponding transcriptomics sample. The quantities for functionality assessment are associated with this LC. For example, we define the internal phase T of a test REV g as the time at which the estimated likelihood function L_g(t) ≈ P(t|g) is maximal i.e., the maximum likelihood estimate. Given T, we define the likelihood ratio function (LRF) as R_g(t) = L_g(t)/L_g(T), i.e. it is the LC but normalised so that the value at the maximum is 1. The internal phase can be compared with the SCT but, as noted above, there may be consistent phenotypic deviations of T from the SCT in genetically heterogeneous populations.

It is important to emphasise that when we analyse test data the results for any test data sample are independent of the results for any other test data sample. This is because the calculation of the likelihood curve of a test data sample only involves the probability model and the test data sample and has nothing to do with the other test samples. Therefore, the result for any test sample will be exactly the same as if it were the only sample in the test dataset.

Fig 2A shows the estimated LRFs for the Zhang et al. microarray data [1]. Each LRF’s highest peak is centered at 12noon to enable visual comparison of many LRFs, a plotting technique used throughout the figures. Many examples of estimated LCs and LRFs can be seen in Figs 1–5 and the S1 Appendix. The resulting predicted timing plotted against the sample time is shown in Fig 2A together with the times corrected to allow for the chronotype explained in the section below on timing. LCs for the Bjarnason et al. human training data [30] are shown in Fig 2C and Fig I in S1 Appendix.

These show the general form of the LCs and demonstrate that one can clearly observe qualitative differences between one individual’s LC and those of the others.

Multiple tools for assessing functionality in test data

Importantly, this consistency of good apparent timing errors, high ML and low Θ was also observed in the analysis of the various WT/control components of the test data sets considered. For example, using the Zhang et al. microarray data [1] for training and intergene normalisation our analysis of the microarray timecourse control dataset created by LeMartelot et al. ([36] and Note A in S1 Appendix) produced a mean absolute error for time estimation of less than one hour and Θ values similar to those found in the training data. Similar results were found for the Acosta-Rodríguez et al. data [37] for ad libitum fed mice using the Zhang et al. RNA-seq data [1] for training and timecourse-matched normalisation, and for liver microarray test data from Hughes et al. [38] after training with the Zhang et al. microarray data [1] (Fig C in S1 Appendix). For the latter we used timecourse normalisation for the training and test data as the microarray platforms are different, demonstrating good results across different platforms.

To further test the use of TimeTeller across different transcriptomics platforms we carried out a cross-validation experiment where we trained TimeTeller on the Zhang et al. microarray data [1] and used this to test the Zhang et al. RNA-seq data [1] and vice-versa (Fig 1A and Fig D in S1 Appendix). This not only tests the robustness of our approach but also examines the effectiveness of timecourse normalisation in allowing us to work across different transcriptomics technologies. The timing results are given in lines F and G in Table 1 with small mean and median errors of a size compatible with the within-dataset leave-one-out analysis. As well as the relatively small timing errors we observe informative visualisation and consistent Θ values.

Across the various datasets we consider, this good timing, high ML and low Θ for WT/control test data almost always differed from that found for the perturbed test data. For example, in Fig 2H–2O we consider the timing, ML and Θ diagnostics for the Fang et al. [32] and Kinouchi et al. test data [34] discussed above and use them to illustrate how to gain more insight into dysfunction. In such an analysis one should start with an assessment of the MLs and an inspection of the LRFs for training, control and test data.

From this one can choose an initial value for the important parameter l_thresh using the approach described in Methods and Note C in S1 Appendix. This parameter truncates the likelihood curves so they do not go below exp(l_thresh) and this plays an important role in ensuring that the incorporation of the local likelihoods into a global one (see Methods) is not wrecked by inaccurate and uninformative exceptionally low local likelihoods.

For the Fang et al. data [32] and the Kinouchi et al. skeletal muscle test data [34] we see that the ML values for the perturbed test data are significantly lower than those for the control test data (Fig 2J and 2N) suggesting that dysfunction of the lowML type is present in the perturbed systems. For the test data from Fang et al. [32] we also observe significant differences between the WT and Nr1d1KO samples for the timing and Θ diagnostics (Fig 2H–2K). The timing T of the KO observations is significantly further from the true timing (Fig 2H). Analysis as in Methods and Note C in S1 Appendix suggests setting l_thresh around -5 but the results are very similar for any value between -4 and -7. However, the centred LRCs also indicate that there is a significant amount of highTvar dysfunction in the KO sample because of the second peaks and increased width of the LRCs near the maximum.

For the Kinouchi et al. skeletal muscle data [34] the FED samples show uniformly small errors in timing T (MAE 0.44h) and uniformly low Θ values (Fig 2L and 2M). In contrast, the timings T of the FAST samples are clustered around ZT 18–24 reflecting the tight clustering seen in the visualisation (Fig 1E) and the MAE is significantly greater at 4.04h. So far as dysfunction is concerned, the situation is somewhat different from the Fang et al. data [32] since, although there are also significant differences in ML and Θ values between FED and FAST (Fig 2I–2L), there are no second peaks in the centred LRFs contributing to Θ (Fig 2L). Consequently, the primary difference between the FED and FAST samples is due to the significant difference in the MLs. Thus, only substantial lowML type dysfunction is present. The stratification by Θ in Fig 2K reflects this.

The Kinouchi et al. data [34] provides a very informative example of how the choice of the parameter l_thresh works because the maximum likelihoods ML for both the FED and FAST liver data are significantly higher than that for the skeletal muscle data discussed above and this means that different values of l_thresh are appropriate. The discussion in Note C and Fig H in S1 Appendix shows that the value for the liver data should be substantially larger at -6 or -7 rather than -12. When this value is chosen the results for the Kinouchi et al. liver data [34] are similar to those above for the skeletal muscle data (Fig 2I–2L).

In order to understand the effect of fasting on the amplitude of core clock components Kinouchi et al. [34] needed to treat the data as though the FAST samples belonged to a continuous time series even though each timepoint was proceeded by 24 hours of starvation. This underlines a significant extra advantage of TimeTeller because the FAST test samples can be considered independently from one another.

Koronowski et al. [39] compared the liver transcriptomes of wild-type (WT) with whole body Arntl deficient mice (KO) or Arntl KO mice with liver-specific Arntl reconstitution (Liver-RE). Their data enables us to test TimeTeller’s sensitivity to not only the substantial KO perturbation but also the much subtler one of the Liver-RE. Our analysis shows statistically significant differences in timing between WT, KO and Liver-RE including the phase advancement noted in [39] of the Liver-RE clock relative to WT (Fig 2M–2P). The MLs (resp. Θs) for the KO data are significantly smaller (resp. larger) than for both the WT and Liver-RE data with no significant difference between WT and Liver-RE (Fig 2N and 2O). However, the LRFs (Fig 2P) suggest a clear difference between WT and Liver-RE data in that, unlike WT, about half of the Liver-RE samples have a significant extra peak suggesting a contribution of highTvar type disruption and a hypothesis that this is causing the observed timing change in the Liver-RE data. Significant second peaks are also observed in about a half of the the KO data suggesting a combination of some highTvar dysfunction combined with the significant lowML dysfunction.

The potential for the use of the Θ stratification to identify differential effects in patients

It is particularly interesting and important to apply TimeTeller to genetically heterogeneous human data because it allows us to test the idea that it can uncover corresponding heterogeneity in the “clock” phenotype or effects on individuals such as patients.

We firstly consider data from a study of the effects of cigarette smoke on the human oral mucosal transcriptome, In this study (Boyle et al. [43]) transcriptomes from buccal biopsies of 39 current smokers (≥ 15 pack-year exposure) and 40 age- and sex-matched never smokers (< 100 cigarettes per lifetime) were analysed and compared. The authors found that smoking altered the expression of numerous genes but none of those found were core clock genes nor did they consider the effect of smoking on the circadian clock. They found smokers had increased expression of genes involved in xenobiotic metabolism, oxidant stress, eicosanoid synthesis, nicotine signalling and cell adhesion and decreases were observed in the genes CCL18, SOX9, IGF2BP3 and LEPR. It has been reported elsewhere that smoking has an impact on multiple sleep parameters and significantly lowers sleep quality [44–46] and this was confirmed in an experimental study which also correlates poor sleep to inflammation [47] while inflammation has been linked to clock disruption. Moreover, CS exposure has been shown to cause circadian disruption in the lungs of WT mice and this is exaggerated in the Nr1d1 knockouts [48] and has a connection to Arntl [49].

Interestingly, when analysed by TimeTeller (Fig 3A–3C) we see a clear and statistically significant difference between the Θ values of the never smoked and smoking individuals (Fig 3A) which is reflected in the 3D visualisation (Fig L(A) in S1 Appendix). Inspection of the LRFs show that the variations in Θ come mainly from second peaks rather than low ML (Fig 3B). Indeed, the ML values for smokers and non-smokers were not significantly different although the smokers had more observations with a very small ML (Fig L(A-C) in S1 Appendix). The lowest values were around e⁻¹¹ suggesting that a l_thresh of about -12 would be appropriate.

A significant proportion of the smokers had Θ values similar to those of the never-smokers but many had much higher values (Fig 3A). Therefore, we asked if we could identify differentially expressed genes (DEGs) between the individuals with high Θ versus those with lower Θ. To do this we tested for differential gene expression between the n worst clocks (defined as the bad clock group (BCG)) and the others (good clock group (GCG)) adjusting the p-value appropriately to allow for the multiple testing. For a fixed l_thresh in the range from 11 to 13 with n between 8 and 20 we found many differentially expressed genes (DEGs) at the appropriately adjusted p = 0.05 level including some clock genes (Fig L(D,E) in S1 Appendix). However, the particular genes found were sensitive to changing the value of l_thresh among the suggested values of -11, -12 or -13 and changing the group size n.

We calculated that the probability of finding such numbers of DEGs by chance is extremely low (Fig L(H) in S1 Appendix) and we noticed significant differences between the behaviour when the BCG size was in the range 8 to 20 from that when it was 30 to 40. Therefore, we estimated by simulation the probability p_rand(m) of finding m or more DEGs by chance when we choose a random group of n individuals for our BCG and compared this to the probability p_Θ(m) of finding m or more DEGs when the stratification by Θ is used to choose the BCG and l_thresh and n are chosen randomly in the ranges -11 to -13 and 8 to 20. We find that uniformly in m, p_Θ(m)/p_rand(m) > 100 (Fig 3C). We get an interestingly different result if we instead let the group size n range between 30 and 40. The probability p_rand(m) behaves in approximately the same way but p_Θ(m) does not (Fig 3C(inset)). For m very small p_Θ(m) is high but as m increases p_Θ(m) rapidly decreases to values much smaller than those for p_rand(m). There is a 34.26% chance of getting no DEGs but when this is not the case there is a more than 99% chance of getting the gene PER3 and a 68% chance of getting NR1D2. Thus, this analysis identifies two interesting groups of individuals with a nontrivial transcriptional phenotype that distinguishes them from the individuals with good clocks. One of these groups appears to be associated with differential expression of PER3 and NR1D2, genes not identified in the original paper where all non-smokers and smokers were compared.

In conclusion, any link between smoking and clock dysfunction is likely to be complex, but these results suggest that in a genetically heterogeneous population where the effects of a perturbation such as smoking are likely to be diverse, TimeTeller’s Θ stratification can help identify individuals or groups where the smoking effect is significant.

As a final example of this section we consider the distribution of Θ values by disease state for the transcriptomic data of healthy or dysplastic oral mucosa and oral squamous cell carcinoma (OSCC) from Feng et al. [50]. Since the lowest ML values were around e⁻¹¹ a l_thresh of -12 was used. The ML values for normals and cancer were not significantly different although the cancer group had more observations with a very small ML (Fig L(F) in S1 Appendix). However, there is a highly significant difference in median Θ values between the cancer group (167 individuals) and the the normal mucosa group (45 individuals) (p < 0.002) (Fig 3D). Moreover, there appears to be significant dysfunction in terms of timing estimation (Fig 3E) that can be significantly ameliorated if the second peaks in the LRF is used for timing when the first peak is clearly misleading (Fig 3F, details below). Inspection of the LRCs (Fig 3G) shows that, as for the Boyle et al. data [43], the variations in Θ come mainly from second peaks rather than low ML.

As for the Boyle et al. data [43] we asked if there are DEGs between the worst clocks in the cancer group (high Θ) and the best clocks within the same group and carried out a similar analysis. For genes in general and BCG sizes n between 12 and 40 we find similar results with p_Θ(m)/p_rand(m) > 100 for the number m of DEGs between 2 and 1200 (Fig 3H). Many of these DEGs are associated with gene signatures such as DNA repair, E2F targets, G2M checkpoint and the mitotic spindle. However, we do not find any groups like that for the Boyle et al. data [43] (with n between 30 and 40) that have very low numbers of specific DEGs.

A study of the estimated timing T for this data (Fig 3E) was very informative. The estimates for the normal data are generally between 7 am and 3 pm. A large number of cancer samples have unlikely times well outside the normal working day and the median is clearly much too early. Interestingly, it appears that the mistimed samples are primarily so because the likelihood curve has a second peak (Fig 3G) and the peak giving an unreasonable timing estimate is slightly higher than one giving the best estimate. In fact, there are 109 samples whose timing T is before 7am and 93 of these have a second peak. if we replace the timing by that given by the second highest peak, the great majority moved to a time firmly in the early afternoon between 12noon and 4pm (Fig 3F). As a result 74% of all samples then fall in this time slot and only 7% remain before 7am. The analysis in the next section indicates that this corrected timing is likely to be the correct time of sampling to within approximately 0.4h.

TimeTeller’s precision on non time-stamped cancer data

The only method currently utilised to estimate the precision of timing/phase algorithms is to use time-stamped data and compare the algorithm’s predicted times T with the SCT time stamps t. However, such a measure of precision is problematic when the individuals, tissues or conditions have a nontrivial molecular chronotype as is the case with the human data considered here and cannot be done if the data is not time-stamped. A related test which avoids these problems is instead to determine the variance or standard deviation of the distribution P(T|g) where T is the predicted time and g is the relevant REV (Note F in S1 Appendix). This addresses the question of how well the estimated timing T is determined by the REV g. Interestingly, we can calculate this precision measure even in some cases where we have no timing data and where there is dysfunction and the Feng et al. data [50] gives a very informative example of this.

To illustrate this we study the 77% (176 samples) of that data for which the estimated timing T after adjustment by second peaks is between 12 noon and 4pm (Fig 3F). We ask if within this data we can see coherent timing structure or not. We can estimate the required standard deviation by carrying out a principal component (PC) analysis of the expression data (see Note F in S1 Appendix) and plotting the projection of these data onto the first PC against the predicted time T (Fig 3I). An upper bound for the standard deviation of P(T|g) can be estimated from this (Note F in S1 Appendix) and we obtain an estimate of less than 0.4 hours. If we consider all the deviations from the mean for the timings T (given by the horizontal deviation of the relevant data point from the red curve in Fig 3I) across all of the REVs in this data we obtain the distribution shown in Fig 3J. Remarkably, although the data is not timestamped and has significant dysfunction giving rise to significant second peaks, TimeTeller is able to accurately measure the internal phase T of the clock as a function of the REV g.

We carried out a similar analysis (see Note F in S1 Appendix) and found similar results but a bigger standard deviation of 0.83h for the large breast cancer dataset analysed by Cadenas et al. in [51]. In this case there is no need for adjustment for second peaks as 86% of the data has its predicted time T between 10am and 8pm (Note F in S1 Appendix). We believe this approach gives a new simple method to assess timing performance.

Comparing clocks across individuals, conditions and tissues

Current analyses comparing the circadian clock across individuals, tissues and conditions such as the three studies we consider below proceed by analysing the behaviour of the individual interesting genes separately. Such analyses tend to focus on the level of expression and do not take into account correlations between related genes. We asked whether using TimeTeller such an analysis could be done in a more integrated way treating the clock as a noisy multigene dynamical system (and hence using correlations) and whether such an approach uncovers some aspects that are hard to see when done gene by gene. The key results here are that it enables us to identify coherent differences in timing across individuals, conditions and tissues and that using these we can determine in a quantifiable way if the timing differences come from a more or less coordinated change in gene phases.

Probing the effect of changes in the core clock on downstream genes

Changes in the core clock will affect the regulation of rhythmic genes that are downstream of it. Current methods allow one to check whether these genes remain rhythmic when the clock is perturbed in some way but TimeTeller also allows examination of the extent to which they maintain their relationship with the clock in a coherent fashion. The way in which the different conditions of the Acosta-Rodríguez et al. mouse data [37] changed the phase of the core clock provides a very interesting example where we can demonstrate such an analysis.

Firstly, we noted that for the genes in the REP, all clock genes displayed approximately linear phase changes while for the other genes (Hlf, Wee1 and Cys1) this was not the case for Cys1. We then used this analysis to look at the effect of the clock phase changes upon some other genes that are rhythmic in the liver of AL fed mice. In particular, we inspected the plots for some cell cycle genes and also a number of the genes identified in Acosta-Rodríguez et al. [37] as affected by the CR conditions or ageing. We find that inspection of the PDP plot for these genes gives clear and significant insight into the level of this coherence which we quantify by the p-value and r² of the PDP plot.

Of the cell cycle genes Wee1, p21, P53, Timeless, CyclinA, CHK2, CyclinB1, CyclinE2 and ATM, it appears that only Wee1, p21 and CyclinE2 are rhythmic in the liver in the AL conditions. These three genes maintain coherence with the clock under the other conditions with Wee1 doing so strongly (r² = 0.96) followed closely by p21 (r² = 0.92). The coherence of CyclinE2 seemed somewhat weaker (r² = 0.59). All of the other genes had r² < 0.4 and appeared incoherent. There is a very strong correlation between the maintenance or absence of coherence and rhythmicity or non-rhythmicity.

In Acosta-Rodríguez et al. [37] a number of genes that were affected by ageing or the CR conditions were highlighted and sorted these into four categories: those susceptible to ageing-related changes under any condition tested, those related to fasting conditions, timing related genes and genes associated with effects on circadian cycling such as rhythmic damping. Our analysis using PDP plots for these genes clearly identifies which of them move coherently with the core clock under the different feeding conditions. None of the timing related genes stayed coherent and, amongst the fasting genes, only Hal1 (r² = 0.66) was. Several ageing genes show some level of coherence (Fig N in S1 Appendix) Serpine1 (r² = 0.66) Adora1 (r² = 0.68) Got1 (r² = 0.65) Lepr (r² = 0.68) Pfkfb5 (r² = 0.88). For the genes affecting circadian cycling. while Gys1 (r² = 0.15) and Per1 were incoherent, the rest were coherent: Arntl (r² = 0.99), Nr1d1 (r² = 0.70), Per1 (r² = 0.69), Per2 (r² = 0.97) and Pck1 (r² = 0.60). For a significant number of the genes affected by ageing or the CR conditions, while the gene is not coherent under all conditions it is coherent under a significant number of the conditions with the less extreme timing deviations. This can be seen from the PCPs and was the case for Per1 which seems coherent under all conditions except the four CR-day conditions.

These results demonstrate that such an analysis can give a novel overview of gene response and whether a given gene maintains coherence with the clock when the clock timing changes. Such coherence is associated with genes that show good linearity with a significant slope in the PDP plots. Consequently, TimeTeller can be used to investigate function and dysfunction in genes controlled by the circadian clock when the clock is perturbed.

Ethics statement

The study associated with the Bjarnason et al. human data [30] was approved by the Sunnybrook Health Sciences Centre Research Ethics Board. Project identification number 396–2004. Written informed consent was obtained from each subject as requested by the Research Ethics Board.

Probability model constructed from training data

The training data will have been collected at sample times t_i, i = 1, …, N_t. In the training data used here the number N_s of samples at each time point is the same. Therefore, if the samples are indexed by j, the G-dimensional REVs with sample time t_i can be labelled by i and j and denoted .

In three of the training datasets the instances j correspond to different tissues (with replicates in one case) and in the other (Bjarnason et al. [30]) to different individuals. Each g_ij is then normalised using timecourse and/or intergene normalisation as described above resulting in vectors that will be used to train TimeTeller. The issue of batch effects is considered in Note A in S1 Appendix.

To construct the probability model we firstly construct one for each timepoint t_i in the training data by using the local statistical structure of the data at that timepoint and then we combine these. Associated with this time t_i is the set of N_sG-dimensional vectors , j = 1, …, N_s. We calculate the principal components U_i,k of this dataset and then use the first d of these to define a projection P_i of the normalised training data into (Note F in S1 Appendix) i.e. where U_d is the matrix made up from the column vectors U_i,k for k = 1, …, d. We then fit a multivariate normal distribution (MVN) to the points . The dimensionality d is chosen so that there are enough vectors in to fit a d-dimensional multivariate Gaussian (using the MATLAB function fitgmdist) while ensuring that most of the variance in the data is captured by the d-dimensional projection (e.g. see Fig F in S1 Appendix). In our case we take d = 3.

Now we fix a time t_i and consider the means μ_j and covariance matrices Σ_j of the distributions . We fit a periodic piecewise cubic hermite interpolating polynomial spline through the μ_j and each of the d(d + 1)/2 entries that determine Σ_j so as to extend μ_j and Σ_j to all times t between the time points checking that the Σ_j are positive definite and moving them to the nearest positive definite matrix if this is not the case. We thus obtain μ_i(t) and Σ_i(t) and thus the associated family of d-dimensional MVN distributions for all times t between the first and last data times. For these splines we use the MATLAB function perpchip as this respects the periodicity in t. Our implementation offers some alternatives to perpchip but these are not used here. This family of MVN distributions indexed by time is what we refer to as the probability model.

The likelihood curve L_g(t) and the log threshold l_thresh

Now we define the likelihood curve L_g(t) where g is a REV from either training or test data. Having calculated the probability model, for a given REV g, for each of the time indices i we define the likelihood curve associated with the ith timepoint using the probability given by the MVNs i.e. where g^norm is the vector obtained after normalising g with the relevant normalisation.

The idea is to obtain log L_g(t) by averaging these individual log likelihoods log L_g,i, i = 1, …, N_t but some modification is needed. We will need to fix a lower threshold l_thresh < 0 and replace each log L_g,i by max{log L_g,i, l_thresh} in the sum so that L_g(t) is defined by .

This truncation is necessary to ensure that this sum is not wrecked by inaccurate exceptionally low values of one L_g,i affecting robust high values of another at the same t. A curve L_g,i may take on very low values away from its maximum and the exact values of these very low probabilities may well be unreliable and inaccurate. If this happens at a t value for which another such curve L_g,j has a high accurate value then this may badly affect the estimate of L_g(t). The way to choose the value of l_thresh is discussed in Note C in S1 Appendix.

Definition of Θ

The clock dysfunction metric Θ is defined to be the proportion of time t where the LRF is greater than C(t|T) = η(1 + ϵ + cos 2π(t − T)/24) which is a scaled cosine function phase shifted so that the maximum is at T (Fig 2D). The parameters must satisfy 0 < ηϵ < η(2 + ϵ) < 1. Although we have experimented with changes, the effect on Θ of changing η and ε is clear (see Fig 2D) and we have seen no reason for changing them from the values we have used here.

Choice of parameter l_thresh

The key considerations underlying the choice of l_thresh are that it should be as large as possible subject to the conditions that (i) very few training and control samples have flat regions that significantly intersect C(t|T) so that they contribute significantly to Θ, and (ii) as many as possible of the test data samples should have MLs above exp(l_thresh).

There are two reasons we do not want l_thresh to be decreased further than necessary. Firstly, the considerations above about the need to protect against inaccurate exceptionally low values of some L_g,i and, secondly, because if l_thresh is reduced too far structure in the LRF at times that are away from the time T is likely to be removed. This happens because, if the L_g,i have their maxima not too far from T then decreasing l_thresh causes a much bigger decrease in the likelihood L_g(t) for t away from T than near to T and therefore decreases the LRF away from T while maintaining the peak structure near T. If the dysfunction is mainly manifested by low ML then decreasing l_thresh by too much moves all the flat regions in the LRFs down below C(t|T) while if it manifested by structures such as second peaks then these are also decreased below C(t|T). In both cases this results in a decrease in Θ. These phenomena are illustrated in Note C and Fig H in S1 Appendix.

If the maximum value log ML of log L_g(t) is only just above l_thresh, then it and the corresponding likelihood ratio curve will have intervals on which they are flat. If this is the case then the length of these flat intervals above the minimum of the curve C(t|T) can contribute to Θ. This contribution has interesting information in it because it is related to how low the maximum value ML of L_g(t) is.

If the criterion (ii) results in too small a value so that too much structure has been removed, it is then generally acceptable to set l_thresh at a higher value provided that the number of training samples violating (i) does not get too large. It so also desirable that the the number of test samples violating (ii) is not too large as otherwise many samples have Θ = 1 meaning that these samples do not have a non-trivial stratification even though they will be distinguished as having higher dysfunction than other samples.