Challenges and opportunities for statistical power and biomarker identification arising from rhythmic variation in proteomics

Investigating time of day variation in different human biofluids and tissues has been a growing focus of research over the past decade¹. This variation can be driven by endogenous circadian rhythms, which are influenced by central and/or peripheral clocks, as well as by exogenous diurnal rhythms resulting from behavioural and environmental cycles. In turn, dysregulation or disruption of these chronobiological variations has been associated to a wide range of diseases and conditions^2,3,4,5,6. These rhythms are well described in terms of metabolomics and transcriptomics^7,8,9,10; however, chronobiological variations in proteomic expression are less well mapped^11,12,13,14. This is in part due to the sheer range of proteins and the relatively high costs of untargeted proteomic experiments; the precise number is unknown but a count of around 10,000 distinct proteins has been suggested for peripheral blood, ranging in concentration over at least 9 orders of magnitude¹⁵. Whilst proteomic studies have shown rhythms in biological pathways involved in immune function, metabolism, and cancer, amongst others¹⁴, there are still substantial gaps in current knowledge about which proteins and pathways exhibit diurnal or endogenous rhythmicity.

The rhythmicity of protein expression is relevant beyond analysis of circadian rhythm disorders, as temporal variation increases the potential for study design error and reduces statistical power. On the first issue, if a protein has a rhythmic component, this creates the potential for confounding if studies are designed without taking rhythmicity into account. This might occur, for example, through selection / sampling bias, such as in a case / control study if all cases were to be measured in a hospital setting during a ‘morning round’, and controls were to be measured by convenience sampling, perhaps later in the working day. This bias would increase the risks of Type I errors, i.e. false positive identification of biomarkers as differentiating between cases and controls, with the risks of bias proportionate to the rhythm amplitude. Rhythmicity also increases the risk of Type II errors, by reducing statistical power. This reduction in statistical power for a given study follows naturally from the increase in variance in the features (here proteins) being measured. These challenges to study design can, however, also offer opportunities. By controlling for chronobiological variation through—for example—controlled time-of-day sampling, variance can be reduced and statistical power improved. Such study design steps are likely to be much more cost-effective than simply increasing the number of participants.

In this work, we provide a case study of the significance of rhythmicity for a small set of high-abundance proteins in human serum, by measuring protein concentration at a two-hourly resolution across 30 sequential hours and identifying rhythmic features alongside their acrophase (peak) values and amplitudes. We use these data to provide a real-world illustration of the challenges to study design that rhythmicity presents, in terms of increased n to maintain statistical power (or the opportunity to reduce n and maintain statistical power, if rhythmicity is controlled for). Finally, we set out four suggested mitigation steps against both Type I and Type II errors driven by chronobiology in proteomic analyses. These include incorporating time of sampling in study design; giving consideration to chronobiological issues in statistical power calculations; reporting time of sampling for cases and controls (with p values) as metadata; and reporting any literature-identified rhythmicity in any biomarkers of interest, to help provide context for readers.

Experimental model and subject details

The samples analysed in this work were selected from a previous study as described in Isherwood et al¹⁶. Briefly, 24 healthy male volunteers were recruited by the Surrey Clinical Research Facility (CRF); the study was given a favourable ethical opinion from the University of Surrey ethics committee and was conducted in accordance with the Declaration of Helsinki; all participants gave written informed consent. Participants’ blood was sampled during constant routine (60 blood draws: every 30 minutes for 30 sequential hours). Inclusion criteria included being male; aged between 18 and 35 years; BMI between 18 and 30 kg/m² (inclusive) at screening; and habitual (at least 5 days per week) hours in bed per night between 7 and 9 hours which included a bedtime between 22:00–01:00 and waketime between 06:00–09:00. An ESS (Epworth Sleepiness Scale) score <9 was treated as indicative of normal range daytime sleepiness, an HÖ (Horne-Östberg; diurnal preference) score between 30 and 70 was taken as the normal range and a PSQI (Pittsburgh Sleep Quality Index) score <5 was considered indicative of satisfactory sleep quality and sleep. Shift workers involved in night work within the past six months and volunteers that had travelled across more than two time zones within a month before the study were excluded. Additional exclusion criteria were, smokers or nicotine products users within the 6 months prior to screening, and any regular use of medication known to influence circadian rhythm. A 10-day pre-laboratory routine, a previously validated standard for our human chronobiology experiments, was employed^17,18. The full details of both the pre-laboratory protocol and the laboratory protocol are set out in Isherwood et al.¹⁶, including allowed sleep windows, motion and light measurement via Actiwatches [Cambridge Neurotechnology, Cambridge, UK] and diet diaries, supplemented by continuous glucose monitors (CGMs) [Freestyle Libre 2, Abbott Laboratories Limited]. Compliance with the pre-laboratory routine was monitored by CGM confirmed compliance to the meal eating times. Upon admission in the CRF in the afternoon of Day 0, the participants were assessed to confirm compliance to the pre-laboratory routine and a review of medication, an alcohol breath test and a urine sample for analysis of cotinine and drugs of abuse were performed. Participants were supervised throughout the laboratory sessions by medical/clinical research staff. Laboratory environmental conditions and meals are again set out in Isherwood et al.¹⁶. During the laboratory period 12 participants consumed hourly small meals throughout the waking period and 12 consumed two large daily meals. Following the diet regimes, all participants underwent a constant routine protocol with 30–60 min blood sampling. From the dataset of 12 individuals on the small-meal regime 11 were selected at random. The samples from the constant routine were taken forward for this proteomics analysis.

