An in silico genome-wide screen for circadian clock strength in human samples

Abstract Motivation Years of time-series gene expression studies have built a strong understanding of clock-controlled pathways across species. However, comparatively little is known about how ‘non-clock’ pathways influence clock function. We need a strong understanding of clock-coupled pathways in human tissues to better appreciate the links between disease and clock function. Results We developed a new computational approach to explore candidate pathways coupled to the clock in human tissues. This method, termed LTM, is an in silico screen to infer genetic influences on circadian clock function. LTM uses natural variation in gene expression in human data and directly links gene expression variation to clock strength independent of longitudinal data. We applied LTM to three human skin and one melanoma datasets and found that the cell cycle is the top candidate clock-coupled pathway in healthy skin. In addition, we applied LTM to thousands of tumor samples from 11 cancer types in the TCGA database and found that extracellular matrix organization-related pathways are tightly associated with the clock strength in humans. Further analysis shows that clock strength in tumor samples is correlated with the proportion of cancer-associated fibroblasts and endothelial cells. Therefore, we show both the power of LTM in predicting clock-coupled pathways and classify factors associated with clock strength in human tissues. Availability and implementation LTM is available on GitHub (https://github.com/gangwug/LTMR) and figshare (https://figshare.com/articles/software/LTMR/21217604) to facilitate its use. Supplementary information Supplementary data are available at Bioinformatics online.


Introduction
The circadian clock network controls $24-h rhythms at the cell, tissue and organismal level. The core network is composed of circadian activators (e.g. BMAL1/CLOCK, RORs and DBP) that bind to E-box, RORE and D-box elements to activate transcription of the circadian repressors (e.g. PERs/CRYs, REV-ERBs, NFIL3, DECs and CHRONO), and other clock regulators (e.g. CSNK1D/E and FBXL3/21) (Takahashi, 2017). This feedback mechanism regulates the expression of hundreds to thousands of genes in the mammalian genome (Zhang et al., 2014). Importantly, increasing evidence suggests that the clock network and its time-keeping properties are influenced by other pathways. For example, hepatocyte-specific ablation of Smo, a key component of Hedgehog signaling, decreased the amplitude of clock genes (Bmal1, Clock and Nr1d2) in mouse liver (Marbach-Breitrü ck et al., 2019). We need a strong understanding of clock-coupled pathways in human tissues to better appreciate the links between disease and clock function.
Previously, efforts were made to explore clock-coupled pathways in human U2OS cells by genome-wide knockdown using small interfering RNAs (Zhang et al., 2009). U2OS cells are widely used circadian reporter model and show robust oscillations of luciferase expression with a period length of 24 h after synchronization. Hundreds of genes had strong circadian phenotypes, including period length changes and/or increases in amplitude in response to knockdown (Zhang et al., 2009). Pathway analysis revealed that insulin and Hedgehog signaling are coupled to the circadian clock in U2OS cells. However, an important limitation of the U2OS circadian reporter cells is that they do not model mammalian circadian systems at the tissue level.
Computational approaches can be used to explore clock-coupled pathways in human tissues. Indeed, the machine-learning-based CYCLOPS was developed to identify cycling genes in a human tissue without the need for collection time of day (Anafi et al., 2017). CYCLOPS was applied to multiple human tissues (Ruben et al., 2018;Wu et al., 2018). However, it has limitations, including optimal seed gene list selection (Ruben et al., 2018;Wu et al., 2020) and challenges in evaluating the quality of ordering without known time-stamped samples (Wu et al., 2018). These issues prevented our ability to order many publicly available datasets. Perhaps most importantly, and a general limitation of computational approaches in this field, is that they focus on identifying cycling genes. However, a gene that is cycling may not directly impact clock function. Conversely, a gene that is coupled to the clock network, but is not cycling may directly impact clock function.
We developed a new computational approach to explore candidate pathways coupled to the clock in human population samples without time information. This method, named LTM, is an in silico screen mimicking genome-wide knockdown experiments done in U2OS cells. As opposed to genetic manipulation, we use the natural variation in expression levels of genes in human population datasets to infer genes' influence on circadian clock function. The idea is that genes whose expression correlates with circadian clock function are clock-coupled candidates. The power of this approach is that it (i) directly links expression variability of a gene to clock function, (ii) does not depend on the longitudinal order of samples by timestamps or machine-learning inference and (iii) predicts pathways coupled to the clock network and/or the factors correlated with clock strength in a large-scale dataset.
We applied LTM to human population datasets to identify clock-coupled pathways in skin. We found that the cell cycle is the top candidate clock-coupled pathway in healthy skin, but not in skin tumors. In addition, we applied LTM to thousands of tumor samples in the TCGA database and found that extracellular matrix organization (ECM)-related pathways are tightly associated with the clock strength in human tumors. We released the LTM on GitHub and figshare to facilitate its use.

Datasets
All datasets used in this study are listed in Supplementary Table S1, including experimental design information, the accession number or download link and number of total samples for each dataset. The raw RSEM counts downloaded from FireBrowse (http://firebrowse. org/) were transformed to TPM values. The refGenome package (version 1.7.7) was used to parse the human gtf file (human genome version: GRCh38). For a gene with multiple transcript spliceisoforms, we selected the longest transcript as the representative gene length.

Main steps of the LTM strategy
We implemented the LTM strategy into the LTMR package. This R package is released on GitHub (https://github.com/gangwug/ LTMR). Major steps of performing LTM analysis include: Normalize the input expression matrix: This step can perform quantile normalization across samples, compute the expression profile for each gene if multiple splice-isoforms exist, remove genes with low expression value across samples and blunt extreme low or high outlier expression values. This step is not essential but is suggested. LTMR::LTMprep function is designed for this step.
Select genes for LTM analysis: This step separates samples into four quantile groups (Q1, Q2, Q3 and Q4) based on the expression level of each gene. It helps further filter out genes with low expression value in the Q1 group, and/or low fold change between Q4 and Q1. The mean expression value in each quantile group progressively increases from Q1 to Q4. LTMR::LTMcut function is designed for this optional step.
Compute the quantitative measures of circadian clock strength for each quantile group: LTMR::LTMheat function is designed for this essential step. Given a targeted number (e.g. 4) of quantile groups, LTMheat ranks samples by the expression level of each gene and separates them into four quantile groups (Q1-Q4) based on their expression level. LTMheat computes the expression correlation matrix of 17 clock and clock-associated genes (Wu et al., 2018) in each quantile group that serves as the query correlation matrix. The correlation matrix of 17 clock and clock-associated genes of mouse circadian atlas data (Zhang et al., 2014) serves as the reference correlation matrix (Wu et al., 2020). The phase relationships of clock genes are conserved from mouse to human (Wu et al., 2018). Mouse data reflects uniform genetics and environment, and, as such, serves as a reference for screening human samples by LTM. LTMheat applies Mantel's test to compute the similarity statistical value (Mantel's zstat value) between the query and reference correlation matrix. LTMheat also computes the nCV values (Wu et al., 2021) of 17 clock and clock-associated genes in each quantile group, and further calculates the mean nCV value of 17 clock and clock-associated genes. The number of quantiles depends on the total number of individual samples in the dataset. Increasing the number of quantiles is helpful for the LTM performance. However, LTMR::LTMheat is a time-consuming step. We also suggest balancing the computational time and LTM performance.
Generate the R(nCV) and R(zstat) for each gene: LTMR::LTMcook function is designed for this essential step. Given a targeted number (e.g. 4) of quantile groups, LTMcook generates a correlation value, R(nCV), between the mean expression value and the mean nCV value of 17 clock and clock-associated genes of four quantile groups. LTMcook also generates the correlation value, R(zstat), between the mean expression value and the Mantel's zstat value of four quantile groups. Both R(nCV) and R(zstat) vary between À1 and 1.
Calculate the LTMabs value for each gene: LTMR::LTMdish function is designed for this essential step. LTMdish averages the mean value of R(nCV) and R(zstat) as LTMpre, which varies between À1 and 1 for all genes. The advantages of averaging R(nCV) and R(zstat) include: (i) up-weighting genes with highly correlated R(nCV) and R(zstat); (ii) easier to rank genes with one integrated value (LTMpre); and (iii) prior knowledge of improved performance when averaging circadian parameters (e.g. phase and period) (Wu et al., 2016). For the same dataset, multiple LTMpre values will be generated when several targeted numbers of quantile groups are tested. The average of multiple LTMpre values for the same gene is calculated as LTMori, and its absolute value is defined as LTMabs for this gene. The LTMabs of all genes vary between 0 and 1.
Integrate the LTMori value of multiple datasets: When multiple datasets are available for the same tissue, this step can integrate multiple LTMori values for each gene. The integrated-LTMabs value for each gene is the absolute value of averaged multiple LTMori values. LTMR::LTMmeta function is designed for this optional step and output the integrated-LTMabs value.

LTM analysis of skin and cancer datasets
The key parameters of LTMR functions used to analyze 3 skin datasets and 11 cancer datasets are listed in Supplementary Table S2. Those cell-cycle regulation genes that encode cyclin and cyclindependent kinase were selected to map the LTMabs value. To compare these selected genes between healthy skin and skin tumors, the LTMabs value of each gene was normalized by the mean LTMabs value of 1000 rounds of randomly sampled genes in each dataset.

Enrichment analysis of LTM ranked genes in skin and cancer datasets
The integrated-LTMabs value of three skin datasets was used to rank the genes. Gene set enrichment analysis (GSEA) of all genes ranked by integrated-LTMabs value was performed using the fgsea package (https://bioconductor.org/packages/release/bioc/html/fgsea. html). The 'c2.cp.v7.1.symbols.gmt' was downloaded from MsigDB (https://www.gsea-msigdb.org/gsea/msigdb). From this file, we only used those REACTOME pathways with gene numbers between 10 and 1000 as the reference gene set of GSEA. Similar procedures were performed on the GSEA of all genes ranked by LTMabs of melanoma or integrated-LTMabs value of 11 cancer types.

Predicted cell fractions of tumor samples from 11 cancer types
The cell fractions of tumor samples were predicted with the immunedeconv package (Sturm et al., 2019), using EPIC method (Racle et al., 2017). The predicted fractions of cancer-associated fibroblast and endothelial cell in each tumor sample were used in this study.

In silico screen for genes associated with circadian clock strength in population data
Our hypothesis is that high amplitude molecular rhythms will lead to high amplitude physiological and behavioral rhythms. A strong clock is defined by high amplitude and conserved phases of clock genes. We evaluated clock strength at the network level through two measures from population samples ( Supplementary Fig. S1A-D): (i) an indirect measure of clock amplitude (nCV) (Wu et al., 2021) and (ii) a measure of conserved phase relationships between clock genes (Mantel's zstat) (Wu et al., 2018(Wu et al., , 2020. Using these two measures of clock strength, we developed a new method, LTM, to screen for genes associated with clock strength in population data. The principle behind LTM is that genes whose variation in expression correlates with variation in clock strength are potentially coupled to the molecular clock. We used human epidermis data (18 076 genes and 298 samples from 239 individual human donors) (Wu et al., 2018) to explain our strategy (Fig. 1). The schematic focuses on one example gene, CCNA2, to demonstrate the major steps of LTM. CCNA2 was chosen because it is one of the top gene by LTM score. First, we ranked all 298 samples by the expression level of CCNA2. We separated the 298 ranked samples into one of four quantile groups (Q1-Q4), and got 75 and 74 samples in Q1/Q3 and Q2/Q4, respectively. Next, we calculated the mean expression value of CCNA2 of each quantile group. We also quantified the clock strength per quantile group using two measures: the mean nCV of 17 clock and clockassociated genes, and the Mantel's zstat value. We calculated the correlation value between the mean expression of CCNA2 and two clock strength measures per quantile group, resulting in the terms R(nCV) and R(zstat). The mean value of R(nCV) and R(zstat) is defined by LTMpre. Then, we repeated this same process twice more. In details, the 298 samples ranked by the expression level of CCNA2 were separated into 7 and 10 quantile groups, respectively. This resulted in two Fig. 1. The schematic diagram of LTM strategy. The expression matrix contains 18 076 genes (one gene per row) and 298 samples (one sample per column). All 298 samples are ranked by the expression level of each gene (e.g. CCNA2) and separated into four quantile groups. There are 75 (Q1 and Q3) or 74 (Q2 and Q4) samples in each quantile group. The mean expression value of all samples in each quartile group is calculated, as indicated by barplot. Two measures of clock strength were calculated for each quantile group. The first is clock robustness that is indicated by the mean nCV of clock genes in Q1-Q4. The second is phase conservation of clock genes, which is indicated by Mantel's zstat value. The correlation matrix of clock genes in mouse circadian atlas and in each quantile group serves as the reference and query matrix, respectively. The zstat value of the Mantel test indicates the similarity between a query matrix and the reference matrix. Higher mean nCV value and Mantel's zstat value means a stronger clock. The correlation value between mean nCV or Mantel's zstat value and mean expression of CCNA2 (ROR1) from Q1 to Q4 is named as R(nCV) or R(zstat). The LTMpre is defined as the mean value of R(nCV) and R(zstat) of CCNA2 (ROR1). The LTMori is averaged LTMpre values for the same gene. The absolute value of LTMori is defined as the LTMabs additional LTMpre values for CCNA2. Finally, we computed an LTMori and LTMabs value for CCNA2. LTMori is the average of all three LTMpre values calculated when separating samples into 4, 7 or 10 quantile groups. A positive and negative LTMori value means that samples in a quantile group with higher (e.g. CCNA2) and lower expression (e.g. ROR1) of this gene has a stronger clock, respectively. LTMabs is the absolute value of LTMori. The LTMabs value makes it easier to explore clock-coupled pathways given that up-or downregulated genes in the same pathway may affect the circadian clock similarly. In sum, genes with higher LTMabs values represent a stronger correlation between expression level and clock strength. Next, we applied LTM to human population skin samples collected from healthy subjects in multiple sources.

Meta-analysis of three human population skin datasets using LTM
We observed inter-individual variation in skin clock strength in our prior study of 20 subjects sampled longitudinally over four timepoints (Wu et al., 2020). Therefore, we expected a range of clock strengths among single skin samples collected from hundreds of human subjects. On the other hand, expression levels of thousands of genes show variation among human skin samples (Kimball et al., 2018). LTM searches for strong correlations between gene expression levels and clock strength, which we applied to three skin datasets: (i) epidermis samples from 239 subjects (PG), (ii) sunexposed whole skin samples from 601 subjects (GTExSE) and (iii) non-sun-exposed whole skin samples from 479 subjects (GTExNSE). The distributions of LTMabs values were similar in three skin datasets ( Fig. 2A). To control for dataset-specific bias, we computed an integrated-LTMabs value to identify genes most strongly correlated with clock strength across datasets. For example, IPO13 ranks among the top integrated-LTMabs values indicating a strong linear correlation between its expression and mean nCV [R(nCV)>0.87 in all three datasets] or Mantel's zstat value [R(zstat)>0.86 in two datasets] (Fig. 2B). IPO13 encodes a member of the importin-beta family of nuclear transport proteins. Importinbeta genes have been implicated in circadian clock protein transport LTM is used to screen three groups of human population skin samples, including human epidermis samples from 239 subjects (PG), sun-exposed whole skin samples from 601 subjects (GTExSE) and not sun-exposed whole skin samples from 479 subjects (GTExNSE). The distribution of LTMabs for PG, GTExSE and GTExNSE is indicated by different color points. 'Integrated' represents LTMabs of a gene by integrating three datasets. A gene with a large LTMabs (e.g. IPO13) means its expression level is strongly correlated with clock strength in skin.

Cell-cycle gene expression correlates with clock strength in healthy human skin
We performed GSEA of genes ranked by integrated-LTMabs. Cellcycle regulation was among the top significantly enriched pathways ( Fig. 2C; ADJ.P < 0.01). Specifically, genes regulating G2/M-phase (e.g. CDK1, CCNB1 and CCNB2) have higher integrated-LTMabs values (i.e. their expression tightly correlates with clock strength) compared to genes regulating G1/S phase (e.g. CDK4/6, CCND1 and CCND2) (Supplementary Fig. S2A and B). CDK1 and cyclin B genes (CCNB1 and CCNB2) are critical components of M-phasepromoting factor (Hara et al., 2012;Nurse, 1990), which is the universal inducer of M-phase in eukaryotic cells. In sum, LTM found a strong correlation between mitotic entry and circadian clock strength in healthy skin samples.
Highly active cell proliferation is common in human tumors. Is mitotic entry coupled to the circadian clock in skin tumors? Using LTM, we screened 472 tumor samples of skin cutaneous melanoma (SKCM) from the TCGA database. Top ranked pathways include processive synthesis on the lagging strand, metabolism of watersoluble vitamins and cofactors and processive synthesis on the Cstrand of the telomere, but none of them were significantly enriched ( Fig. 2D; ADJ.P > 0.01). A previous study showed that DNA synthesis was required for circadian clock function in Neurospora (Liu et al., 2017). These results indicate that DNA synthesis rather than mitotic entry per se may be important to circadian clock function in skin tumors. We further compared the LTMabs value of cell-cycle regulators between healthy skin and skin tumors. Unlike what we observed in healthy skin samples, the LTMabs values of cell-cycle regulators were much lower in tumors ( Fig. 2E and Supplementary Fig. S2C), indicating a weak correlation between their expression levels and clock strength. This weak correlation may be caused by a dysregulation of cell cycle, circadian clock or coupling between them in tumors. Further experiments are required to clarify this. 3.4 Pan-cancer analysis using LTM Next, we asked which pathway is related to the clock strength in tumor samples. We extended LTM screening on tumor samples from other 10 cancer types. Cancer types included urothelial bladder carcinoma (BLCA), breast invasive carcinoma (BRCA), head-neck squamous cell carcinoma (HNSC), clear cell renal cell carcinoma (KIRC), lung adenocarcinoma (LUAD), low grade glioma (LGG), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), stomach adenocarcinoma (STAD) and thyroid carcinoma (THCA). We computed integrated-LTMabs values from total 11 cancer types (Fig. 3A; red point line). PER1 was among the top ranked genes, with lower expression correlating with weaker clock function in tumors [ Fig. 3B; R(nCV) and R(zstat) above 0.5 in 10 of 11 cancer types]. PER1 downregulation has been observed in many tumors, including breast, prostate, glioma and colorectal cancers (Savvidis and Koutsilieris, 2012), giving us confidence in this pancancer analysis. We performed GSEA on all genes ranked by their integrated-LTMabs values. Top significantly enriched pathways include ECM, degradation of ECM and ECM proteoglycans ( Fig. 3C; ADJ.P < 0.01). These results suggest that ECM-related pathways are associated with the clock strength in tumor samples across cancer types.

Clock strength in tumor samples is correlated with fractions of fibroblast and endothelial cells
Those LTM selected genes in ECM-related pathways are primarily collagen genes. We mapped integrated-LTMabs values to the human collagen gene family, which includes 33 genes in 8 subfamilies (Gelse et al., 2003). The collagen genes with high integrated-LTMabs values (LTMabs!0.5) are from the fibril-forming, basement membrane, microfibrillar and multiplexin subfamilies (Fig. 4A). In sum, LTM identified a strong correlation between collagen gene expression and circadian clock strength in tumors.
To understand the basis for this correlation, we asked how tumor samples with low collagen gene expression differ from tumor samples with high collagen gene expression. It is known that Type I collagen genes (e.g. COL1A1 and COL1A2) are highly expressed in human lung fibroblasts (Habermann et al., 2020). This encouraged us to test the correlation between COL1A2 expression and the proportion of fibroblasts across lung tumor samples. Using EPIC (Racle et al., 2017), we computed the fraction of cancer-associated fibroblasts in two lung cancer types, LUAD and LUSC. Indeed, we found a strong positive correlation between COL1A2 expression and fibroblast proportion in both LUAD and LUSC (Fig. 4B). We extended this analysis to test all the seven top collagen genes (LTMabs!0.5) among 11 cancer types. We found a consistent positive correlation between fibroblast proportions and expression of collagen genes belonging to fibril-forming (COL1A2) and microfibrillar (COL6A1, COL6A2 and COL6A3) subfamilies across nine cancer types, the exception being KIRC and LGG (Fig. 4C). Overall, this suggests that variation in clock strength between tumors is related to the proportion of fibroblasts in those samples.
On further inspection, we identified two Type-IV collagen genes (COL4A1 and COL4A2) whose expression levels are negatively correlated with fibroblast proportion in THCA and LGG. It is known that Type-IV collagen plays an important role in endothelial cell adhesion and migration (Herbst et al., 1988). Thus, we measured the correlation between COL4A1 expression and endothelial cell proportions in THCA and LGG samples, and found a strong positive correlation (Fig. 4D). The extended analysis of top seven collagen genes among 11 cancer types shows a consistent positive correlation between the fractions of endothelial cells and expression of basement membrane (COL4A1 and COL4A2) and multiplexin (COL15A1) collagen subfamily (Fig. 4E). In sum, variation in clock strength among tumor samples is related to the proportion of fibroblasts and endothelial cells. This makes sense considering non- What causes variation in the relative abundance of fibroblasts and endothelial cells between tumor samples? One possibility is technical bias whereby surgical resection captures more or less surrounding, non-tumor tissue. A second possibility is that tumors vary in cell composition. It is known that the immune cell composition of the tumor microenvironment varies among cancer types and between cancer patients (Ansell and Vonderheide, 2013). If true, perhaps an important question to ask is, do cancer cells surrounded by more fibroblast and endothelial cells have stronger clocks than other cancer cells? More importantly, do clock strength differences in cancer cells impact prognosis? We will leave these questions for later study.

Discussion
The circadian field has built a strong understanding of clockcontrolled pathways from years of time-series gene expression studies across species. However, comparatively little is known about how 'non-clock' pathways influence circadian clock function, especially in humans. We developed LTM as a window of insight into clock-coupled pathways in human population scale data.
LTM applied to thousands of healthy human skin samples revealed specific cell-cycle genes that are strongly correlated with clock strength. While it is known that the circadian clock regulates progression of the cell cycle (Kowalska et al., 2013;Matsuo et al., 2003), it is far less clear if and how the cell cycle regulates progression of the circadian clock (Droin et al., 2019;O'Neill and Hastings, 2008). We found that expression levels of cell-cycle genes (e.g. CDK1 and CCNB1) controlling the G2/M-phase are tightly correlated with clock strength in skin samples. What is the benefit of having a strong coupling between circadian clock and the cell cycle? It is known that the circadian clock regulates DNA repair (Gaddameedhi et al., 2011) and the cell cycle (Geyfman et al., 2012) in mouse skin. To minimize DNA damage induced by light exposure, a stronger circadian clock may direct more cells to divide at the right time. More cells dividing at the right time may add additional synchronization signals and further strengthen the skin clock. Therefore, the mutual regulation between the circadian clock and cell cycle may have evolutionary advantages in mammalian skin.
There is increasing evidence that circadian dysregulation contributes to cancer initiation and progression (Hadadi et al., 2020;Kettner et al., 2016;Papagiannakopoulos et al., 2016). This may be caused by less robust clock gene oscillations in tumor samples (Wu et al., 2021). A weak clock may release its temporal gating of the cell cycle (Kowalska et al., 2013;Matsuo et al., 2003), which may contribute to tumor initiation. The release of circadian gating on cell proliferation may in turn damage the coupling between the circadian clock and the cell cycle. Consistent with this idea, LTM on tumor samples found that expression levels of cell-cycle regulators are less correlated with clock strength in skin cancer than in healthy skin.
LTM revealed that the expression levels of collagen genes are tightly correlated with clock strength in tumors from 11 cancer types. Our subsequent analysis identified a strong positive correlation between expression levels of these collagen genes and the proportion of fibroblasts and endothelial cells in tumor samples. However, LTM cannot tell us anything about the reason why cell proportions vary among tumor samples. One possibility is that tumor microenvironment differs between cancer patients. For example, cancer cell composition may differ in patients. This may relate to the diversity of molecular clock strength between patients. Patients with a stronger molecular circadian clock may have tumors with higher fractions of fibroblasts and endothelial cells. An alternative explanation is that surgical resection obtained more or less tumor, unintentionally capturing margins. This provides an example that the LTM screened genes in the complex tissues (e.g. tumors) may pick other stronger factors (e.g. fractions of non-cancer cells) than clock-coupled pathways. Considering multiple factors may influence the clock strength in human population datasets, nonbiological factors (e.g. batch effects) may also affect the interpretation of LTM screening results. It is important to know biological context of the dataset before applying LTM.
LTM has other limitations. First, the genes identified by LTM are based on correlations. The strong correlation is not the direct evidence of causal effect. Therefore, the LTM predicted clockcoupled pathways need experimental validation and/or literature evidence. Second, LTM cannot screen genes correlated with the period and phase variation of the circadian clock. Third, ranking samples by the gene expression are strongly influenced by the circadian phase for strong cyclers, which makes it difficult to accurately quantify the clock strength when separating samples by circadian phase. So LTM may be biased to non-cycling genes. It may need to check the phase distribution bias if clock genes are enriched in the top LTM selected genes. To fully explore the power of LTM by the research community and fixing its limitations in the future, we release the LTM method into GitHub.
LTM is an in silico genome-wide screen of genes whose expression correlates with clock strength. This method takes advantage of the natural expression variation of genes in human population datasets. Compared to current rhythmic detection methods, the power of LTM is its ability to directly link expression variability of a gene to clock network properties, allowing for the detection of clockcoupled pathways without requiring known or predicted sampling times. An improved knowledge of clock-coupled pathways is important for the clinical application of circadian medicine. For example, a recent study shows that disruption of the clock gene BMAL1 impacts insulin sensitivity and liver disease (Jouffe et al., 2022). Another showed that targeting the Hedgehog signaling pathway affects the liver clock (Marbach-Breitrü ck et al., 2019). Therefore, when applying drugs targeting core clock genes or Hedgehog signaling pathway, the mutual influence needs to be considered. Outside of circadian medicine, there are well-known clinical risks of intervention when mutual regulation of pathways is poorly understood. For example, Vioxx (a COX2 inhibitor) is designed to inhibit pain and inflammation. However, Vioxx also inhibits the cardioprotection intermediates produced by COX2 and increased cardiovascular risk in patients (FitzGerald, 2004). This led to its withdrawal from the market, just 5 years after FDA approval and after many patients' deaths. LTM provides a chance of understanding clock coupling mechanisms at the tissue level, and the translation of circadian medicine in clinical application.