Evidence for widespread dysregulation of circadian clock progression in human cancer

Study design

The main goals of this study were to develop a co-expression-based signature of the circadian clock, to validate the signature in healthy human organs, and to use the signature to infer the extent of normal circadian clock progression in human cancer. We focused on transcriptome data because of its wide availability. We selected the datasets of circadian gene expression in mice (both for defining the reference signature and for comparing clock gene knockouts to wild-type) to represent multiple organs, light-dark regimens, and microarray platforms.

We selected the datasets from healthy human organs to include as many circadian studies as possible and to include a range of organs. The dataset from human skin consisted of samples taken at only three time-points for each subject (9:30 am, 2:30 pm, and 7:30 pm). Datasets from human blood consisted of multiple samples taken throughout the 24-h cycle for each subject. Datasets from human brain were based on postmortem tissue from multiple anatomical areas, and zeitgeber time for each sample was calculated using the respective donor’s date and time of death and geographic location. Datasets from cells cultured in vitro were based on time-courses following synchronization by dexamethasone, serum, or alternating temperature cycles. For GSE45642 (human brain), we only included samples from control subjects (i.e., we excluded subjects with major depressive disorder).

For datasets of human cancer, we analyzed data from the Cancer Genome Atlas (TCGA), as it contains samples from various cancer types. To complement the breadth of cancer types in TCGA data, we further included human liver cancer and lung cancer datasets from NCBI GEO. We analyzed human cancer datasets that had at least 30 non-tumor samples. When analyzing clock gene expression in human cancer, unless otherwise noted, we included all tumor and non-tumor samples, not just those from patients from whom both non-tumor and tumor samples were collected. For details of the datasets, all of which are publicly available, see Table S1.

Processing the gene expression data

For TCGA samples, we obtained the processed RNA-seq data (in units of transcripts per million, TPM, on a gene-level basis) and the corresponding metadata (cancer type, patient ID, etc.) from GSE62944 (Rahman et al., 2015). For E-MTAB-3428, we downloaded the RNA-seq read files from the European Nucleotide Archive, used Salmon to quantify transcript-level abundances in units of TPM (Patro et al., 2017), then used the mapping between Ensembl Transcript IDs and Entrez Gene IDs to calculate gene-level abundances.

For the remaining datasets, raw (in the case of Affymetrix) or processed microarray data were obtained from NCBI GEO and processed using MetaPredict, which maps probes to Entrez Gene IDs and performs intra-study normalization and log-transformation (Hughey & Butte, 2015). MetaPredict processes raw Affymetrix data using RMA and customCDFs (Irizarry et al., 2003; Dai et al., 2005). As in our previous study, we used ComBat to reduce batch effects between anatomical areas in human brain and between subjects in human blood (Johnson, Li & Rabinovic, 2007; Hughey & Butte, 2016). We also used ComBat to normalize expression between between subjects in the circadian human skin dataset (Fig. S8). ComBat is an empirical Bayes method to standardize the mean and variance of gene expression values between groups of samples.

Analyzing clock gene expression and co-expression

We first focused on the expression of 12 genes that are considered part of the core clock or are directly controlled by the clock and that show strong, consistently phased rhythms in multiple mouse organs (Zhang et al., 2014; Hughey, Hastie & Butte, 2016). We calculated times of peak expression (Fig. S2) and strengths of circadian rhythmicity of expression in wild-type and knockout mice (Fig. S15) using ZeitZeiger (Hughey, Hastie & Butte, 2016), with three knots for the periodic smoothing splines (Helwig & Ma, 2014). ZeitZeiger uses the time of day information for each sample.

We quantified the relationship between expression values of pairs of genes using the Spearman correlation (rho), which is rank-based and therefore invariant to monotonic transformations such as the logarithm and less sensitive to outliers than the Pearson correlation. Using the biweight midcorrelation (Song, Langfelder & Horvath, 2012), which is also robust to outliers, gave very similar results. All heatmaps of clock gene co-expression have the same mapping of correlation value to color, so they are directly visually comparable.

For a group of samples from a given dataset, we estimated the strength of clock gene co-expression (Fig. S7) as the difference between the 95th and 5th percentiles of the distribution of Spearman correlations between pairs of the 12 clock genes.