Sample preparation

Human sera (50 µL) were aliquoted into fresh Eppendorf tubes prior to denaturing with 5 µL of 0.1% (w/v) RapiGest™ (Waters Corporation, Milford, MA) in 50 mM ammonium bicarbonate and incubated at 80 °C for 45 min. Following incubation, 100 mM DTT (3 µL) was added and incubated for a further 30 mins at 60 °C to reduce the proteins, before being alkylated with 200 mM iodoacetamide (3 µL) at room temperature for 30 min. Trypsin 1:50 (w/w) (Gold Mass Spectrometry grade, Promega, Madison, WI, USA) was added to each sample for proteolytic digestion and left incubating overnight at 37 °C. TFA was added to a final concentration of 0.5% (v/v) to hydrolyse the RapiGest and heated for a further 45 min at 37 °C, before centrifuging for 25 min at 18,000 g. The supernatant was collected and 3 µL aliquoted for LC-MS analysis. Aliquoted samples were diluted 1:250 (v/v) with 750 µL of 0.1% FA (v/v) and 19 µL of MassPREP™ Digestion Standard Mix 1 (Waters Corporation, Milford, MA) was added as an internal reference.

LC-MS analyses

Peptides resulting from the tryptic digests were analyzed using the Evosep One EV-1000 (Odense, Denmark) coupled to a SYNAPT™ XS mass spectrometer (Waters Corp., Wilmslow, UK). Samples were loaded onto the Evosep tips as per the manufacturer’s instructions. Peptides were separated using the Evosep 60 SPD method, configured with a EV-1064 column. MS data were collected with the SYNAPT XS mass spectrometer, operated in positive electrospray ionisation (ESI) mode with a nominal resolution of 25,000 FWHM (V optics). The capillary voltage was 3.2 kV, the cone voltage was 35 V and the source temperature was set at 100 °C. Data were acquired over 50–2000 Da mass range with a scan time of 0.5 s. All mass spectral data were acquired in continuum mode using HDMSE to obtain fragmentation data simultaneously¹⁹. Function one (low energy) data were collected using a constant trap and transfer energy of 6 eV, while the second (high energy) function consisted of a transfer collision energy ramp of 19 to 45 eV. For mass accuracy, [Glu1]-fibrinopeptide (m/z = 785.8426) was acquired as lock mass at a concentration of 100 fmol/µL (in 50:50 CH₃CN/H₂O, 0.1% formic acid). Lock mass scans were collected every 60 s and averaged over 3 scans to perform mass correction. The time-of-flight was externally calibrated over the acquisition mass range (50–2000 Da) before analysis with a NaCsI mixture. These data were collected using MassLynx v 4.1 software (Waters Corp., Wilmslow, UK) in a randomized order with three technical replicates acquired per sample. Lock mass consisting of [Glu1]-Fibrinopeptide was delivered to the reference sprayer of the MS source using the M-Class Auxillary Solvent Manager with a flow rate of 1 µL/min.

Data analysis

The times of sampling were mapped to each participant’s dim light melatonin onset (DLMO) time by a previously described method²⁰. Progenesis QI for Proteomics (Nonlinear Dynamics, Newcastle upon Tyne, UK) was used to process all LC-MS data. Retention time alignment, peak picking and normalization were conducted to produce peak intensities for retention time and m/z data pairs. Data were searched against reviewed entries of a Homo sapiens UniProt database (20,435 reviewed entries, release 2024) to provide protein identifications with a false discovery rate (FDR) of 1%. A decoy database was generated as previously described²¹, allowing for protein/peptide identification rates to be determined. Peptide and fragment ion tolerances were determined automatically, and searches allowed for one missed cleavage site. Carbamidomethyl of cysteines was applied as a fixed modification, whilst oxidation of methionines and deamidation of asparagine/glutamine were set as variable modifications. Following this process, a dataset in the form of an array comprising protein concentrations across features by time by participant was generated. These concentrations were then normalised by standard scaling, i.e., dividing by the per-participant standard deviation to express each concentration as a participant z-score. This generated an array of scaled protein concentrations with dimensions of participant n by DLMO time t by protein identifier p.

