Circadian Regulation Patterns With Distinct Immune Landscapes in Gliomas Aid in the Development of a Risk Model to Predict Prognosis and Therapeutic Response

The RNA sequencing data of TCGA-LGG (Lower Grade Glioma; WHO grade II–III), TCGA-GBM (glioblastoma multiforme; WHO grade IV), and GTEx-brain in transcripts per kilobase million (TPM) format uniformly processed by TOIL were downloaded from UCSC xena (https://xenabrowser.net/datapage/↗). The corresponding clinicopathological data were obtained from the cBioPortal website (http://www.cbioportal.org/↗). Patients without survival information were excluded in a further study. The genomic mutation data (including somatic mutation and copy number variation) of TCGA-LGG and TCGA-GBM were downloaded from the Genomic Data Commons (https://portal.gdc.cancer.gov/↗) using the “TCGAbiolinks” R package (21). Somatic mutation data were analyzed using “maftools” R package, and significant amplifications or deletions of copy number were detected using GISTIC 2.0. To obtain the CGGA validation set, the RNA sequencing data and related clinicopathological data were downloaded from the CGGA website (https://www.cgga.org.cn↗). Two independent datasets, including one with 325 individuals and another with 693 individuals, were acquired. Batch effects from non-biological technical biases were corrected using the “ComBat” algorithm of the “sva” R package (22).

Six CR gene sets were concluded from the Molecular Signatures Database (MSigDB) (23), KEGG pathways (24), REACTOME database (25), GO biological processes (GO: 0007623 CIRCADIAN RHYTHM), WikiPathways (26), and review literatures (27–29). The gene annotated in not less than two CR gene sets was defined as a CR gene for further integrated analysis and modeling. A total of 91 acknowledged CR genes were curated and analyzed to identify distinct CR patterns (Table S1 and Figure S1). The CR gene widely described in reviews is defined as a core CR gene. Unsupervised clustering analysis was applied to identify CR patterns based on the expression of 91 CR genes from TCGA or CGGA and classify patients for further analysis. A consensus clustering algorithm was performed using the “ConsensuClusterPlus” R package and was repeated 1,000 times for guaranteeing the stability of classification (30).

To investigate the difference in biological process between CR patterns, we performed GSVA enrichment analysis using the “GSVA” R package (31). The gene sets of “h.all.v7.4.symbols” and “c2.cp.kegg.v7.4.symbols” were downloaded from the MSigDB database to run GSVA analysis (http://www.gsea-msigdb.org/gsea/downloads.jsp↗). The GSVA scores among different circadian patterns were determined using the “limma” R package. The threshold was set as the adjusted p-value < 0.05 and the absolute (log₂ Fold change, FC) > 0.1. The “clusterProfiler” R package was used to perform functional annotation for CR-related genes, and the top 20 GO annotations were shown (32).

As we did in our previous study, the “ESTIMATE” algorithm was used to calculate Tumor Purity, ESTIMATE Score, Immune Score, and Stromal Score (33, 34). Furthermore, we used the ssGSEA (Single-sample Gene Set Enrichment Analysis) algorithm, CIBERSORT algorithm (35), and xCELL algorithm (36) to quantify the relative abundance of each type of cells infiltrating the glioma TME. The gene sets of “TME-associated signatures” were downloaded and analyzed using the “IOBR” R package (37). The ssGSEA method was chosen in the process of signature score evaluation. TIP (Tracking Tumor Immunophenotype) is a meta-server that systematically integrates two existing third-party methods “ssGSEA” and “CIBERSORT” for tracking, analyzing, and visualizing the status of anti-cancer immunity and the proportion of tumor-infiltrating immune cells across a seven-step cancer-immunity cycle using RNA-seq or microarray data (38). The correlations between the CRscore and the steps of the cancer-immunity cycle were analyzed using the “ggcor” R package.

The construction of CRscore was performed as follows: firstly, the differentially expressed genes (DEGs) associated with the CR cluster phenotype were determined using the “limma” R package. Specifically, gene expression data were normalized by voom and then fed to lmFit and eBayes functions to calculate the differentially expressed statistics. The significance filtering criteria of DEGs were set as an adjusted p-value < 0.05 and |log₂(FC)| ≥1.5. Then, the prognostic analysis was performed for each DEGs using a univariate Cox regression model. A total of 719 genes with significant prognosis were extracted for further analysis using the “randomForestSRC” R package. We also used the “randomSurvivalForest” algorithm to rank the importance of 719 genes, and the gene with a relative importance >0.2 was selected to construct the CR signature (39). We then curated the expression profile of the final 50 determined genes to perform PCA, and principal components 1 and 2 were extracted and served as the signature score. We then adopted a formula similar to previous studies (40), CRscore = Σ(PC1_i + PC2_i), where i is the expression of 50 prognostic-related CR cluster DEGs.

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (41) and Immune Cell Abundance Identifier (ImmuCellAI) algorithm (42) were used to predict response to ICB therapy. We also used immunotherapeutic cohorts with complete clinical information to predict patients’ response to ICB therapy. In the GSE78220↗ cohort, the data of patients with metastatic melanoma treated with pembrolizumab, an anti-PD-1 antibody, were downloaded from Gene Expression Omnibus (GEO), and the FPKM data of gene expression profiles were converted to the more comparable TPM value among samples. We also used the “Prophetic” R package to predict the chemotherapy response of each sample based on Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/↗). The correlation between estimated IC₅₀ and CRscore was established using Spearman correlation analysis with thresholds as |Spearman R| > 0.5 and p-value < 0.05.

For comparison between two groups, statistical significance for normally distributed variables was estimated by Student’s t-tests, and nonnormally distributed variables were analyzed by the Wilcoxon rank-sum test. For comparison among three or more groups, one-way ANOVA and Kruskal–Wallis test were used for normally or nonnormally distributed variables, respectively. Fisher’s exact test was performed for categorical data. Correlation coefficients were computed by Spearman and distance correlation analyses. The “surv-cutpoint” function from the “survival” R package was applied to stratify samples into CR-high and -low groups. Kaplan–Meier curves and the Log-Rank test were adopted to assess whether there were differences in overall survival among groups. Univariate and multivariate Cox regression analyses were utilized to evaluate the independent prognostic value of the CRscore regarding OS, which were revealed by forest maps. We next performed multivariate Cox regression to establish a nomogram; the survival predictive accuracy of prognostic models was assessed based on the calibration plots. All statistical p-values were two-sided, with p < 0.05 considered as statistically significant. All statistical tests were performed in R statistical software (v4.05, R Core Team, R Foundation for Statistical Computing, Vienna, Austria).

To investigate whether there were specific expression patterns of CR genes in gliomas, we determined several differentially expressed CR genes between glioma and normal samples using the “limma” R package (Figure 1A). The genes involved in negative transcription regulation (ID3, ID2, and KLF10), cell metabolism (NACLU and AHCY), and cell survival (MAGED1 and TIMELESS) were upregulated in tumors when |log₂FC| ≥ 1 and adj.p ≥ 0.05 were set as thresholds. Furthermore, the comprehensive landscape of interaction among CR genes and feature groups to which they belong was depicted in the Figure 1B network, and we noticed that most of the CR genes were positively correlated with each other. To investigate whether genomic variation contributed to abnormal CR gene expression in gliomas, we described mutations and copy number variations (CNVs) of CR genes, respectively. In the top 10 mutations of CR genes, only PTEN exerted 11% mutant frequency, suggesting that genetic mutations of CR genes might not be the main cause of rhythm perturbations (Figure 1C). The vast majority were mutations in PETN, a common tumor suppressor gene, which could be included here because of its expression governed by circadian rhythms. Besides, the CNV frequency showed that CLOCK, the positive arm of TTFL, gained copy number, while PER3, a part of negative arm of TTFL, had a frequency of CNV deletion (Figure 1D). Thus, we have come to know that CNV partly conduced to disturbed rhythms of tumors. The location of CR genes with CNVs on chromosomes is shown in Figure 1E. To explore whether the core regulators of circadian TTFLs had unique expression patterns in gliomas, we analyzed the mRNA levels of core CR genes among normal, LGG, and GBM samples (Figure 1F). It was worth noting that major factors (CRYs and PERs) of negative arms were upregulated in LGG compared to normal, but with relatively low expression in GBM. Given all this, the disruption of circadian rhythms due to the high heterogeneity of genomic variations and transcription profile covertly contributed to tumorigenesis and progression of gliomas.

Given the aberrant expression of certain circadian regulators in gliomas, we wondered whether there were distinctively different CR patterns, on the whole, among individual tumors. The unsupervised hierarchical cluster analysis was used to uncover circadian patterns, and the TCGA cohort population was separated into two distant clusters, termed CR cluster-A and -B, by the “ConsensusClusterPlus” R package (Figures 2A, S2A–C). Additionally, similar results obtained in the CGGA cohort verified our findings above (Figures S3A–C). Moreover, principal component analysis confirmed that CR genes could distinguish two CR clusters perfectly, in both TCGA and CGGA cohorts (Figure 2B). Obvious discrepancy between two CR clusters was embodied in histopathology features and patients’ prognoses. Patients in CR cluster-A were characterized by lower grades and harmless molecular genetic features, which coincided with particularly prominent survival advantage (Figures 2C, S2D, E). The same results were obtained in LGG and GBM groups, respectively (Figures S2D, S3D). We also added the unsupervised clustering results in LGG and GBM, respectively. Unsupervised clustering could also divide LGG patients into two groups, and there was a significant difference in prognosis between the two groups (Figures S2F–I). As for core CR genes, the distribution of their transcriptome levels with remarkable differences between two clusters corresponded with the previous clustering (Figure 2D). ARNTL (BMAL1) and CLOCK together constitute the positive arm of TTFLs, but they took advantages of expression level in different patterns. Similarly, the components of the negative arm were expressed preferentially.

Next, to explore the global functions in two CR patterns beyond the individual and single gene expression, we performed GSVA enrichment analyses, which showed the activation states of Hallmark pathways (MSigDB) in TCGA (Figure 2E) and CGGA (Figure 2F) cohorts. Hedgehog Signaling activation, Pancreas beta cells, and Kras Signaling Downregulation were observed in CR cluster-A, while immune response-related pathways consisting of “Interferon Response”, “Interleukin Signaling”, and “Tnfα Signaling via Nfκb” and pathways associated with malignant biological behaviors including “Epithelial-Mesenchymal Transition” and “Angiogenesis” were activated in CR cluster-B. KEGG pathway enrichment analyses in TCGA (Figure S3F) and CGGA (Figure S3G) cohorts also supported the activated immune reaction-related pathways. From above, we identified two internally gathered circadian patterns and functional enrichment analyses recognized CR cluster-B closely related to the immune response.

To describe immune status in different CR patterns, we utilized four different methods (ESTIMATE, ssGSEA, CIBERSORT, and xCELL) for evaluating immune cell infiltration (Figure 3A). ESTIMATE algorithm showed that CR cluster-B exhibited high Stromal and Immune Scores and low Tumor Purity, which represented a significantly increased immune cell infiltration and generally indicated a relatively hot tumor immune microenvironment. However, in gliomas, components of immune system enriched in CR cluster-B were macrophages (ssGSEA, CIBERSORT, and xCELL), dendritic cells (ssGSEA and xCELL), neutrophils (ssGSEA, CIBERSORT, and xCELL), endothelial cells (xCELL), and Th2 cells (ssGSEA and xCELL), generally considered as immunosuppressive components in tumors. Additionally, CR cluster-B was mainly characterized by relatively high expressions of human leukocyte antigens (HLA) and immune checkpoint molecules (Figure 3B). On the one hand, class I HLA, expressed on the surface of almost all nucleated cells, could increase the possibility of presenting more immunogenic antigens and of benefiting from ICB. High immune checkpoint expression also implied effective anti-immune checkpoint therapy. On the other hand, in line with the characteristics of immune cell infiltration, class II HLA molecules, mainly expressed on the surface of antigen-presenting cells (dendritic cells, B cells and macrophages), and many inhibitory molecules (IL-4, IL-10, and VEGF) were unregulated in CR cluster-B. We then examined the specific correlation between each TME infiltration cell type and each core CR gene using Spearman’s correlation analyses (Figure 3C). The correlation focused on the relationship between immunosuppressive cells (macrophages, neutrophils, and T helper cells) and the core TTFL components (ARNTL, BHLHE40/41, CRY2, and PER2/3). Furthermore, CR clusters could be distinguished by different TME relevant signatures and TME score (Figure 3D). The above results revealed that immunosuppressive components enriched in CR cluster-B suggested an immunosuppressive microenvironment and coincided with poor prognoses of patients.

To further validate CR regulation patterns and investigate the potential biological behavior of CR patterns, we determined 719 genes, which were DEGs between CR cluster-A and -B and also related to prognosis of glioma patients using the “limma” R package and univariate Cox regression analysis (Figure 4A). Then, we performed unsupervised clustering analyses based on the obtained CR-related DEGs. Consistent with the clustering grouping of previous CR patterns, the unsupervised clustering algorithm also revealed two distinct CR-related signatures termed CR gene cluster-A and -B, respectively (Figure 4B, Figures S4A, B). This confirmed that two distinct CR patterns did exist in gliomas. Survival analysis for CR gene clusters in the TCGA cohort using Kaplan–Meier curves suggested survival advantages in CR gene cluster-A no matter in all glioma samples (Figure 4C) or in TCGA-LGG and TCGA-GBM, separately (Figure S4C). Next, the “clusterProfiler” R package was used to perform GO enrichment analysis for Gene Set 1 that were upregulated in gene cluster-B (Figure 4D), and Gene Set 2 upregulated in gene cluster-A (Figure 4E). GO’s annotations for Gene Set 1 fell into several categories including antigen presentation, inflammation, and immune responses, while the annotations for Gene Set 2 focused on neurotransmitter transmission and other major duties of CR genes controlling the brain. These analyses confirmed that gene cluster-B was significantly associated with immune-relevant signatures, whereas gene cluster-A was associated with regulations of synaptic transmission. Besides, the expression of core CR genes (Figure S4D), the infiltration of immune cells (Figure S4E), and Hallmark (Figure S4F) as well as KEGG (Figure S4G) pathways in two CR gene clusters were consistent with the previous clustering results. The results of secondary clustering verified that there were indeed two CR-based patterns in gliomas.

Considering the intratumoral heterogeneity and complexity of circadian disruption, based on 50 prognostic-related CR cluster DEGs, we developed a score termed CRscore using the PCA algorithm. The general distribution of CRscore, CR clusters, gene clusters, and molecular pathological parameters in the TCGA cohort were described in the heatmap (Figure 5A). Low CRscore crowded in CR cluster-B and gene cluster-B (Figure 5B). Considering that the immune microenvironment of LGG and GBM might differ, we explored the immunological relevance of CRscore separately in LGG and GBM. The correlation was supported by TME-related pathways and immune cell infiltration analyses and even more pronounced in LGG (Figures 5C, D). Next, we explored the indicating effect of CRscore in the anti-tumor immune process (Figure 5E) and discovered that CRscore was likely to be associated with the recruitment of immune cells. Some steps including priming and activation of immune cells; recruiting Th1, Th22 cells, and macrophages; and tumor killing effect were positively correlated with CRscore in both LGG and GBM. Peculiarly, in steps significantly correlated with the CRscore, a low score in GBM is related to recruitment of immunosuppressive cells (neutrophils, eosinophils, basophils, Tregs, and MDSCs). A low score in LGG indicated enhanced recruitment of anti-tumor immune cells (CD8+ cells, Th1 cells, and DCs).

After having identified the CRscore as an intrinsic pattern closely linked to the immune process, we sought to determine whether the CRscore could accurately predict outcomes of glioma patients. Firstly, the overlap among tissue types, CR clusters, gene clusters, and CR groups was described in Figure 6A. The 596 patients in the TCGA cohort were assigned to two groups (CR-high and -low) based on CRscores using the cutoff value (0.4) obtained with the “survminer” R package. Survival analysis indicated that the CR-high group had prolonged survival time in the TCGA cohort as well as in LGG or GBM (Figure S5A), which was further validated in the CGGA cohort (Figures 6B, S5B). The area under the ROC curves (AUC) at 5 years (Figure 6C) and time-dependent AUC (Figure 6D) of the CRscore compared with other histological or molecular indicators in the TCGA and CGGA cohort, respectively, verified the predictive efficiency of the CRscore. When the CRscore was evaluated as a categorical variable (high or low CRscore) in the Cox regression model, the CRscore was determined to be an independent and robust prognostic factor. Univariate and multivariate analyses of clinicopathological characteristics and CRscore in the TCGA and CGGA cohort are shown in Figures 6E, F, respectively. Furthermore, we established a prognostic nomogram to predict 3-, 5-, and 10-year overall survival based on the stepwise Cox regression model. WHO grade, CRscore, and age were included in the prediction model (Figure 6G). The C-index of the nomogram was 0.864 (95% CI, 0.840–0.888). Nomogram prediction and actual observation in the TCGA cohort reached an excellent agreement at the 3-, 5-, and 10-year survival probability after calibration (Figure 6H). Net decision curve and the net reduction analyses demonstrated the superiority of this nomogram in predicting prognosis (Figures 6I, J).

Given that the aim of exploring the value of the CR pattern on treatment and the response to ICB are closely related to somatic mutation, we decided to analyze differences in somatic mutation between CR-high and -low groups in the TCGA cohort using the “maftools” R package. As shown in Figure 7A, the CR-low group presented more extensive tumor mutation burden (TMB), on the whole, than the high-score group. Differences in CRscores between wild-type and mutant groups of each gene inversely proved the negative correlation between mutation and score (Figures 7B, C). The mutation occurrence varied between two groups. In detail, genetic mutations of IDH1, TP53, ATRX, and CIC mainly appeared in the CR-high group, while TTN, PTEN, and EGFR appeared in the CR-low group (Figure 7D).

Newly identified predictors, such as TIDE (Figure 8A) and ImmuCellA (Figure 8B), are widely used to evaluate the immune response. Our analysis revealed that a lower CRscore was not only associated with a poorer response to ICB but also more prone to immune escape. In anti-PD-1 immunotherapy cohort (GSE78220↗), a survival benefit trend was observed in patients with high CRscore (Figures 8C–E). To further understand the effects of the CRscore on drug response, we assessed the association between the CRscore and the response to drugs in cancer cell lines. Using the Spearman correlation analysis, we identified twenty-two drugs targeting various pathways (Figure 8G) significantly correlated between CRscores in the Genomics of Drug Sensitivity in Cancer (GDSC) database (Figure 8F). Among them, twenty agents with lower IC₅₀ in the CR-low group might be beneficial to low CRscore patients, while two agents might be beneficial to patients with high scores (Figure 8H). Together, these results implied that CR pattern played a crucial role in mediating the immune response and was also correlated with drug sensitivity. Thus, the CRscore might be a potential biomarker for establishing appropriate treatment strategies.

Publicly available datasets were analyzed in this study. These data can be found here: TCGA (The Cancer Genome Atlas); CGGA ( Chinese Glioma Genome Atlas).