We calculated the reference Spearman correlation for each pair of clock genes (Table S2) using a fixed-effects meta-analysis (Hedges & Olkin, 1985) of the eight mouse datasets shown in Fig. S1. First, we applied the Fisher z-transformation (arctanh) to the correlations from each dataset. Then we calculated a weighted average of the transformed correlations, where the weight for dataset i was n_i − 3 (corresponding to the inverse variance of the transformed correlation), where n_i is the number of samples in dataset i. Finally, we applied the inverse transformation (tanh) to the weighted average.

To quantify the similarity in clock gene co-expression between the reference and a group of samples from another transcriptome dataset, we calculated the Euclidean distance between the respective Spearman correlation vectors, which contain all values in the strictly lower (or strictly upper) triangular part of the correlation matrix. We call this distance the clock correlation distance (CCD). A smaller CCD indicates greater similarity in clock gene co-expression. Although here we used Euclidean distance, other distance metrics could be used as well.

To estimate the probability that the observed CCD is less than would be expected by chance, we randomly selected 12 genes (from all genes measured in that dataset) 1,000 times, in order to estimate a null distribution of CCDs. We then calculated the exact one-sided p-value using the number of random samplings that gave a CCD less than or equal to the observed CCD and the method of Phipson and Smyth in the statmod R package (Phipson & Smyth, 2010). For the eight mouse datasets that form the reference, we calculated the CCD between each unique pair of datasets, essentially taking one dataset in each pair to be the reference.

To quantify the extent to which clock gene co-expression is perturbed in samples from one condition (e.g., tumor) relative to samples from another condition (e.g., non-tumor) in a given dataset, we calculated the CCD for each condition relative to the mouse reference, then calculated the difference between those CCDs. We call this difference the delta clock correlation distance (ΔCCD). In this study, a positive ΔCCD indicates that clock gene co-expression in tumor (or knockout) samples is less similar to the reference than is clock gene co-expression in non-tumor (or wild-type) samples.

To evaluate the probability that the ΔCCD is greater than would be expected by chance, we permuted the condition labels of the samples and recalculated the ΔCCD 1,000 times, always keeping the reference fixed. We then calculated the exact one-sided p-value using the number of permutations that gave a ΔCCD greater than or equal to the observed ΔCCD and the method of Phipson and Smyth in the statmod R package (Phipson & Smyth, 2010). Since we used the one-sided p-value, the alternative hypothesis was that non-tumor (or wild-type) is more similar to the reference than is tumor (or knockout).

To calculate the ΔCCD for individual tumor grades in TCGA data, we used the clinical metadata provided in GSE62944. We analyzed all combinations of TGCA cancer type and tumor grade that included at least 50 tumor samples. In each case, we calculated the ΔCCD using all non-tumor samples of the respective cancer type.

To compare ΔCCD and tumor purity, we used published consensus purity estimates for TCGA tumor samples (Aran, Sirota & Butte, 2015). The estimates are based on DNA methylation, somatic copy number variation, and the expression of immune genes and stromal genes (none of which are clock genes).

To quantify differential expression between tumor and non-tumor samples and between knockout and wild-type samples, we used limma and voom (Smyth, 2004; Law et al., 2014). To ensure a fair comparison between human and mouse data, we ignored time of day information in the mouse samples. We quantified the variation in expression of clock genes in each dataset and condition using the median absolute deviation (MAD), which is less sensitive to outliers than the standard deviation.

Analyzing circadian gene expression in human lung cancer

We obtained the set of rhythmic transcripts (identified by microarray probe ID) inferred by CYCLOPS to be rhythmic in healthy human lung (Anafi et al., 2017). This set contains transcripts whose abundance was well described by a sinusoidal function of CYCLOPS phase in samples from both the Laval and GRNG sites of GSE23546 (Bonferroni-corrected P < 0.05 and peak/trough ratio >2) and whose orthologs were rhythmic in mouse lung. We mapped the microarray probes to Entrez Gene IDs and calculated the acrophase for each gene as the circular mean of the provided acrophase for the corresponding probes. Phase values inferred by CYCLOPS are relative, so Anafi et al. adopted the convention of setting phase π(in radians) to the average acrophase of the PAR bZip transcription factors (DBP, HLF, and TEF) (Anafi et al., 2017).

Because some of the module preservation statistics calculated by NetRep use the expression matrix directly, we used ComBat (Johnson, Li & Rabinovic, 2007) to merge the expression data from the Laval and Groningen sites of GSE23546. We then followed WGCNA’s recommended procedure for identifying gene modules, i.e., groups of genes with correlated expression across samples (Langfelder & Horvath, 2008). Briefly, starting with the merged expression data, we calculated the Spearman correlation matrix for the 1,292 genes, then used the signed method with a soft thresholding power of 12 to calculate the adjacency matrix, then calculated the topological overlap dissimilarity matrix. We used hierarchical clustering of the dissimilarity matrix and adaptive branch pruning to define modules of genes, followed by a procedure to merge closely related modules. We then used the DAVID web application (version 6.8) to discover functional categories enriched in each module (Huang, Sherman & Lempicki, 2009).