Cosinor analysis for rhythm detection was performed using CosinorPy (version 3.0) in the Python programming language (version 3.9.18), using the Spyder IDE (version 5.4.3)^22,23. All code used in the preparation of this manuscript used open-access libraries without modification. A single-component cosinor model was used in this work; using CosinorPy’s population-mean module, a cosinor model was fitted to each individual. The response of the whole population was described and analysed as a mean of all individual cosinor models. CosinorPy includes functionality for multiple components, and also for the inclusion of non-linear factors, but given the small sample size employed here, complex models were rejected due to the risk of overfitting. Proteins that were identified as rhythmic with p value < 0.05 and amplitude > 0.1 were collated and processed for pathway upregulation / downregulation using the STRING online platform^24,25, in order to identify clusters of proteins that were significantly altered (enrichment p value < 0.05). Since Benjamini-Hochberg assumes independence of variables, and will be too conservative when proteins are correlated within pathways²⁶, FDR correction was applied at the pathway level, rather than at the individual protein level.

To assess the impact of rhythmicity on biomarker identification and statistical power, we start with the assumption that the variance of a biomarker measured in a cohort controlled for chronobiological variation is ({sigma }_{{base}}^{2}). To this base variation, an uncontrolled chronobiological variation represented by a cosine function would add additional variance ({sigma }_{{cosine}}^{2}), plus the covariance ({sigma }_{{bc}}^{2}), as shown in Eq. 1. If there is no relationship between base variation and variation due to chronobiology, the covariance term is zero. Furthermore, the chronobiological variance could be expressed solely in terms of the amplitude A of the cosine function, yielding Eq. 2. For biomarkers with multiple components, the total variance will be the sum of the base variance (excluding any rhythmic component) plus the sum of ({A}^{2}/2) for each additional cosinor component (by the variance sum law)²⁷.

Because variance increases when adding uncorrelated cosinor components, the number of participants n required to maintain statistical power would be increased for given critical values Z and effect size d, as illustrated in Eq. 3. The critical values used here were based on α = 0.05 and β = 0.20.

For each rhythmic protein identified here, the increase in n required to maintain statistical power is reported. Additionally, the impact on power (Z_β and therefore β) is derived from Eq. 4).

Cohort analysis

11 participants were randomly selected from the small meal cohort of 12 from the original study (1 was excluded due to constraints on LC-MS runtime); 1 additional participant was excluded from subsequent analysis due to failed LC-MS injections. The key characteristics of the 10 study participants analysed in this work are summarised in Table 1. All recruited participants were male.

Rhythm analysis

Proteins (n = 281) were identified across the triplicate injections. All proteins with greater than 10% missing values were excluded from further analysis, leaving 73 proteins for rhythm analysis. Of these, 13 proteins were identified as including a 24-hour or 12-hour rhythm based on two criteria: a p-value when fitted to a single component cosinor of <0.05 and an amplitude > 0.1. Of the 13 proteins, 9 showed a 12-hour rhythm and 6 showed a 24-hour rhythm, with an overlap of 2 proteins. The individual proteins identified as rhythmic are shown in Table 2 with their acrophases and amplitudes. These amplitudes are shown in terms of the intra-individual z-scores, i.e. relative to the standard deviation exhibited by each individual.

