Identification of the oscillating transcriptome

Averaged FPKM values at each time were used for cycling transcripts identification. Lowly expressed transcripts were removed by calculating the background expression in each participant using the average expression of a panel of 62 genes known not to be expressed in peripheral blood cells (Table S2). Temporal transcriptomes were participant to linear detrend prior to identification of oscillations by either the eigenvalue/pencil or RAIN methods. For the eigenvalue/pencil method [^{10, 12}], a maximum of four oscillations were identified for each gene. Criterion for circadian genes were: period between 20h to 25h for first and second participants and 24h to 30h for the third participant, decay rate between 0.8 and 1.2; for ~12h genes: period between 9.6h to 13.6h for the second and third participants and 10h to 14h for the first participant, decay rate between 0.8 and 1.2; for ~8h genes: period between 6h to 8h for the first participant and 7h to 9h for the second participant, decay rate between 0.8 and 1.2; for ~16h genes; period between 14h to 18h for the third participant. The relative amplitude was calculated by dividing the amplitude by the mean expression value for each gene. To determine FDR, we used a permutation-based method that randomly shuffles the time label of gene expression data and subjected each permutation dataset to the eigenvalue/pencil method applied with the same criterion [²⁹]. These permutation tests were run 5,000 times, and FDR was estimated by taking the ratio between the mean number of rhythmic profiles identified in the permutated samples (false positives) and the number of rhythmic profiles identified in the original data. Analyses were performed in MatlabR2017A. RAIN analysis was performed as previously described in Bioconductor (3.4) (http://www.bioconductor.org/packages/release/bioc/html/rain.html) with either 48h continuous data or 24h data with biological duplicates as input [¹¹]. We included genes exhibiting a period between 10h and 14h with a p value less than 0.05 as having ~12h expression in all three participants. FDR was calculated using the Benjamini-Hochberg procedure. Heat maps were generated with Gene Cluster 3.0 and TreeView 3.0 alpha 3.0 using Z score normalized values.

For meta-analysis, we used Fisher’s method, which combines extreme value probabilities from each test, commonly known as “p-values”, into one test statistic (X2) using the formula

where pi is the p-value for the ith hypothesis test.

For RAIN analysis on temporal IR events and splicing gene expression, raw data was subjected to polynomial detrend (n=2) before RAIN analysis.