We used NetRep (Ritchie et al., 2016) to calculate seven module preservation statistics between healthy human lung and the non-tumor and tumor samples from five datasets of human lung cancer (LUAD and LUSC from TCGA and GSE19188, GSE19804, and GSE32863 from NCBI GEO). We excluded GSE10072, because it included expression levels for only 876 of the 1,292 genes (the other five datasets included at least 1,199 of the 1,292 genes). When calculating module preservation, NetRep automatically removed any genes that were not measured in both healthy human lung and the group of samples to which it was being compared. For non-tumor and tumor samples from each dataset, we calculated the Spearman correlation matrix and adjacency matrix using the same procedure as for healthy human lung.

To evaluate the statistical significance of the difference in each module preservation statistic between non-tumor and tumor samples, we performed permutation testing and calculated one-sided p-values similarly to the procedure for ΔCCD, permuting the sample labels (non-tumor or tumor) in the test dataset and recalculating module preservation 1,000 times. Note this is different from the standard way WGCNA and NetRep perform permutations, which is to shuffle the gene labels.

Periodic smoothing splines of log₂ fold-change as a function of acrophase in healthy lung were fit using ZeitZeiger (Hughey, Hastie & Butte, 2016). P-values of non-zero amplitude of the spline fits were calculated using 1,000 permutations of the relationship between acrophase and log₂ fold-change and the method of Phipson and Smyth (Phipson & Smyth, 2010). For Phase Set Enrichment Analysis (Zhang et al., 2016), we used the canonical pathways gene sets from MSigDB (Liberzon et al., 2011) and looked for gene sets with a q-value ≤0.1 against a uniform distribution and vector-average value within 0.4 radians of either 0.07 π or 1.06 π (the mean phases of the peak and trough of the spline fits, respectively). No gene sets met the criteria for the former.

Consistent co-expression of clock genes in mice

The progression of the mammalian circadian clock is marked by daily rhythms in gene expression throughout the body (Zhang et al., 2014). We hypothesized that the relative phasing of different genes’ rhythms would result in those genes having a characteristic pattern of co-expression. Such a pattern could be used to infer the progression of the clock within a group of samples, even if the samples are not labeled with time of day, as long as they were taken from multiple parts of the circadian cycle.

Although most circadian gene expression is tissue-specific, the rhythms of a select group of genes show consistent relative phasing in multiple organs in mice and humans (Hughey & Butte, 2016). We focused on 12 genes from this group, each of which is part of the core clock or directly controlled by the clock, in order to define a general signature of circadian clock progression. For the rest of the paper, we refer to these 12 genes as “clock genes.”

As a first test of our hypothesis, we assembled eight publicly available datasets of genome-wide circadian gene expression from various healthy mouse organs (Oster et al., 2006; Hoogerwerf et al., 2008; Hughes et al., 2009; Negoro et al., 2012; Geyfman et al., 2012; Haspel et al., 2014; Zhang et al., 2014) (Table S1). The organs were collected from mice under both constant darkness and alternating light-dark cycles. For each dataset, we calculated the Spearman correlation between expression values (over all samples) of each pair of clock genes.

The pattern of co-expression was highly consistent across datasets (for 59 of 66 gene pairs, the sign of the Spearman correlation was the same in at least six of eight datasets; Fig. S1). Co-expression revealed two groups of genes, where the genes within a group tended to be positively correlated with each other and negatively correlated with genes in the other group. Genes in the first group (Arntl, Npas2, and Clock), which form the positive arm of the clock (Partch, Green & Takahashi, 2014), peaked in expression shortly before zeitgeber time 0 (ZT0, which corresponds to time of lights on or sunrise; Fig. S2). Genes in the second group (Cry2, Nr1d1, Nr1d2, Per1, Per2, Per3, Dbp, and Tef), which form the negative arms of the clock, peaked in expression near ZT10. Cry1, which was part of the first group in some datasets and the second group in others, tended to peak in expression near ZT18. These results indicate that the progression of the circadian clock in various mouse organs produces a consistent pattern of co-expression between clock genes.

