Diagnosis and treatment of circadian rhythm sleep-wake disorders both require assessment of circadian phase of the brain’s circadian pacemaker. The gold-standard univariate method is based on collection of a 24-hr time series of plasma melatonin, a suprachiasmatic nucleus-driven pineal hormone. We developed and validated a multivariate whole-blood mRNA-based predictor of melatonin phase which requires few samples. Transcriptome data were collected under normal, sleep-deprivation and abnormal sleep-timing conditions to assess robustness of the predictor. Partial least square regression (PLSR), applied to the transcriptome, identified a set of 100 biomarkers primarily related to glucocorticoid signaling and immune function. Validation showed that PLSR-based predictors outperform published blood-derived circadian phase predictors. When given one sample as input, the R2 of predicted vs observed phase was 0.74, whereas for two samples taken 12 hr apart, R2 was 0.90. This blood transcriptome-based model enables assessment of circadian phase from a few samples.https://doi.org/10.7554/eLife.20214.001
Many physiological and molecular variables exhibit 24 hr rhythmicity generated by central and peripheral circadian oscillators, as well as environmental and behavioral influences (Morris et al., 2012; Mohawk et al., 2012). The master circadian pacemaker located in the suprachiasmatic nucleus (SCN) of the hypothalamus drives 24 hr rhythmicity in endocrine and neural factors such as cortisol, melatonin and the autonomic nervous system (Brancaccio et al., 2014; Kalsbeek and Fliers, 2013). Under normal conditions, the SCN is synchronized to the 24 hr light-dark (L-D) cycle by a light-input pathway (Lucas et al., 2014). The effects of light on the SCN depend on the timing, that is, the circadian phase of the SCN, at which light is administered (St Hilaire et al., 2012). The SCN is also sensitive to feedback from the pineal melatonin rhythm which impinges on the SCN through its melatonin receptors (Pevet and Challet, 2011). Exogenous melatonin can phase-shift and entrain the SCN, but its effects depend on the circadian phase of melatonin administration (Lewy et al., 2005).
Synchronization of the SCN to the L-D cycle and synchronization of the sleep-wake cycle to the SCN are disrupted during shift work and jet lag and in circadian rhythm sleep-wake disorders (CRSWDs). Examples of CRSWDs are delayed sleep-wake phase disorder, which is prevalent in adolescents (Micic et al., 2016), and non-24-hr sleep-wake disorder, such as occurs in blind individuals (Sack et al., 1992; Sack and Lewy, 2001). Treatment of these circadian disorders by melatonin or light exploits the circadian phase-dependent response of the SCN to these stimuli (Auger et al., 2015). Optimal treatment requires that the circadian phase of the SCN is known and is thereby one example of chronotherapy, an increasingly important field in drug development and medicine (Dallmann et al., 2016). Knowledge of the circadian phase of the SCN is also informative for the diagnosis of CRSWDs, assessment of adaptation to new time zones, adaptation to shift work (Mullington et al., 2016) and identification of circadian rhythm disturbances in psychiatric conditions. Conventionally, circadian phase of the SCN is assessed through long-term collection of biomarkers in the field or in the clinic. Examples are repeated urinary 6-sulfatoxy melatonin assessments (Flynn-Evans et al., 2014; Lockley et al., 2015), long-term core body temperature recordings (Klein et al., 1993), or frequent blood or saliva sampling in dim light, for the assessment of melatonin rhythm (Klerman et al., 2002; Danilenko et al., 2014). The plasma rhythm of melatonin is driven by rhythmic melatonin synthesis in the pineal and under control of the SCN (Arendt, 2005). In human clinical research, melatonin is therefore considered a gold-standard proxy for SCN phase. However, circadian phase assessments based on the plasma melatonin rhythm are time consuming, burdensome and costly. This would be greatly reduced if circadian phase could be assessed from only one or few blood samples.
At the molecular level, circadian oscillators consist of a network of genes that, through a transcriptional-translational feedback loop, generates tissue-specific circadian rhythmicity in transcription and translation (Partch et al., 2014; Koike et al., 2012). In human blood, 6.4–8.8% of the transcriptome exhibits circadian rhythmicity when sleep occurs in phase with the melatonin cycle and in the absence of a sleep-wake cycle (Archer et al., 2014; Möller-Levet et al., 2013). This peripheral circadian rhythmicity is, however, markedly disrupted when sleep occurs out of phase with the melatonin rhythm. Under these conditions, only 1% of transcripts are rhythmic (Archer et al., 2014). The human blood transcriptome is also affected by both chronic sleep restriction and acute total sleep deprivation (Möller-Levet et al., 2013). This sensitivity of the peripheral circadian transcriptome to factors such as sleep timing and sleep restriction poses a significant challenge for the assessment of circadian phase from the blood transcriptome. In many target conditions and patient populations, for example, shift work and non-24-hr sleep-wake disorder, sleep timing and duration are altered. Nevertheless, one or two blood transcriptome mRNA abundance profiles could potentially contain information about the circadian phase at which these samples were taken under all these conditions.
Phase prediction methods based on omics data from a few samples of blood or other tissues have been developed. Some methods selected the transcripts (features) to be used as predictors based on existing knowledge of their rhythmic expression and/or involvement in circadian processes (Lech et al., 2016), that is, a transcriptome wide search was not performed. Other methods required predictors to be rhythmic. These methods first fitted a rhythmic function (e.g. cosine-wave [Ueda et al., 2004; Minami et al., 2009] or rhythmic spline [Hughey et al., 2016]) to time-series of transcriptomic or metabolomic data to identify (24 hr) rhythmically expressed features, which were then used as predictors. None of these methods, either applied to the transcriptomic profiles of different mouse organs (Ueda et al., 2004; Minami et al., 2009), or the metabolomic (Minami et al., 2009) or transcriptomic profiles (Lech et al., 2016) of human blood, specifically aimed to predict the phase of the melatonin rhythm. Moreover, the performance of these methods was not assessed in conditions of desynchrony of the sleep-wake cycle and the melatonin rhythm.
Here, our aims are (1) to develop and validate a circadian phase prediction method by a transcriptome wide search for a set of features that is capable of predicting the circadian phase of the plasma melatonin rhythm, regardless of sleep history and timing, using a few blood samples; (2) to compare the performance of this method to published approaches; and (3) identify the processes associated with the set of identified biomarkers.
mRNA abundance features and melatonin concentrations were assessed in blood samples collected under four sleep timing conditions (Archer et al., 2014; Möller-Levet et al., 2013) (Figure 1d,h): (i) Sleep in phase with melatonin; (ii) Sleep out of phase with melatonin; (iii) total sleep deprivation, no prior sleep debt and (iv) total sleep deprivation, with prior sleep debt. The melatonin rhythm persists in the absence of sleep and when sleep occurs during the daytime provided that light levels are low, as was the case in these experiments (Figure 1d,h). Each of 678 genome-wide mRNA abundance profiles were assigned a circadian melatonin phase by referencing it to the onset of melatonin secretion (0 degrees) in the corresponding individual time series (Archer et al., 2014; Möller-Levet et al., 2013). Melatonin phases of mRNA samples spanned the entire circadian cycle (Figure 1c,g). To develop practical models, we divided our data set into a training set comprising 329 samples and a validation set comprising 349 samples (see Materials and methods).
Existing methods, such as the molecular timetable method (Ueda et al., 2004) rely on identifying features that correlate well with a 24 hr cosine wave. We estimated rhythmicity of the transcriptome by computing the correlation between a cosine wave and the time series of all features in each of our four data sets and across all four data sets separately for the training and validation set (see Materials and methods). Plotting the r2 value against the rank of the correlation shows that overall 24 hr rhythmicity in the transcriptome is most prominent in the ‘sleep in phase with melatonin’ condition and smallest in the ‘sleep out of phase with melatonin’ condition (Figure 2a,d).
To test the performance of the timetable method (Ueda et al., 2004), we identified features that were robustly rhythmic across all conditions in the training set. Seventy-three features satisfied the criterion of a minimum correlation of 0.3 in each condition and across all four conditions (see Materials and methods and Figure 2—figure supplement 1 for threshold derivation; Figure 2—figure supplement 2 for temporal profiles of genes comprising the timetable model derived from our data). The ‘rhythmicity’ of these 73 features was compared to the ‘rhythmicity’ of previously published phase marker gene sets of (Lech et al., 2016; Hughey et al., 2016) by rank ordering their correlation coefficients. The ‘rhythmicity’ of the 73 identified features is in general, more robust than the rhythmicity of the genes used in previous approaches (Figure 2b,e). In particular, it is the condition ‘sleeping out of phase with melatonin’ (condition ii in Figure 1) (Figure 2b,c,e,f) that appears to exclude many of the previously published phase marker genes as good candidates for features to be included in a timetable model (see Figure 2—figure supplement 3, Figure 2—figure supplement 4, for temporal profiles of genes in the previously reported phase marker gene sets; Supplementary file 1A for their corresponding correlation r coefficients). Indeed, the set of 73 features does not significantly overlap with previously reported sets of phase predictors (Supplementary file 1B for comparison of models’ gene lists). Gene ontology (GO) functional enrichment analysis of the 73 features showed that they were significantly (Benjamini and Hochberg adjusted p-value<0.05) associated with the following biological processes: T cell receptor V(D)J recombination, somatic diversification of T cell receptor genes and somatic recombination of T cell receptor gene segments (Supplementary file 2).
The performance of our one-sample timetable model based on the 73 identified features is summarized and compared to other models in Table 1, Figure 3 and Figure 3—figure supplement 1. The timetable model using these 73 features identified in the training set predicts the circadian phase of samples in the validation set on average within 9 min of the observed melatonin phase which is comparable to the 15 min and −2 min bias when the original gene lists of (Lech et al., 2016) and (Hughey et al., 2016) were used, respectively (Table 1, Figure 3c,f,i). Any identified systematic error in a model can be corrected for when optimizing a final model for future applications and the relevant measure of performance is the standard deviation of the error. The molecular timetable model comprising the 73 features we identified had the least variation of prediction error. The errors of the 73-feature model varied across the circadian cycle (Table 1, Figure 3, Figure 3—figure supplement 1b,e,h). Our one-sample molecular timetable model outperformed the models constructed from published gene sets of Lech et al. (2016) and Hughey et al. (2016), when considering the ‘proportion of samples with ≤2 hr error’ (Table 1, Figure 3j) and the R2 (Table 1, Figure 3—figure supplement 1a,d,g). Nevertheless, it is clear that each of these one-sample molecular timetable models perform poorly when applied to our human transcriptome data.
We next used ‘Zeitzeiger’ (Hughey et al., 2016), which can be considered a machine learning approach to the molecular timetable method, to identify circadian phase biomarkers in our data sets. In Zeitzeiger, a high-dimensional transcriptomic profile sample is first projected into a low-dimensional space using the maximum likelihood principle before assigning the circadian time to that sample. Our implementation of Zeitzeiger identified 107 features as the optimal number with which to predict the circadian phase of one sample (see Materials and methods and Figure 3—figure supplement 2). Only three genes, PER1, NR1D1 and NR1D2, from the original Zeitzeiger gene list identified in animal tissues (Hughey et al., 2016) were also identified as predictors when applied to our human blood transcriptome data. Of the 107 features we identified using Zeitzeiger, 23 were in common with the molecular timetable method (see Supplementary file 1B). The 107 features identified by the Zeitzeiger model did not show any significant functional GO enrichment (Supplementary file 2). The performance of the Zeitzeiger model, when applied to our validation set, is somewhat similar to that of the molecular timetable model across all but two of the performance metrics (Table 1, Figure 3a–f and j). Both methods are similar in their approach, constructing a model based on features identified as being rhythmic within a data set, so this is perhaps expected. Unlike the molecular timetable, errors did not vary significantly across circadian bins. The reduced error of the Zeitzeiger model may be due to inclusion of a greater number of features, their maxima phases providing greater coverage of 360 degrees with which to more accurately predict circadian time (see Figure 3—figure supplement 3).
In view of the suboptimal performance of the approaches described above, we applied another method taking into account the following considerations. We have a much larger number of variables (each of our transcriptomic profiles comprises N (~26 K) quantitative features) than available observations (melatonin phases, one per transcriptomic profile). Furthermore, the variables are not necessarily independent, for example, co-expression of transcripts and/or some features targeting the same gene (transcript). In this instance, a suitable approach to developing a predictive model is to use Partial Least Squares Regression (PLSR) (Boulesteix and Strimmer, 2007). In PLSR, the dimensionality of the predictor set (transcriptome profile x) is reduced by projecting both the predictor set and response variable (melatonin phase y) into orthogonal latent (hidden) spaces, referenced by latent factors, that are relevant for predicting the response variable, without directly prioritizing any underlying time-dependency within the data set (see Materials and methods). Factor loadings, the correlation between each feature and each factor, can then be used to select features a posteriori to produce the circadian phase prediction method. Features with the most weight (loading) across latent factors are of interest and/or the better predictors.
For predictions based on one blood sample, the optimal PLSR model parameters were a combination of five latent factors and 100 mRNA abundance features (see Figure 3—figure supplement 4a and b and Materials and methods; see Supplementary file 1C for the list of features and Supplementary file 1B for the list of 80 corresponding genes). The representation of a given transcriptomic sample in a given latent factor (space) is provided by the latent factor score of that sample, a linear combination of feature values multiplied by the corresponding latent factor loadings. The pseudo time-course of factor scores of the five latent factors (presented in Figure 3—figure supplement 4b) are rhythmic in nature (R2 of fit to a cosine wave range of 0.44–0.51), with phases of the maxima covering 360 degrees (58°, 250°, 56°, 310° and 280°, respectively) and all have a similar amplitude (range of 0.016–0.018). Each latent factor was characterized via functional enrichment analysis of the 200 most weighted features (see Supplementary file 2 for lists).
Latent factor 1, the factor that has the greatest covariation of x and y, gives greatest weight to features significantly enriched for biological processes linked to immune and inflammatory response (e.g. TCF7, LEF1, NFAM1, TFE3, GIMAP5, PRKCA, GAB2, CD28), muscle hypertrophy and cytokine production (e.g. IL13RA1, IL11RA, CCR7, IL1RAP, IL6ST).
Latent factor 2: not significantly enriched for any biological processes.
Latent factor 3: regulation of phosphatidylinositol 3-kinase (ERBB2, KIT, F2R, FLT3, PDGFRB), regulation of interleukin 4 and 6 (GFI1, IL6ST, TCF7, LEF1, XCL1), cell activation (e.g. TTN, SH2D1B, NCR1, PTPN22, PRKCA), and T cell activation (e.g. LEF1, IL6ST, TBX21, VSIG4).
Latent factor 4: circadian rhythm (most significant enrichment and includes TIMELESS, ID2, PROK2, GHRL, PER1, PER2, PER3), response to estrogen (CITED4, GHRL, SOCS1, SOCS3, BGLAP, LDLR, IL1B, DUSP1), response to growth hormone (IRS2, GHRL, SOCS1, SOCS3) and response to corticosteroid stimulus (ANXA3, IL1RN, IL1B, BGLAP, DUSP1, SOCS3). In addition to being enriched for genes involved in circadian rhythmicity, latent factor 4 also contains the temperature-sensitive RNA binding proteins CIRBP and RBM3, which are known to regulate the translation of circadian transcripts (Liu et al., 2013).
Latent factor 5: regulation of plasma lipoprotein clearance (HNRNPK, CD36, LDLR, APOM), regulation of interleukin-12 production (TLR2, TNFSF4, CD36, TIGIT), and regulation of immune effector process (TLR2, CADM1, SHSD1B, TNFSF4, CD36, CLCF1).
Functional enrichment analysis of the 80 genes targeted by the 100 mRNA abundance features identified significantly enriched biological processes associated with leukocyte activation (HSPH1, SH2D1B, IL6ST, CCND3, CD1C, NCR1, TNFSF4, FLT3, IL1B, TXK) and inflammatory response (TLR2, TNFSF4, GHRL, IL1B).
The one-sample PLSR-based model outperformed all other models discussed so far, for most performance metrics (Table 1, Figure 3). The average error is comparable to other models but the distribution of error is narrower and the errors did not vary significantly across circadian bins. Fifty-four percent of samples fell within ≤2 hr of the true circadian phase and the R2 of 0.74 also compares well to other models. Finally, the prediction error of the one-sample PLSR model did not differ significantly across the four conditions. Nevertheless, we aimed to produce a better performing model.
One expects the accuracy of a prediction model to increase as the number of input samples increases. We thus applied the one-sample models to pairs and triplets of consecutive samples. Two-samples performed better than one-sample, and three-samples performed better than two (Table 1, Figure 3j,k and l, Figure 3—figure supplement 5). The PLSR-based model gave the best performance with 71% of samples having a prediction error of ≤2 hr and a R2 of 0.88 using three consecutive samples. However, when taking three consecutive samples as input, the model performed significantly worse for samples derived from the condition ‘sleeping out of phase with melatonin’ (Tukey HSD post-hoc test, p-value<0.05).
An alternative approach is to develop a model which takes multiple samples as input and utilizes the difference in mRNA abundance between input samples to predict the phases instead of averaging independent predictions. Thus, it is the difference in the ‘state’ of the transcriptome, without the significant ‘trait-like’ characteristic of the transcriptome (De Boever et al., 2014) that predicts the phase of an individual. Here, we developed a two-sample based model that takes two-samples taken 12 hr apart in the same participant during the same condition as input. We did this for both the molecular timetable and PLSR approaches (see Materials and methods).
Our trained two-sample differential PLSR model comprises 5 latent factors and 100 features (see Supplementary file 1D for the list of features and Supplementary file 1B for the list of corresponding genes). The pseudo time-course of factor scores of the five latent factors (Figure 3—figure supplement 4d) shows that all are rhythmic (R2 of fit to a cosine wave range of 0.71–0.8), with similar amplitude (range of 0.021–0.022) and maxima phases spanning 360 degrees (250°, 310°, 330°, 320°, 44°, respectively). Functional enrichment analyses of the top 200 features that carry the most weight in each latent factor identified significantly enriched biological processes (see Supplementary file 2) including:
Latent factor 1: no significant enrichment.
Latent factor 2: regulation of metabolic process (e.g. HSPH1, ANXA3, PPARGC1B, LDLR, FOSL2, CTNNBIP1, BCL6, PIK3IP1), regulation of acute inflammatory response (PTGS2, IL1B, IL6ST, OSM), response to glucocorticoid stimulus (ANXA3, PTGS2, IL1RN, IL1B, PPARGC1B, JUNB, SOCS3) and circadian rhythm (NR1D1, ARNTL, PER1, PER2, PER3).
Latent factor 3: negative regulation of cellular and metabolic processes (e.g. DYNLL1, ID2, NEDD4L, DHCR24, DUSP1, NR2F1, AVIL, TLR2), lymphocyte activation (IRS2, CXCR5, HSPH1, KIF13B, CCND3, ID2, BCL3, NDFIP1, FLT3, IL1B, TXK, CD83), rhythmic process (GHRL, TEF, ID2, TIMELESS, ADAMTS1, PROK2, PER1, PER3) and apoptotic signaling in response to DNA damage (HIC1, FNIP2, BCL3, DDIT4).
Latent factor 4: positive regulation of cellular processes (e.g. CLEC4E, INSL3, SLC1A3, ARG1, RBM14, RBM3, TRAF4), developmental process (e.g. CLIC5, SSH2, SPIB, RHBDD1, KLF9, FKBP4, KCNE3, SPATA6), response to lipid (FADS1, TLR2, SLC30A4, CITED4, GHRL, BGLAP, NFKBIA, IRAK3, LDLR, GATA3, IL1RN, REST, ARG1, IL1B, DUSP1) and circadian rhythm (TIMELESS, GHRL, PTEN, PER1, PER2, PER3).
Latent factor 5: regulation of response to stimulus (e.g. CDKN1C, INPP5B, OPA1, MLC1, DVL2, APCDD1, VSIG4, CACNG4), regulation of immune system process (KLRD1, KIR2DL2, SH2D1B, A2M, CMKLR1, NCR1, PDE5A, KIR2DS2, ADA, CCL4, ADORA2A, TBX21, VSIG4, KLRG1, RELA, NR1D1, APOBEC3G, FCN2), rhythmic process (ADAMTS1, FANCG, NR1D1, NPAS2, BHLHE41, PER1, ADA) and calcium-mediated signaling (NMUR1, CHPS, CCL4, ADA, CMKLR1).
For the two-sample differential PLSR, four of the five latent factors contained genes enriched for circadian rhythm or rhythmic process. Each factor contained at least three clock genes, with factor two containing six genes (ARNTL, NR1D1, NR1D2, PER1, PER2, PER3) together with CIRBP and RBM3. Aspects of immune function are also common to the factors, but factor five in particular showed enrichment for immune-cell-specific genes (e.g. immune killer cell genes KIR2DL2, KIR2DL5A, KIR2DS2, KIR2DS4, KLRD1, KLRF1, KLRG1), in addition to the interesting inclusion of adenosine deaminase (ADA) and the adenosine A2a receptor (ADORA2A).
Forty-seven (62%) of the 76 unique genes targeted by the features identified using the two-sample approach also appear in the set of 80 unique genes identified by our one-sample approach (Supplementary file 1B). Functional enrichment analysis of the 76 genes identified significantly enriched biological processes including; immune response, circadian and rhythmic processes and regulation of cellular processes. Unlike the one-sample method, the 76-gene list comprising the two-sample differential PLSR model contained all three PER genes (together with NR1D2), and several genes with immune-cell-specific functions (e.g. KLRF1, KLRD1, CD83, IL1RN) .
Both two-sample differential models outperform their preceding one-sample models and multiple consecutive-sample applications (Figure 4, Figure 3—figure supplement 5, Table 1). The PLSR-based model performs well, now able to predict within 20 min on average, with an SD of 1 hr 41 min, without any significant difference in error distribution across circadian phase bins. Significantly, the two-sample differential model is able to correctly predict 82% of samples within 2 hr, with an R2 of 0.9. Similar to the one-sample model, we find that the two-sample differential model is able to predict baseline conditions with higher accuracy than the molecular timetable equivalent. However, unlike the one-sample PLSR model, we find a significant difference in error distribution between conditions, where ‘Sleep out of phase with melatonin’ has a significantly larger mean error than the other three conditions (Tukey-HSD post-hoc test p-value<0.05). Nevertheless, the difference in mean error between conditions for our two-sample differential model is less than that observed for the one-sample PLSR model when applied to two or three consecutive samples.
The human blood transcriptome is in part under circadian control and potentially contains information about the phase of circadian rhythmicity in many organs. The analyses presented here show that a few transcriptome samples contain enough information to predict the phase of the melatonin rhythm. This remains true even under conditions in which the sleep-wake cycle, which also influences the blood transcriptome, is desynchronized from the melatonin rhythm. Melatonin phase and shifts in melatonin phase correlate well with other variables driven by the SCN, for example, cortisol (Dijk et al., 2012). In clinical practice, it is therefore assumed that the melatonin phase reflects the status of the SCN with minimal error, although the actual error in SCN phase measurement by melatonin is not known. The methods developed here may aid the routine assessment of circadian melatonin phase in clinical practice and benefit the diagnosis and treatment of sleep and circadian disorders (Mullington et al., 2016).
We consider our PLSR-based approach the first method for developing a melatonin phase predictor without a priori selection of features or feature characteristics: no a priori information is used to supervise model development. We do not select for 24 hr rhythmically expressed genes that are either known in the literature (Lech et al., 2016) or specifically identified within the data (Ueda et al., 2004; Hughey et al., 2016). The validation analyses show that our model comprising features selected a posteriori outperforms ‘a priori’ approaches. Directly comparing our ‘a posteriori’ method with ‘a priori’ methods (Ueda et al., 2004; Hughey et al., 2016), by deploying and testing each using the same training and validation data, shows that our method is valid and capable of generating superior circadian phase prediction models. It should be noted that many of the genes identified by PLSR are not classical clock genes or clock-controlled genes (see Supplementary file 1B). Nevertheless, 68% of the genes forming the one-sample, and 63% genes forming the two-sample differential PLSR-based model, we have previously identified as being rhythmically expressed in at least one of the four conditions we have analyzed (Archer et al., 2014; Möller-Levet et al., 2013; Laing et al., 2015). Only 18% and 12% of the one and two sample(s)’ genes (respectively) were found to be rhythmically expressed in the condition ‘sleep out of phase with melatonin’. In our previous analyses, we used a prevalence-based classification which is a fundamentally different method to the rhythmic detection strategies of the molecular timetable and Zeitzeiger methods (Ueda et al., 2004; Hughey et al., 2016). However, comparing across prediction models it is clear that >30% (34% and 31% for the one-sample and two-sample differential models, respectively) of the gene lists forming the PLSR-based models are unique to the approach (compared to the molecular timetable and Zeitzeiger approaches applied to the same training data set) (Figure 5).
We used OmicCircos (Hu et al., 2014) to visualize and compare the overlap of gene features that were identified by each of the four models trained on the same data (Figure 5). Six genes were identified in all four models: PER1 (circadian clock gene, induced by glucocorticoid [Cuesta et al., 2015]), SMAP2 (involved in vesicle trafficking in the trans-golgi network [Funaki et al., 2013]), DDIT4 (part of mTOR signaling pathway [Polman et al., 2012]), NCR1 (natural killer cell receptor [Lai and Mager, 2012]), ZBTB16 (transcription factor that regulates glucocorticoid-induced gene expression [Wasim et al., 2010]) and FKBP5 (chaperone that plays a key role in the regulation of glucocorticoid receptor signaling [Binder, 2009]). Apart from SMAP2 and NCR1, all these genes are involved in glucocorticoid signaling or are regulated by glucocorticoid. Immune function and glucocorticoid signaling are themes that link all four models. More than one third of all the genes in the four models are associated with glucocorticoid signaling (molecular timetable 44%, Zeitzeiger 36%, one-sample PLSR 36%, two-sample differential PLSR 38%) (Nehmé et al., 2009).
The gene lists of the molecular timetable and Zeitzeiger models show little overlap with those of the PLSR models, which have a high level of similarity (62%). Fifteen of the top-25 genes (by weighting) in the PLSR gene lists are common to both models (PER1, DDIT4, FKBP5, ZNF397, TSPAN4, C20orf112, MIA, AVIL, ZBTB16, SPATA6, A_24_P178167, BQ706262, SMAP2, FLT3, HSPH1). Several of these genes have related functions: MIA and TSPAN4 bind integrin and regulate cell motility (Bosserhoff, 2005; Kashef et al., 2013), while AVIL binds actin and regulates cell outgrowth (Tümer et al., 2002). Twelve of the genes in the top-25 lists are associated with glucocorticoid signaling and form part of a linked glucocorticoid signaling network (Figure 5—figure supplement 1). The top four features in the PLSR two-sample differential model gene list are present in all four models (DDIT4, PER1, FKBP5 and ZBTB16).
As discussed, the PSLR approach does not pre-select features with respect to rhythmicity, and indeed cannot be used to detect rhythmicity. Only the feature lists identified by the two sample differential PSLR model appear enriched for rhythmic, circadian processes. This may be because in the differential two-sample approach individual differences in average levels of expression (i.e. trait like characteristics of the transcriptome [De Boever et al., 2014]) are removed. This apparently emphasizes features with a 24 hr variation around a baseline.
Both the cortisol and the melatonin rhythm are driven by the SCN. The prominence of glucocorticoid signaling in the lists may reflect a more extensive influence of cortisol on the peripheral transcriptome. Overall the features of the PLSR models seem to focus more on genes with functions specifically related to blood such as glucocorticoid signaling and immune cell-specific roles. Thus, the features within the blood transcriptome that are most predictive of the melatonin phase are closely related to blood-cell-specific functions. We have previously shown that the influence of the sleep-wake cycle and circadian rhythmicity varies across the transcriptome transcripts in blood (Archer et al., 2014). Whether and how these blood-cell-specific functions remain in phase with the SCN, even when the sleep-wake cycle is shifted, remains unclear. It may be that those that predict SCN phase are regulated by local circadian clocks, which in turn are synchronized by systemic cues from the central circadian clock, such as cortisol. Alternatively, these SCN phase predictors in blood may be directly driven by systemic cues such as cortisol or melatonin. As far as melatonin is concerned, only very few melatonin-related genes appeared in any of the models. The circadian expression of pineal HNRNP proteins regulates the expression of the rate-limiting melatonin synthesis enzyme AANAT (Kim et al., 2005). In this study, HNRPDL is a feature of each model except the PLSR two-sample differential model and was previously identified as remaining rhythmic in all four of the sleep conditions that have contributed data to this study (Laing et al., 2015). Thus, while glucocorticoid signaling features strongly in the models developed here, it is also possible that melatonin signaling contributes to phase entrainment in whole blood, possibly also locally within blood cells.
The enhanced power of our PLSR-based models compared to others cannot be due to the unique combination of experimental conditions that we have included in our training and validation sets because all models were tested on the same data sets. Our training and validation data are limited; all samples were collected in darkness and dim light and from young healthy individuals. Comparing the performance of the various models across conditions emphasizes the importance of inclusion of conditions in which synchrony between various organs is challenged. Several models that performed well under normal conditions performed less well under desynchronized conditions. Indeed, we expect the conditional dependency of performance for the two-sample differential models to be due to the significant temporal disruption evoked by sleeping out of phase with melatonin, where the difference in mRNA abundance between two samples taken 12 hr apart may not be so well preserved as in baseline and/or sleep deprivation conditions. Thus, inclusion of similar conditions will be important for the validation for any blood-based biomarker of organ-specific rhythmicity and the generation of clinically relevant applications in ‘real-world’ scenarios. Patients may be suffering from sleep deprivation or sleep and food-intake desynchronized from the circadian system, conditions in which rhythmic molecular profiles are disrupted. A ‘final’ model to be applied in the clinic should also be trained on real conditions in which light exposure and food availability are not controlled to ensure its validity to other sections of the population. This could include not only circadian rhythm sleep-wake disorder patients but also other populations in which melatonin phase and or cortisol phase are disrupted, such as in Smith-Magenis syndrome (De Leersnyder, 2013) or in Addison’s disease (Oster et al., 2016).
Several non-omics-based methods to assess circadian phase are currently in use. In the Dim-Light-Melatonin Onset (DLMO) Protocol, many samples (6–18) are collected during an approximately 6–12 hr period timed to be close to the expected onset of melatonin secretion (Burgess et al., 2016). Thus, unlike the method presented here, the DLMO protocol requires prior knowledge of the approximate circadian phase of the pacemaker to identify the phase of the SCN (with unknown error). In fact, one immediate application of the current model is to use it to identify the time window during which a DLMO assessment should be scheduled. Alternative assessments of SCN circadian phase, when there is no prior knowledge about the likely circadian phase, requires measurements covering at least one circadian cycle. In field situations, this can be accomplished by the scheduled collection of urine for determination of 6-sulfatoxy-melatonin, where the 95% CIs of such assessment are thought to be ±1 hr 30 min (Lushington et al., 1996). Yet, even when urine samples are collected very frequently (e.g. every 4 hr), the correlation between the urine-based phase marker 6-sulfatoxy-melatonin and the gold-standard plasma-based phase marker melatonin is only 0.75 (Gunn et al., 2016). Similarly, accurate assessment of circadian phase from core body temperature in the presence of a sleep-wake cycle is problematic because sleep and wakefulness have a direct effect on core body temperature, which depends on circadian phase. Therefore, extraction of phase of the circadian pacemaker from core body temperature requires either recordings over several days or weeks, or 24–40 hr T recordings in specialized laboratory conditions.
From a practical perspective, a circadian phase predictor that requires only one sample is ideal, minimizing the total time in the clinic and overall cost. However, the marked increase in performance offered by taking multiple samples, either two or three consecutive samples 3 or 4 hr apart, or two-samples 12 hr apart for differential analysis, should also be considered. Certainly, the multiple-sample models/applications we tested here rely on far fewer samples than 6-sulfatoxy-melatonin and/or core body temperature, even plasma melatonin assessments, and can be completed during one ‘out-patient’ clinic visit. Our PLSR-based models and applications can predict more than 80% of circadian phases within 2 hr of the observed phase when using the differential analysis of two-samples as input.
The current models are based on whole-blood transcriptomic profiles comprising mRNA abundance features measured by microarrays. However, the methodology we present can be readily applied to map any high-dimensional molecular abundance ‘feature’ profile, such as those generated by microarray, RNA sequencing and metabolomics technologies, to any specified measure of circadian phase. Although we used melatonin/SCN phase as an example with obvious applications in circadian rhythm sleep-wake disorders, one can imagine in the era of precision medicine the development of specific PLSR-based circadian phase prediction models for the liver and/or other specific tissues and organs. All these organs may be targeted by chronotherapy (Zhang et al., 2014). All that is required is the identification of specific phase markers for the tissue/organ of clinical interest.
Transcriptome-based measures appear to suffer from less noise than other ‘-omes’ (Su et al., 2011; Tachibana, 2014). The 100 transcriptomic features we have identified as predictors of SCN circadian phase using whole-blood can be simultaneously measured by a cost-effective low-resolution molecular technique. Therefore, the PLSR-based models we present here offer a viable solution for a clinical platform (Mullington et al., 2016); even more so when redundant features that target the same gene/transcript, are removed.
Data were collected during two large-scale circadian experiments to investigate the effects of mistimed sleep and total and repeated partial sleep deprivation on the temporal organization of the human blood transcriptome (Archer et al., 2014; Möller-Levet et al., 2013). The protocols received a favorable opinion from the University of Surrey’s Ethics committee and all participants provided written informed consent. In each experiment, we obtained from the same participant, blood samples (every 3 or 4 hr) covering at least one circadian cycle to measure mRNA abundance and hourly blood samples for measuring melatonin rhythms. This allowed assignment of a circadian melatonin phase to each mRNA sample. Melatonin concentrations were determined with a radioimmunoassay (Stockgrand Ltd, Guildford, UK [Fraser et al., 1983]). Circadian phase was defined as the time when melatonin levels reached 25% of the melatonin amplitude and is referred to as the Dim Light Melatonin Onset (DLMO). DLMO = (((melatonin amplitude - baseline melatonin level)*25)/100) + baseline melatonin level (Hasan et al., 2012; Lo et al., 2012).
We collected 678 samples for which we have both the melatonin phase and mRNA abundance profile from 53 participants. Fifty-nine percent (397) of these samples are from male participants (see Supplementary file 3 for demographics, including ethnicity). The participant population was stratified according to the PER3 rs57875989 genotype (38% PER34/4; 35% PER35/5; 27% PER34/5). The complete time-series of samples from homozygote participants were hybridized to microarrays (i.e. 14–20 samples per participant; 7–10 samples per sleep condition). Select RNA samples (one to three samples per sleep condition) from the complete time-series of samples for heterozygote individuals were hybridized to microarrays. All transcriptomic data are accessible from the Gene Expression Omnibus (Barrett et al., 2013) (RRID:SCR_007303, via the Accession numbers: GSE48113 [PER3 homozygote, PER34/4 and PER35/5] and GSE82113 [PER3 heterozygote, PER34/5] data for mistimed sleep, and GSE39445 (homozygote) and GSE82114 (heterozygote) data for sleep deprivation and training) and the data sources at http://sleep-sysbio.fhms.surrey.ac.uk/PLSR_16/.
Samples were split into a training set to develop the models, and a validation set to test the predictive performance of the models. Factors that were balanced between the training and validations sets were: the experimental protocol (mistimed sleep or sleep deprivation), PER3 genotype, sex and circadian phase (as indicated in Supplementary file 3). All samples pertaining to a participant were assigned to the same data set. Data of 26 participants, comprising 329 samples were randomly assigned to the training set and the data of 27 participants, comprising 349 samples to the validation set. To facilitate the development and benchmarking of future models, the data matrices for the training and validation data sets are provided at: http://sleep-sysbio.fhms.surrey.ac.uk/PLSR_16/.
For models developed and tested for a ‘one sample’ application, including the prediction of consecutive samples requiring the application of a one-sample model multiple times, all 329 and 349 samples of the training and validation sets were used. For ‘two-sample differential’ models the mRNA abundance values for a pair of samples from the same participant during the same condition, separated by 12 hr, were reduced to a single observation. The single observation represents the difference in the mRNA abundance, that is, a vector obtained as a result of computing Sample one minus Sample 2 (taken 12 hr later than Sample 1). In total, the ‘two sample differential’ training and validation data sets comprised 163 and 180 mRNA abundance difference ‘samples’, respectively. The data matrices for the training and validation data sets are provided at: http://sleep-sysbio.fhms.surrey.ac.uk/PLSR_16/. Note that the two-sample differential molecular timetable model is based on the one-sample molecular timetable. For this two-sample model, the difference in the timetable templates with a 12 hr time lapse is used to predict the difference in mRNA abundance between two samples taken 12 hr apart.
For the training set the Log2 mRNA abundance values were processed by applying quantile-normalization using the R Bioconductor package limma (Gentleman et al., 2004; Smyth, 2005) (RRID:SCR_010943) followed by batch correction using the ComBat function of the R package sva (Leek et al., 2012; Johnson et al., 2007) (RRID:SCR_012836). For the validation set, the Log2 mRNA abundance values were quantile-normalized using the same reference array of empirical quantiles used to normalize the samples included in the training set. No batch correction was applied to the validation data set to reflect the type of data that would be provided in a clinical setting (i.e. one or two samples processed in different batches to the original training data used to build the model). Non-control technically replicated probes (mRNA abundance features), along with their corresponding Agilent feature flags (Agilent Feature Extraction Software [v10.7], Reference Guide, publication number G4460-90036) were averaged. Individual samples with more than 30% of features flagged were excluded and features (probes) flagged in more than 10% of the samples were filtered out.
Assuming angles lay on a unit circle (radius = 1), an angle can be converted to Cartesian coordinates using the trigonometric functions of and . The Cartesian coordinates can be converted back to an angle using the arctan function for in quadrant I of the coordinate plane, for in quadrant II or III, and for in quadrant IV of the coordinate plane. In this way, the mean value () of a set of n angles () can be calculated as
where and if is in quadrant I of the coordinate plane, if is in quadrant II or III and if is in quadrant IV.
The difference between angles and can be calculated using
For all operations, angles were converted to Cartesian coordinates (cosine , sine ) to perform arithmetic operations.
To assess the overall performance and clinical relevance of our models, we used several metrics:
The frequency distribution of error between the predicted circadian phase and the experimentally observed circadian phase , when binned in 15° (60 min) circadian phase bins. The distribution is summarized by the mean and standard deviation.
The fit of predicted phases to the observed phase via R2. If is the observed circadian phase and the predicted circadian phase of the th sample, R2 is defined as:(2)
where is the directional mean of the observed circadian phases, RSS: regression sum of squares; TSS: total sum of squares. Thus, R2 compares the variance in prediction error to the variance in the observed circadian phases, giving a fraction of unexplained variance. When subtracted from 1, the R2 value represents the variance that can be explained by the model. Where RSS (variance in prediction error) is greater than TSS, the resultant R2 will have a negative value, corresponding to a model that performs worse than fitting a horizontal line.
The error distribution across the sampling period of 24 hr (360 degrees), where a better performing model is one that exhibits no (or least) dependency on the time at which the sample is taken. Differences in error across time was assessed by ANOVA (Analysis of Variance) analysis for the factor of ‘sampling time bin’ using the lm and anova functions in R (RRID:SCR_001905); a p-value<0.05 suggests that the prediction accuracy of the model is dependent on the circadian phase at which the sample is taken.
The cumulative frequency distribution of absolute error, to identify the proportion of samples that can be predicted within a given error range, where the highest proportion combined with smallest error range denotes the better prediction model.
To identify parameters and/or thresholds with which to build a model and/or assess the overall performance of the final model, we used the performance metric R2. The exception to this is the Zeitzeiger model which relies on the Mean Average Error to identify parameters (see below).
Differences in circadian phase prediction errors across the four conditions were assessed using the ‘aov’ function of R followed by Tukey Honest Significant Difference (HSD) post-hoc test (R function TukeyHSD) (RRID:SCR_001905).
To avoid over-fitting during model development, we performed leave-one-participant-out cross-validation. For p participants in the original one or two-sample training set, all samples associated with p-1 participant were used for initial training of the model, and the phase of all samples associated with the excluded participant were used to test the model. Once all p-1 sets were used to create and test the model the predicted phases were used to determine the predictive power of the cross-validated model, as described above.
All code used to build and test the models, calculate the given statistics and generate the provided figures, along with a description on their usage can be found at: http://sleep-sysbio.fhms.surrey.ac.uk/PLSR_16/.
We implemented the molecular timetable method described in Ueda et al. (2004). Instead of using a summary metric (e.g. mean or median) or model to represent time-series across participants, we used each individual sample within the four experimental conditions in the training set to generate a single, densely-sampled time-series per mRNA abundance feature per experimental condition. We then created a fifth, ‘all conditions’ time-series, by using each individual sample in the training set irrespective of experimental condition. Subsequently, we selected the mRNA abundance features with correlation values to the cosine templates above the cutoff correlation value of 0.3 for all five time-series. The value of 0.3 was identified as the optimal threshold as it produced the highest R2 value following leave-one-participant-out cross-validation using different correlation thresholds. Figure 2—figure supplement 1 depicts the change in model performance (R2 values) at different correlation thresholds. It was not necessary to use the coefficient of variation to further filter the time-indicating features given that only a small fraction of all features have a detectable circadian oscillation under all experimental conditions included in our analyses. Where multiple mRNA features targeting the same gene met our correlation threshold, the best performing (most correlated to a cosine) mRNA feature was selected to represent a gene. Subsequently, the method identified 73 mRNA abundance features, which target 73 unique genes, as time-indicating features.
The timetable model is essentially a look-up table. To predict the phase of a single sample, the set of mean-centred mRNA abundance values for the 73 features comprising the model are extracted and compared to the expected set of values at different phases within the molecular timetable model (look-up table). The phase producing the highest correlation is the sample’s predicted melatonin phase. In a similar approach, the resulting molecular timetable was further used to construct a ‘two-sample’ molecular timetable, where the difference in mRNA abundance of two samples taken 12 hr apart is compared to the difference in the timetable templates with the same time lapse.
We used the Zeitzeiger method of (Hughey et al., 2016) implemented in the zeitzeiger R package, downloaded from https://github.com/jakejh/zeitzeiger (zeitzeiger_0.0.0.9000, RRID:SCR_014791). To build a Zeitzeiger predictor, given a training data set, there are two parameters that need to be optimized: ‘sumabsv’, which controls for the overfitting of the model by controlling the number of features (probes measuring mRNA abundance) to be used to form principal components, and ‘nSPC’, which controls the number of principal components that describe how the features change over time. To identify the best parameters for training the predictor, we ran leave-one-participant-out cross-validation of the training set. The cross-validation was run over seven values of sumabsv (1, 1.5, 2, 2.5, 3, 3.5 and 4) and five values of nSPC (1, 2, 3, 4 and 5). To compare the performance of the predictors generated from each set of parameter values, we calculated the Mean Average Error (MAE). As seen in Figure 3—figure supplement 2, the smallest MAE is produced by the combination of sumabsv = 4 and nSpc = 3. These optimal values were then used to produce a Zeitzeiger based model to predict the circadian phase of samples in the validation set, with which we could assess the overall performance of the model and directly compare with other models. The 107 mRNA abundance features (probes) that can predict circadian phase based on the Zeitzeiger approach were identified as those carrying weight in at least one of the three principal components used to construct the model.
In this work, is a matrix of samples by mRNA abundance features and the vectors are the corresponding observed melatonin phase of each sample represented in Cartesian coordinates. Partial Least Squares Regression (PLSR) is a linear regression model that regresses on to by projecting both and into a reduced-dimensional space (so called latent factor space) that maximizes their covariance before performing a least squares regression. This projection to a latent space enables PLSR to handle the observed co-linearity between mRNA abundance features derived from the same gene, that is, multiple probes targeting the same transcript and/or co-expressed ‘modules’. To build a PLSR model, given a training data set, there are two parameters that need to be optimized: (1) the number of latent factors, that is dimensions, and (2) the number of mRNA abundance features . If the values of these parameters are set too high (e.g. maximum value), there is a danger of over-fitting the model to the training set, and if the parameters are set too low, the model would not be sensitive nor specific enough to provide high accuracy. Similar to the training in the molecular timetable and Zeitzeiger methods, above, we used leave-one-participant-out cross-validation of the training set to identify the best combination of model parameters, in this case the number of latent factors and number of mRNA abundance features. For our one-sample PLSR approach, we compared the cross-validated performance of 136 combinations of 5 to 40 latent factors (in increments of 5) and 100 to 5000 features with which to build a model (see Figure 3—figure supplement 4a). For our two-sample PLSR approach, we compared the cross-validated performance of 100 combinations of four to eight latent factors and 50 to 1000 features (Figure 3—figure supplement 4c). To generate each model, the plsregress function of Matlab (184.108.40.2061 [R2013b]) was invoked twice. The first call to the plsregress function generates a PLSR-based model given the z-scored mRNA abundance feature values (for all features) of the training set (where samples for a single-participant have been removed), the corresponding melatonin phase values represented in Cartesian coordinates, and the number of specified latent factors as input. The absolute loadings (weights) of this model are then used to identify the top features to create a subset of the training set comprising the features. This subset of the training set used in the first call is then used as input to the subsequent second call of the plsregress function given the same (as the preceding call to the function) melatonin phase values represented in Cartesian coordinates, and the number of specified latent factors as input. The model derived from the second call to the plsregress function is then used to predict the melatonin phase of the left out (cross-validation) samples. Note that each p-trained model during cross-validation may have a different set of ‘top’ (most weighted) features as feature order was not fixed in the cross-validation. The performance of a combination of latent factors and features is then assessed by R2 between the predicted and true melatonin phase values across all ‘leave-one-out’ predictions from models comprising latent factors and features. The overall optimal combination of parameters was defined as the combination that produced an R2 >= 70% using the least number of features and latent factors. Figure 3—figure supplement 4 depicts the changes in performance (R2) of the model as the number of features and latent factors change when using one sample (Figure 3—figure supplement 4a) or two samples (Figure 3—figure supplement 4c) as input. As can be seen in Figure 3—figure supplement 4a, when predicting circadian phase using one transcriptomic sample one of the highest R2 values, at 73%, occurs when using 100 mRNA features and five latent factors. The maximum R2 value of 76% is obtained by increasing the features to 2000 and 10 latent factors. However, given that this number of features could induce overfitting and affect clinical applications (e.g. too many features to measure in a single platform) without a significant increase in R2 we chose the combination of 100 mRNA features and five latent factors as an acceptable combination of parameters with which to build the final model. Using the same approach, for two samples (see Figure 3—figure supplement 4c), we identified 100 features and five latent factors (R2 value of 93%) as acceptable parameter values with which to train the final predictive model.
Provided that a future user has blood transcriptome mRNA abundance values from the same microarray platform (GEO [Barrett et al., 2013; RRID:SCR_007303] accession number GPL15331) and/or features, our models can be used to predict the circadian phase of the participant from which the sample(s) were taken. The raw Log2 mRNA abundance values of the test sample(s) will need to be quantile-normalized using the same reference array of empirical quantiles used to normalize the samples included in our training set (see data at http://sleep-sysbio.fhms.surrey.ac.uk/PLSR_16/).
The PLSR method projects both the predicted and the observed variables onto a new space to find a linear regression model. The linear regression model is used to predict circadian phase () from a set of mRNA abundance values (), where is calculated from the predicted Cartesian coordinates of (). Thus, the model parameters are the regression coefficients () and the intercept value () for the cosine of the predicted phase (), such that , and the regression coefficients () and the intercept value () for the sine of the predicted phase (), such that The predicted phase is then , where if () is in quadrant I of the coordinate plane, if () is in quadrant II or III and if () is in quadrant IV.
A circadian phase prediction can therefore be obtained by computing the sum of the products of the normalized z-scored mRNA abundance values and their corresponding regression coefficients and adding the intercept values (see data at http://sleep-sysbio.fhms.surrey.ac.uk/PLSR_16/).
The mean true circadian phase of two or three consecutive samples was directly compared with the mean predicted circadian phase of two or three consecutive samples, as calculated using a one-sample model, for example, molecular timetable, Zeitzeiger and PLSR.
Functional enrichment analysis was performed using the online tool Webgestalt (Wang et al., 2013) (RRID:SCR_006786). Here, a given mRNA abundance feature list was converted to a list of unique HGNC gene symbols and uploaded to the Webgestalt tool. Gene Ontology (GO) enrichment analysis was performed using the reference human genome as the background and selecting for significant (Benjamini and Hochberg adjusted p-value<0.05) GO biological processes and GO molecular functions.
Dual-specificity phosphatase 1: a critical regulator of innate immune responsesBiochemical Society Transactions 34:1018–1023.https://doi.org/10.1042/BST0341018
Melatonin: characteristics, concerns, and prospectsJournal of Biological Rhythms 20:291–303.https://doi.org/10.1177/0748730405277492
FMS-like tyrosine kinase 3 interacts with the glucocorticoid receptor complex and affects glucocorticoid dependent signalingBiochemical and Biophysical Research Communications 368:569–574.https://doi.org/10.1016/j.bbrc.2008.01.146
Clinical practice guideline for the treatment of intrinsic circadian rhythm Sleep-Wake disorders: Advanced Sleep-Wake phase disorder (ASWPD), Delayed Sleep-Wake Phase Disorder (DSWPD), Non-24-Hour Sleep-Wake Rhythm Disorder (N24SWD), and Irregular Sleep-Wake Rhythm Disorder (ISWRD). An Update for 2015: An American Academy of Sleep Medicine Clinical Practice GuidelineJournal of Clinical Sleep Medicine 11:1199–1236.https://doi.org/10.5664/jcsm.5100
NCBI GEO: archive for functional genomics data sets--updateNucleic Acids Research 41:D991–D995.https://doi.org/10.1093/nar/gks1193
Glucocorticoids increase C/EBPbeta activity in the lung epithelium via phosphorylationBiochemical and Biophysical Research Communications 334:638–645.https://doi.org/10.1016/j.bbrc.2005.06.146
The role of FKBP5, a co-chaperone of the glucocorticoid receptor in the pathogenesis and therapy of affective and anxiety disordersPsychoneuroendocrinology 34 Suppl 1:S186–S195.https://doi.org/10.1016/j.psyneuen.2009.05.021
Melanoma inhibitory activity (MIA): an important molecule in melanoma development and progressionPigment Cell Research 18:051012082332003–051012082332416.https://doi.org/10.1111/j.1600-0749.2005.00274.x
Partial least squares: a versatile tool for the analysis of high-dimensional genomic dataBriefings in Bioinformatics 8:32–44.https://doi.org/10.1093/bib/bbl016
Home dim light melatonin onsets with measures of compliance in delayed sleep phase disorderJournal of Sleep Research 25:314–317.https://doi.org/10.1111/jsr.12384
Glucocorticoids entrain molecular clock components in human peripheral cellsThe FASEB Journal 29:1360–1370.https://doi.org/10.1096/fj.14-265686
Dosing-Time makes the poison: Circadian regulation and pharmacotherapyTrends in Molecular Medicine 22:430–445.https://doi.org/10.1016/j.molmed.2016.03.004
The hockey-stick method to estimate evening dim light melatonin onset (DLMO) in humansChronobiology International 31:349–355.https://doi.org/10.3109/07420528.2013.855226
Smith-Magenis syndromeHandbook of Clinical Neurology 111:295–296.https://doi.org/10.1016/B978-0-444-52891-9.00034-8
Circadian rhythm disorders and melatonin production in 127 blind women with and without light perceptionJournal of Biological Rhythms 29:215–224.https://doi.org/10.1177/0748730414536852
Rhythmic serotonin N-acetyltransferase mRNA degradation is essential for the maintenance of its circadian oscillationMolecular and Cellular Biology 25:3232–3246.https://doi.org/10.1128/MCB.25.8.3232-3246.2005
Comparisons of the variability of three markers of the human circadian pacemakerJournal of Biological Rhythms 17:181–193.https://doi.org/10.1177/074873002129002474
Evaluation of mRNA markers for estimating blood deposition time: Towards alibi testing from human forensic stains with rhythmic biomarkersForensic Science International: Genetics 21:119–125.https://doi.org/10.1016/j.fsigen.2015.12.008
Adrenal clocks and the role of adrenal hormones in the regulation of circadian physiologyJournal of Biological Rhythms 30:20–34.https://doi.org/10.1177/0748730414553971
Melatonin entrains free-running blind people according to a physiological dose-response curveChronobiology International 22:1093–1106.https://doi.org/10.1080/07420520500398064
Urinary 6-sulfatoxymelatonin cycle-to-cycle variabilityChronobiology International 13:411–421.https://doi.org/10.3109/07420529609020912
Central and peripheral circadian clocks in mammalsAnnual Review of Neuroscience 35:445–462.https://doi.org/10.1146/annurev-neuro-060909-153128
Dexamethasone inhibits chondrocyte differentiation by suppression of wnt/β-catenin signaling in the chondrogenic cell line ATDC5Histochemistry and Cell Biology 144:261–272.https://doi.org/10.1007/s00418-015-1334-2
Molecular architecture of the mammalian circadian clockTrends in Cell Biology 24:90–99.https://doi.org/10.1016/j.tcb.2013.07.002
Melatonin: both master clock output and internal time-giver in the circadian clocks networkJournal of Physiology-Paris 105:170–182.https://doi.org/10.1016/j.jphysparis.2011.07.001
Circadian rhythm abnormalities in totally blind people: incidence and clinical significanceThe Journal of clinical endocrinology and metabolism 75:127–134.https://doi.org/10.1210/jcem.75.1.1619000
Circadian rhythm sleep disorders: lessons from the blindSleep Medicine Reviews 5:189–206.https://doi.org/10.1053/smrv.2000.0147
Bioinformatics and Computational Biology Solutions Using R and Bioconductor (Statistics for Biology and Health)Springer.
Human phase response curve to a 1 H pulse of bright white lightThe Journal of Physiology 590:3035–3045.https://doi.org/10.1113/jphysiol.2012.227892
Integrated metabolome and transcriptome analysis of the NCI60 datasetBMC Bioinformatics 12 Suppl 1:S36.https://doi.org/10.1186/1471-2105-12-S1-S36
WEB-based GEne SeT AnaLysis toolkit (WebGestalt): update 2013Nucleic Acids Research 41:W77–W83.https://doi.org/10.1093/nar/gkt439
PLZF/ZBTB16, a glucocorticoid response gene in acute lymphoblastic leukemia, interferes with glucocorticoid-induced apoptosisThe Journal of Steroid Biochemistry and Molecular Biology 120:218–227.https://doi.org/10.1016/j.jsbmb.2010.04.019
Joseph S TakahashiReviewing Editor; Howard Hughes Medical Institute, University of Texas Southwestern Medical Center, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Blood transcriptome based biomarkers for human circadian brain pacemaker time" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Pål O Westermark (Reviewer #1); John B Hogenesch (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Laing and co-authors developed an unbiased circadian phase prediction method using partial least squares regression (PLSR) to model the relationship between whole-blood mRNA abundance profiles and their corresponding melatonin phase. Knowing body time is increasingly important for drug development and delivering medicines at optimal times for treatment. In this paper, the authors collected published datasets studying circadian transcription and plasma melatonin levels in human blood samples under different sleep conditions. Using these data, the authors compared published methods to predict body time with transcripts. The authors also developed their own partial least squares regression (PLSR) method to predict body time. Compared with the Molecular timetable and ZeitZeiger methods, PLSR better predicts circadian phase with blood transcripts.
In principle, the topic of the paper is important and the findings are novel and PLSR seems to be an advance; however, there was no statement regarding the availability of the PLSR code or information on what language it was written. Without this, it will be difficult to independently confirm the authors' work, or have the larger research community make use of the method. This reduced the enthusiasm for the paper among the reviewers.
1) Because the paper is long and so verbose, the main point of this paper is hard to understand. All of the main texts and figures should be more concise. Many figures are of technical nature and should perhaps be shown in the supplement. In the present manuscript version, the important messages drown in a wealth of material of more technical nature.
2) The SCN phase is not the same as the melatonin phase, the authors describe "predicting the circadian phase of the SCN" throughout the manuscript. As the authors described in Discussion, actual error in SCN phase measurement by melatonin is not known. The authors should use "melatonin phase" instead of "SCN phase". In addition, knowing the melatonin phase is not a panacea for delivering drugs at appropriate times for all diseases (for example, knowing the phase of the liver for a liver-specific disease may be more relevant than knowing the melatonin phase). The authors need to add discussion for this point.
3) The authors should make the code available (e.g. an R package, e.g. ZeitZeiger, or Python or whatever). In addition to making the code available, a more detailed description of PLSR should be included.
4) The authors did a good job of pointing to the primary data. However, they should also document the code/scripts used to run the Molecular timetable method, ZeitZeiger, and PLSR so that others can replicate this analysis.
5) What influence does the size of the training set have on PLSR's accuracy?
6) Can PLSR be used to rule out rhythmicity (ala Molecular timetable and ZeitZeiger)?
7) The Molecular timetable method has also been used with metabolite data. Does PLSR work with metabolite data as well?
9) The Agostinelli et al. paper is cited. BIO_CLOCK should be added to the algorithms benchmarked here.
10) This model improves upon previous models, but what advantage does this improved model have over previous models? Or, what advantage does this model have over just measuring melatonin? For example, what disease would benefit from having a model that has 20% error compared to 10% error? Of course, having a more accurate model is better, but it's not clear what problems or challenges that this new model will solve. Or what issues could be addressed that were not addressed by previous models?
11) The SD of the prediction error is on the order of 2-3 hrs. Is this likely to be enough for the method to be practically useful in e.g. clinical applications?
12) Authors should explain why they think only PER1 was identified in both molecular timetable methods (subsection “Molecular Timetable and Zeitzeiger methods”, second paragraph).
13) Authors should explain why they think those particular 4 genes were identified as circadian in all four conditions (subsection “Molecular Timetable and Zeitzeiger methods”, second paragraph). Because other hnRNPs can directly regulate AANAT, the rate-limiting step in melatonin biosynthesis, is it possible that hnRNPDL is also directly regulating AANAT (as opposed to, or in addition to, NRF as the authors suggest in their model)? And could this regulation occur via locally produced melatonin from peripheral blood mononuclear leucocytes?
14) In Figure 2, the cutoff correlation value of 0.3 is slightly arbitrary because there is so much variation in model performance with slightly different correlation values. For example, a cut-off correlation value of 0.4 performs almost as well as 0.3. Authors should discuss about this point.https://doi.org/10.7554/eLife.20214.034
- Emma E Laing
- Simon N Archer
- Derk-Jan Dijk
- Simon N Archer
- Derk-Jan Dijk
- Norman Poh
- Derk-Jan Dijk
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
We thank Drs Alpar Lazar, June Lo, Sibah Hasan, Emma Arbon, research nurses, research officers, program managers, and study physicians of the Surrey Clinical Research Centre for data acquisition and clinical support. We thank Malcolm von Schantz, Jonathan Johnston, Colin Smith, Giselda Bucca and John Groeger for helpful discussions and Benita Middleton for melatonin assays.
Human subjects: The protocols used to produce the data used in this work received a favourable opinion from the University of Surrey's Ethics committee and all participants provided written informed consent.
- Joseph S Takahashi, Reviewing Editor, Howard Hughes Medical Institute, University of Texas Southwestern Medical Center, United States
© 2017, Laing et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.