Both the 24-hour and 12-hour sets of proteins were subjected to separate pathway analysis using STRING. For the 12-hour set of proteins, one cluster of proteins was identified (Fig. 1A) with a p value of < 0.001 for the overall level of enrichment. The proteins within this cluster were investigated for functional enrichment. Different databases may produce different results, but here showed considerable overlap; enrichment according to GO Biological Processes was significant for fibrinolysis (p = 0.020), negative regulation of blood coagulation (p = 0.004) and blood coagulation (p = 0.020); whilst according to KEGG Pathways, functional enrichment was significant for the complement and coagulation cascades (p < 0.001). The full lists of pathways are shown in Table S1, Supplementary Material. Comparison of acrophases for the proteins identified also showed clustering for PLG, KNG1 and AZGP1 around 5 hours after DLMO ( = 03.77 (dec. h, clock time) (Fig. 1B).

A Node and edge relationships of proteins identified as having statistically significant 12-hour rhythms (p value < 0.05) and forming a functionally enriched cluster. Network nodes represent proteins and edges represent protein-protein interactions, line thickness indicates the strength of data support (B) polar plot of acrophases of proteins shown in Fig. 2A, amplitude is shown in terms of z-score along the radius, DLMO time is shown around the perimeter (DLMO 0 = clock time 22:46 h:min). All analyses based on a cohort of n = 10.

Full size image

For the 24-hour set of rhythmic proteins, one cluster of proteins was identified (Fig. 2A) with a p value of p < 0.0001 for the overall level of enrichment. The proteins within this cluster were investigated for functional enrichment in the same way as for the 12-hour set of proteins. Enrichment according to GO Biological Processes was significant for pathways including plasminogen activation (p = 0.004), negative regulation of triglyceride metabolism (p = 0.004); while for KEGG pathways, functional enrichment was again significant for the complement and coagulation cascades (p < 0.001), for cholesterol metabolism (p = 0.016) and platelet activation (p = 0.049). The full lists of pathways containing rhythmic proteins are shown in Table S2, Supplementary Material (proteins not found to be rhythmic in this work are also listed for completeness, in Table S3, Supplementary Material). Comparison of peak timings showed that complement and fibrinolysis-related proteins all had an acrophase between a DLMO time of -2 h and 0 hours (clock time 20.77–22.77 h)(Fig. 2B).

A Node and edge relationships of proteins identified as having statistically significant 24-hour rhythms (p value < 0.05) and forming a functionally enriched cluster. Network nodes represent proteins and edges represent protein-protein interactions, line thickness indicates the strength of data support (B) polar plot of acrophases of proteins shown in Fig. 3A, amplitude is shown in terms of z-score along the radius, DLMO time is shown around the perimeter (DLMO 0 = clock time 22:46 h:min). All analyses based on a cohort of n = 10.

Full size image

Finally, the potential impact of rhythmicity on biomarker identification and pathway contribution is summarised in Table 3. This table shows for each protein the amplitude identified in this work, and the increase in n required to maintain statistical power if rhythm is not controlled, or conversely the decrease in n achievable for a given Type II error rate if rhythm is controlled.

Power and implications for statistical errors

In addition to the specific proteins analysed here, the data also provide a case study for calculating the impact of biological rhythmicity on statistical power. Both the Type II error rate and the number of participants required increase as the amplitude of the rhythmic component increases. The overall relationship is shown in Fig. 3. For any given amplitude, there will be a required increase in n needed to maintain statistical power (Fig. 3A), or an increased Type II error rate if n is not increased (Fig. 3B). Conversely, controlling for rhythmicity allows for the maintenance of statistical power without increasing n. Both functions (for n and for the Type II error rate) vary with the square of the amplitude, so both charts show the same exponential relationship.

A Relationship between required n and amplitude of rhythm to achieve a 20% Type II error rate (B) Type II error rate with fixed n (number of participants) and an increasing independent cosinor rhythmic component

Full size image

Temporal variation is a key aspect of physiology. Daily variation in the transcriptome and metabolome of human tissues and fluids have been described^28,29,30, but rhythms of human proteomic data are poorly understood. Here we provide an analysis of the influence of temporal variation in the human serum proteome. Overall, 13 out of the 73 proteins meeting QC thresholds were identified as rhythmic, 6 being statistically significant at 24 hours, and 9 being statistically significant at 12 hours (including 2 rhythmic at both). The identification of two groups of proteins with different cycle durations is consistent with recent work on 12 hour ultradian rhythms, including in humans³¹, in the transcriptome. The generation of such rhythms in cell culture is poorly understood^32,33, but has recently been suggested to be linked to degradation, rather than transcriptional activation³⁴. Ultradian rhythms have been hypothesised to be relevant to the maintenance of metabolic homeostasis, but the regulatory and synchronising mechanisms of these 12 hour rhythms remain unknown^35,36.

The overall proportion of proteins identified as rhythmic was 18%. This is broadly consistent with a recent larger proteomics study which found that 15% of analysed proteins exhibited significant daily rhythmicity¹³, and is also concordant with metabolomics studies suggesting a range of 15% to 20% of features exhibiting circadian rhythmicity in constant routine conditions (in nonconstant / entrained conditions this can be higher)⁹. Rhythmicity of many of the biological functions captured here is already well-described, with the immunoglobulins identified here having previously been described as having a circadian component^37,38,39. As well as immunoglobulins, this work also highlights a number of protein clusters and pathways as rhythmic, such as the coagulation and the complement cascade as well as apolipoproteins. Several of the individual proteins here have previously been identified as rhythmic, in particular the apolipoproteins APOA2, APOC3 and APOE^13,40. Within the complement and coagulation cascades, FIBA, FIBB and KNG1 have also previously been reported as rhythmic⁴⁰. From the total of 13 rhythmic features analysed here, 4 proteins – AZGP1, ITIH2, PLG and CFH – are to our knowledge identified for the first time as rhythmic in this work.

Circadian variation of fibrinolytic and complement activity in blood is a well-described phenomenon⁴¹, and has for example been reported to contribute to increased risk of cardiovascular events in the morning⁴². The proteins identified here have also been closely associated with conditions known to dysregulate circadian rhythms, such as Alzheimer’s disease. For example, CFH within the complement cascade has been identified as a biomarker of Alzheimer’s^43,44, as have the fibrinolysis pathway proteins FGA and FGB^45,46, and indeed as has AZGP1⁴⁷. Three apolipoproteins were also identified as rhythmic; the genes governing the expression of many of this family of proteins have previously been demonstrated to show rhythmic expression⁴⁸ and these results are also consistent with lipid metabolism more broadly having a circadian component^49,47; Apolipoproteins have also been directly linked to Alzheimer’s disease; especially with regard to APOE4^50,51. It should be noted that in some cases, where rhythmic proteins are identified as biomarkers of a condition, there is a risk that they are in fact biomarkers of generic circadian dysregulation, and are not specific to the disease in question. This is a known issue in multivariable analyses such as proteomics or metabolomics, especially when machine learning algorithms are trained on idealised case-control datasets of disease versus healthy participants⁵².

In terms of limitations, the work presented here is a small case study measuring only a limited number of proteins. It should be noted that recent analyses of proteomic rhythmicity have identified a much wider range of proteins with a circadian component, due to utilising methods with greater depth of proteomic sampling, such as aptamer-based multiplexed platforms^13,40. Nonetheless, the statistical impact on biomarker or pathway identification shown in this case study is relevant for any rhythmic protein, transcript, metabolite, or other chronobiological biomarker of interest. Therefore, we can view chronobiology as representing a potential ‘omics confounder across a variety of fields, as well as those directly related to chronobiological research, such as chronotherapy and systems biology modelling of circadian rhythms. As presented here, controlling for rhythmic variation allows for a reduction of sample size n whilst maintaining statistical power, whilst conversely failing to control for rhythmicity ceterus paribus increases Type II error rates. Therefore, rhythmicity represents both a meaningful opportunity (for improved statistical power at relatively low cost) in biomarker research, as well as a cost (in the form of increased error rates and reduced reproducibility), when not controlled for.

We see the following suggested measures as helpful in aiding the reader’s understanding of biomarker research in the context of chronobiology, reducing the risks of Type I and Type II errors, and improving reproducibility.

1. 1.

   The ideal is to control for rhythmicity as part of the study design, especially when the biomarkers being investigated have a known rhythmic component or in untargeted studies that may have a rhythmic component. At the simplest level, time of day of sampling could be controlled; a better approach is to adjust for circadian phase with the gold-standard being the alignment to DLMO, albeit this will be beyond the scope of the vast majority of projects. Alternatively, single point-in-time methods could be used to calculate and adjust to a DLMO time⁵³.

2. 2.

   If rhythmicity cannot be controlled for during sampling, statistical power calculations during the study design phase should acknowledge the impact of rhythmicity, and the necessary increase in n or the reduction in β.

3. 3.

   We also recommend that time of day of sampling be considered as part of a study’s metadata, and should be reported, alongside for example, the p-value for any difference between case and controls’ time of sampling. This would also be of benefit in biobank data, as researchers would have the option of selecting time-of-day matched data to control for rhythmicity.

4. 4.

   Finally, where rhythmicity in any identified biomarkers has previously been reported, this should be noted, especially in studies reviewing biomarkers, pathways, or conditions, which are known to be influenced by the human timing system.

In conclusion, the rhythmic protein expressions shown in this work demonstrate the potential for confoundment; biological rhythmicity therefore presents a study design problem for case-control experiments in proteomics and other ‘omics platforms. Whilst aligning rhythms to an individual’s circadian phase using DLMO to remove the confounder is likely to be too costly for the majority of studies, adoption of the recommendations suggested in this work would mitigate against the risks of Type I and Type II errors and improve reproducibility of biomarker identification.