To construct a single reference pattern of clock gene co-expression, we combined the eight datasets in a fixed-effects meta-analysis (Figs. 1A–1B and Materials and Methods). To compare clock gene co-expression between the reference and a group of samples from another dataset, we developed a metric we call the clock correlation distance (CCD), which corresponds to the Euclidean distance between the respective correlation vectors (Fig. 1A). A smaller CCD indicates a greater similarity in clock gene co-expression. To estimate the probability that the observed CCD is less than would be expected by chance, we estimated the distribution of CCD that would result from using 12 randomly selected genes instead of the 12 clock genes. As a positive control, we calculated the Euclidean distance between the correlation vectors corresponding to all pairs of the eight mouse datasets (P ≤ 0.001 by permutation test for each pair; Fig. 1C and Fig. S3).

We used the reference pattern and the CCD metric to verify that clock gene co-expression arises from the phasing of the genes’ rhythms relative to each other, not their phasing relative to time of day. Thus, the pattern is not affected by phase differences between groups of samples, such as those caused by daytime feeding in mice (Vollmers et al., 2009) (P ≤ 0.001 for CCD relative to reference; Fig. S4).

Most computational methods for quantifying circadian rhythmicity and inferring the status of the clock require that samples be acquired over the entire 24-h cycle. Because our approach does not attempt to infer oscillations, we wondered if it would be robust to partial coverage of the 24-h cycle. We therefore analyzed clock gene co-expression in three of the original eight datasets, in samples acquired during the first 8 h of the day (or subjective day) or the first 8 h of the night (or subjective night). In each dataset, clock gene co-expression was preserved in both daytime and nighttime samples (P ≤ 0.001; Fig. S5). Thus, our approach can detect the signature of clock progression in groups of samples without using time of day information, even if the samples’ coverage of the 24-h cycle is incomplete.

Validation of the co-expression pattern in human tissues

We next examined clock gene co-expression in nine publicly available datasets of circadian gene expression from human tissues: one from skin, two from brain, three from blood, and three from cells cultured in vitro (Hughes et al., 2009; Spörl et al., 2012; Möller-Levet et al., 2013; Li et al., 2013; Janich et al., 2013; Archer et al., 2014; Hoffmann et al., 2014; Arnardottir et al., 2014; Chen et al., 2016) (Table S1). Samples from these datasets are labeled with time of day (or in the case of cells in vitro, time since synchronization) and were taken from multiple time-points in the circadian cycle. Other than the two datasets of U2OS cells, which are derived from an osteosarcoma but show normal rhythms of the core clock genes (Vollmers, Panda & DiTacchio, 2008; Maier et al., 2009), all samples were from healthy tissue.

Although gene expression rhythms in mice and humans show different phasing relative to the light-dark cycle (Hughey & Butte, 2016), the pattern of clock gene co-expression in the human datasets significantly resembled the reference pattern from mice (P ≤ 0.001 for three of three in vitro datasets and four of six in vivo datasets; Figs. 1C–1D and Fig. S6). The CCDs of the datasets varied: skin and in vitro cells had the smallest (i.e., the strongest similarity to mice), brain was intermediate, and blood had the largest (Fig. 1C). Correlations between clock genes in the latter two organs were generally weaker than those in mice (Fig. S7), consistent with noisier circadian rhythms in those datasets (Hughey & Butte, 2016). The co-expression pattern in the human skin dataset, for which samples were acquired from each individual at three time-points (9:30 am, 2:30 pm, and 7:30 pm), was due to co-expression both across individuals at each time-point and across time-points within each individual (Fig. S8).

To confirm our findings in additional human organs, we analyzed eight transcriptome datasets from healthy human lung, liver, skin, and adipose tissue (Schadt et al., 2008; Innocenti et al., 2011; Bossé et al., 2012; Grundberg et al., 2012; Bonder et al., 2014) (Table S1). Samples from these datasets were not collected for the purpose of studying circadian rhythms. Thus, the samples are not labeled with time of day and may not cover the entire circadian cycle. Nonetheless, clock gene co-expression in these datasets was similar to that of mice (P ≤ 0.001 for all 8 datasets; Figs. 1C, 1E, and Fig. S9), with CCDs comparable to those of datasets designed to quantify circadian rhythms in humans in vivo. The generally weaker co-expression patterns and higher CCDs in human in vivo datasets compared to mouse datasets is likely due to data collected from humans having more sources of variation (including technical, environmental, and genetic variation). We conclude that our approach can detect the signature of circadian clock progression in publicly available transcriptome data from various human tissues, even in datasets not designed to quantify circadian rhythms.

Aberrant patterns of clock gene co-expression in human cancer

We then turned our attention to clock gene co-expression in human cancer. We started with RNA-seq data collected by The Cancer Genome Atlas (TCGA) and reprocessed using the Rsubread package (Rahman et al., 2015). TCGA samples are from surgical resections performed prior to neoadjuvant treatment. We analyzed TCGA data from the 12 cancer types that included at least 30 samples from adjacent non-tumor tissue (Table S1). To study human hepatocellular carcinoma and non-small cell lung cancer in more depth, we also analyzed eight datasets from the NCBI Gene Expression Omnibus (GEO), each of which had tumor and adjacent non-tumor samples (Landi et al., 2008; Hou et al., 2010; Lu et al., 2010; Roessler et al., 2010; Lamb et al., 2011; Selamat et al., 2012; Lim et al., 2013; Villa et al., 2016) (Table S1). For all data from TCGA and GEO, the times of day of surgery are not available. Presumably, most to all of the samples were acquired during surgical working hours (6:00 am to 6:00 pm).

To quantify the extent to which clock gene co-expression is perturbed in tumor samples compared to non-tumor samples in a given dataset, we calculated the difference between the CCDs (relative to the mouse reference) for each condition. We call this difference the delta clock correlation distance (ΔCCD), where a ΔCCD >0 indicates that clock gene co-expression in tumor samples is less similar to the reference than that of non-tumor samples (Fig. 2A). We continued to use the reference from mouse data, instead of one based on human data, in order to maintain a strict separation of training and validation data and because in the mouse data, the effects of other sources of variation and confounding are minimized. Thus, the mouse data offered the clearest window into the clock gene co-expression that results from circadian clock progression in vivo. To evaluate the probability that the ΔCCD is greater than would be expected by chance, we permuted the condition labels of the samples and recalculated the ΔCCD 1,000 times.

Clock gene co-expression in non-tumor samples was similar to that observed in solid healthy human tissues (P ≤ 0.005 for CCD relative to the mouse reference in 14 of 20 datasets; Figs. 2B–2C and Fig. S10). In contrast, co-expression in tumor samples often lacked a significant similarity to the expected pattern (P > 0.05 for CCD in 10 of 20 datasets; Figs. 2B–2C and Fig. S10). Consequently, tumor CCDs were significantly higher than non-tumor CCDs (P = 1.4⋅10⁻⁸ by paired t-test) and each human cancer dataset had a ΔCCD >0 (P ≤ 0.005 in 15 of 20 datasets; Fig. 2D and Fig. S11A). Among the three TCGA cancer types with the lowest ΔCCD, prostate adenocarcinoma had a relatively high non-tumor CCD, whereas renal clear cell carcinoma and thyroid carcinoma each had a relatively low tumor CCD (Fig. S11B). Interestingly, we observed no clear trend between ΔCCD and histological tumor grade (Fig. S12). Overall, these results suggest that circadian clock progression is perturbed in a range of human cancers.

Tumors are a complex mixture of cancer cells and various non-cancerous cell types. The proportion of cancer cells in a tumor sample is called the tumor purity and is an important factor to consider in genomic analyses of bulk tumors (Aran, Sirota & Butte, 2015). We therefore examined the relationship between ΔCCD and tumor purity in the TCGA data. With the exception of thyroid carcinoma and prostate adenocarcinoma, ΔCCD and median tumor purity across TCGA cancer types were positively correlated (Fig. S13; Spearman correlation = 0.67, P = 0.059 by exact test). These findings suggest that at least in some tumors, clock progression is perturbed more strongly in cancer cells than in non-cancerous cells.

Distinct patterns of clock gene expression in human cancer and clock gene knockouts in mice

To better understand the nature of the dysregulation of clock progression, in particular to determine if the pattern of clock gene co-expression in human tumors could be explained by inactivation of core clock genes, we compared clock gene expression in human cancer and in clock gene knockouts in mice. We assembled seven datasets of circadian gene expression, measured in various tissues, that included samples from wild-type mice and from mice in which at least one core clock gene was knocked out, either in the entire animal or in a specific cell type (Vollmers et al., 2009; Cho et al., 2012; Nikolaeva et al., 2012; Paschos et al., 2012; Dyar et al., 2014; Young et al., 2014; Dudek et al., 2016) (Table S1). For each dataset, we calculated clock gene co-expression in wild-type and knockout samples (Fig. S14).

The two datasets with the highest ΔCCD (>50% higher than any ΔCCD we observed in human cancer) were those in which the knockout mice lacked not one, but two components of the clock (Cry1 and Cry2 in GSE13093; Nr1d1 and Nr1d2 in GSE34018; Fig. 3A). The ΔCCDs for the other five datasets were similar to or somewhat lower than the ΔCCDs we observed in human cancer. Due to the smaller sample sizes compared to the human cancer datasets, the ΔCCDs for four of seven datasets of clock gene knockouts were not significantly greater than zero (P > 0.05).

To further compare clock gene expression in human cancer and clock gene knockouts, we calculated differential expression of the clock genes between non-tumor and tumor samples and between wild-type and knockout samples (Fig. 3B). Differential expression in the knockouts was largely consistent with current understanding of the core clock. For example, knockout of Arntl (the primary transcriptional activator) tended to cause reduced expression of Nr1d1, Nr1d2, Per1, Per2, Per3, Dbp, and Tef, and increased or unchanged expression of the other clock genes. In the double knockout of Cry1 and Cry2 (two negative regulators of the clock), this pattern was reversed. Interestingly, neither of these patterns of differential expression was apparent in human cancer.

In the clock gene knockouts, rhythmic expression of most clock genes was reduced or lost (calculated using each sample’s time of day; Fig. S15). Without time of day information in the human cancer datasets, it was not possible to directly quantify rhythms of clock genes in non-tumor and tumor samples. However, we reasoned that a proxy for rhythmicity could be the magnitude of variation in expression. Therefore, for each human cancer dataset, we calculated the median absolute deviation (MAD) in expression of the clock genes in non-tumor and tumor samples. We then compared the log₂ ratios of MAD between tumor and non-tumor samples to the log₂ ratios of MAD between knockout and wild-type samples from the clock gene knockout data (Fig. 3C). As expected, samples from clock gene knockouts showed widespread reductions in MAD compared to samples from wild-type mice. In contrast, human tumor samples tended to show similar or somewhat higher MAD compared to non-tumor samples. The distinct patterns of differential expression and differential variability between tumor and non-tumor samples, compared to knockout and wild-type samples, imply that dysregulation of clock progression in human cancer is not due solely to the inactivation of one or more clock genes.

Dysregulation of broader circadian gene expression in human lung cancer

Finally, we explored how dysregulation of clock progression might affect the circadian transcriptome. We obtained a set of 1,292 genes that were inferred by the CYCLOPS method to be rhythmic in a dataset from healthy human lung (Anafi et al., 2017). By applying WGCNA (Langfelder & Horvath, 2008) to the same dataset, we grouped the 1,292 genes into five modules according to their co-expression (Fig. 4A). Based on DAVID (Huang, Sherman & Lempicki, 2009) and consistent with the analysis of Anafi et al. (2017), the five modules were enriched for genes involved in various biological processes, including angiogenesis, phosphoprotein signaling, and alternative splicing (Table S3). We then used NetRep (Ritchie et al., 2016) to determine the extent to which each module was preserved in non-tumor and tumor samples from five datasets of human lung cancer (two from TCGA and three from NCBI GEO). Notably, all five datasets had ≤65 non-tumor samples and the three datasets from NCBI GEO had ≤91 tumor samples, too few for CYCLOPS to produce a reliable ordering. Based on this analysis, four of five modules (comprising 1,256 of 1,292 rhythmic genes) were significantly more strongly preserved in non-tumor compared to tumor samples (Fig. 4B, P ≤ 0.001 for at least two of seven preservation statistics in at least four of five datasets).

We also quantified differential expression of the 1,292 normally rhythmic genes between non-tumor and tumor samples in the five lung cancer datasets. We then examined the relationship between each gene’s log₂ fold-change (tumor vs. non-tumor) and its phase of peak expression (acrophase) in healthy lung as inferred by CYCLOPS. Remarkably, we observed a consistent pattern, in which genes with an acrophase near πtended to show the strongest decrease in expression in tumor samples (Fig. 4C, P ≤ 0.001 of non-zero amplitude by permutation testing for each dataset). Using Phase Set Enrichment Analysis (Zhang et al., 2016), we found that genes with an acrophase near π were particularly enriched for involvement in G protein-coupled receptor signaling (Table S4). Taken together, these findings suggest that perturbed clock progression in human lung cancer is accompanied by systematic changes in broader circadian gene expression.