Identification of a Novel Nomogram to Predict Progression Based on the Circadian Clock and Insights Into the Tumor Immune Microenvironment in Prostate Cancer

All somatic mutation data, raw counts of RNA-sequencing data (level 3) and corresponding clinical information of PCA were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/↗), and the method of acquisition and application complied with the guidelines and policies. The RNA-seq data in the format of fregments per kilobase per million (FPKM) were converted to the format of transcripts per million reads (TPM), and log2 conversion was performed (24). The prognostic data were integrated into the analyzed dataset (25). Genes with a false discovery rate (FDR)-adjusted p value < 0.001 and an absolute value of log2 (fold change) > 1 were considered differentially expressed genes (DEGs). Genes associated with progression-free interval (PFI) were defined as an adjusted p value < 0.05 through Cox regression analysis. We used PFI as a prognostic factor, which was defined as the period during and after treatment in which a participant was living with a disease that did not worsen. Usually, it is the period from date of diagnosis until 1) locoregional or systemic recurrence, 2) second malignancy, or 3) death from any cause; late deaths not related to cancer or its treatment are excluded (25). The CIC-related genes were obtained from Genecards (https://www.genecards.org/↗) (26). Somatic mutations included nonsilence mutations, nonsynonymous mutations, deletion mutations, frameshift mutations, insertion mutations and so on. Tumor mutation burden (TMB) is defined by the number of somatic mutations (in addition to synonymous mutations and intron mutations) per genome area (38 Mb) for target sequencing (27). Microsatellite instability (MSI) is a pattern of hypermutation that occurs at genomic microsatellites and is caused by defects in the mismatch repair system (28). MSI data were obtained from R package “TCGAbiolinks”. The specific classification and description of markers in 24 immune cells were performed by the previous article (29).

Validation of candidate genes at the protein level was performed through the HPA database (https://www.proteinatlas.org/↗) (30–32). The genetic alteration and mutation location of the ten genes were derived from cBioPortal (http://www.cbioportal.org/↗; TCGA, PanCancer Atlas; RNA Seq V2 RSEM) (33, 34). The GSCALite database was used to further analyze mutation information, immune infiltration, and gene set variation analysis (GSVA) of enrolled genes (http://bioinfo.life.hust.edu.cn/web/GSCALite/↗) (35). In addition, we explored the Genomics of Drug Sensitivity in Cancer (GDSC) and the Cancer Therapeutics Response Portal (CTRP) drug sensitivity in pan cancer through GSCALite (35).

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) analyses of each candidate gene and gene set were conducted to explore possible biological functions and signaling pathways using the package R “clusterProfiler”. GO analysis included biological process (BP), cell composition (CC) and molecular function (MF) (P<0.05 was statistically significant).

All analyses were conducted with R version 3.6.3 (https://www.r-project.org/↗) and its suitable packages. Principal component analysis was performed by the R package “Factoextra” to nonlinear dimensionality reduction and cluster analysis. To make reliable immune infiltration estimations, we utilized the R package “immunedeconv”, including xCell, MCP-counter, and CIBERSORT (36–42). The ssGSEA algorithm included in the R package “GSVA” was also used to assess the immune infiltration level (43). The immune score was imputed by using the “estimate” package. SIGLEC15, TIGIT, CD274, HAVCR2, PDCD1, CTLA4, LAG3, and PDCD1LG2 were selected as immune checkpoint-relevant transcripts, and the expression values of these eight genes were extracted (44–47). Lasso Cox regression modeling was conducted by using the “glmnet” package. A predictive nomogram was further constructed based on Cox regression analysis. The chi-square test was used to assess differences between categorical variables, and t tests or paired sample t tests were used for continuous variables. The two-gene correlation map was realized by the R package “ggstatsplot”, and the multigene correlation map was displayed by the R package “pheatmap”. We used Spearman’s correlation analysis to describe the correlation between quantitative variables without a normal distribution. The survival analysis was conducted through Kaplan–Meier curves and log-rank tests. All the statistical tests mentioned above are two-sided. P values of < 0.05 were considered statistically significant. Distinctive mark: no significance (ns), p≥0.05; *, p< 0.05; **, p<0.01; ***, p<0.001.

A total of 498 tumor and 52 adjacent normal tissues in the TCGA database were analyzed. Sixty-eight genes were identified after the intersection of CIC-related genes in Genecards (26) and differential and progression-related genes in TCGA. Subsequently, Lasso regression analysis was used to identify 20 candidate genes. Finally, 10 genes were used to establish a predictive model by validating the consistency of DEGs and prognosis. A flow chart of this study was provided in Figure 1. The results of principal component analysis, indicating consistency of gene function and risk group were presented (Figure 2A). We calculated the area under the receiver operating characteristic curve (ROC) of the ten genes, indicating certain accuracy of most genes in predicting tumors and normal tissues (Figure 2B). The DEGs at the mRNA level and at the protein level were presented in Figures 2C–E. In predicting progress, the risk score of the predictive model showed certain accuracy and stability (Figures 2F, G). Moreover, we detected that high-risk patients were more prone to progress than low-risk patients (hazard ratio (HR): 4.11, 95% CI: 2.66-6.37; risk score cut-off: 1.194; Figure 2H). CDKN3, AK5, SLC25A27, TUBB3, and ALB showed better diagnostic values and stability in predicting progression as a result of time-dependent ROC curves (Supplementary Figure 1). The diagnostic and prognostic values and subgroup survival analysis of the ten genes were detailed in the Supplementary Material. The enrolled genes functioned in protein coding. The basic gene information was depicted in Table 1. Likewise, the baseline data of the ten genes in PCA patients based on TCGA database were presented in Table 2. Overall, the differential expression of most genes was related to T stage and Gleason score.

We established a nomogram to predict PFI in PCA patients (concordance (also named C-index): 0.734; 95% CI: 0.707-0.760; Figure 3A). Figures 3B, C showed the calibration (C-index: 0.734; se: 0.026) and decision curve analysis (C-index: 0.734; 95% CI: 0.707-0.760). The risk score plot indicated that high-risk patients were susceptible to progression compared to low-risk patients (Figure 3D). T stage, PSA, and Gleason score were regarded as clinically independent progress factors through the univariable and multivariable Cox regression analyses (Supplementary Table 1). We subsequently integrated the clinical factors into the gene nomogram (C-index: 0.777; 95% CI: 0.752-0.802; Figure 3E), and no significant improvement was observed in predictive ability.

Figure 4A showed good correlations between genes. Additionally, TMB and MSI were positively associated with oncogenes but negatively related to anti-oncogenes (Figure 4A). Correlations among genes, risk score and clinical parameters were summarized in Figure 4B. Patients with higher risk scores were prone to biochemical recurrence, higher Gleason scores, positive residual tumors, positive N stages, higher T stages and older age (Figures 4C–H). CDKN3 and PLP1 presented a better clinical relationship in the subgroup analysis of PFI. The mRNA expression of CDKN3 tended to rise with increasing Gleason score, N stage, T stage, age, residual tumor, and deterioration of therapy outcome. In contrast, PLP1 mRNA expression decreased with increasing age, Gleason score, and T stage. The specific clinical correlations of the included genes were shown in Supplementary Figure 2.

Figure 4I showed the heatmap of the most positively and negatively relevant genes of the enrolled genes. The results of GO and KEGG analysis of each gene were shown in Figure 4J. The enriched pathways of the gene set included choline metabolism in cancer, the Hippo signaling pathway, and the TGF-beta signaling pathway, and GSVA analysis of the top 20 co-expressed genes of the enrolled genes indicated that PCA was positively related to hormone AR and was negatively associated with epithelial-mesenchymal transition (Supplementary Figure 3).

The genetic alteration and mutation location of the ten genes from cBioPortal (TCGA, PanCancer Atlas; RNA seq V2; 493 patients/samples) were shown in Figures 5A–C. Patients in the altered group experienced shorter disease-free survival and PFI than those in the unaltered group (Figures 5D, E). Figure 6 presented the mutation status of genes and its correlation with overall survival (OS) and PFI in PCA patients at the levels of copy number variation (CNV), single nucleotide variation (SNV), and methylation from GSCALite (35). A positive relationship between CNV and mRNA expression was detected in PLP1, ALB and CDKN3 in PCA patients (Figure 6D), and CNV of PLP1 and PTGS2 had shorter OS and PFI than wild-type (WT) PLP1 and PTGS2 (Figure 6E). The SNVs of the investigated genes were shown in Figures 6F–J. GUCY2C, PTGS2, and OPCML showed more effective SNV mutations (Figure 6F), and PTGS2 SNV had a higher risk of OS and progression-free survival (PFS) than PTGS2 WT (Figure 6J). Figure 6K summarized the methylation difference between tumor and normal samples of inputted genes in PCA patients. Higher methylation of DRD5 and PTGS2 and lower methylation of TUBB3, ALB and SLC25A27 were obviously detected. In addition, the mRNA expression of the enrolled genes was negatively associated with methylation (Figure 6L). These results were consistent with the functions of oncogenes and anti-oncogenes. Methylation of PTGS2 and demethylation of SLC25A27 and CDKN3 showed lower PFS than their counterparts in PCA patients (Figure 6M). In terms of the ten-gene signature, the SNV gene set had a higher risk of PFS (Figure 6N) than the WT group. Moreover, the gene set had a higher GSVA score and was prone to progress when compared to a lower GSVA score (Figure 6O).

Figures 7A–C presented the immune score distribution in PCA tissues and normal tissues using three algorithms (Cibersort, Xcell and Mcpcounter). In general, tumor samples presented higher infiltration levels of macrophages, T cells and myeloid dendritic cells than normal samples. In addition, tumor samples had higher immune scores, lower stroma scores and lower microenvironment scores than normal samples. Correlations between genes and immune infiltration in PCA patients from GSCALite (35) were presented at the mRNA level (Figure 7D), CNV level (Figure 7E), and methylation level (Figure 7F). Overall, natural killer T cells, natural killer (NK) cells and cytotoxic T cells were positively related to tumor suppressor genes, while dendritic cells, macrophages, and monocytes were positively associated with oncogenes (Figure 7D). Notably, a positive correlation was observed between natural killer T cells and the mRNA expression of DRD5, OPCML and PLP1. Both dendritic cells and macrophages were positively related to the mRNA expression of AK5 and CDKN3. No obvious relationship was detected between the CNV gene and immune infiltration (Figure 7E). Gene methylation had the opposite effect on immune infiltration compared to gene mRNA expression (Figure 7F), which was consistent with the correlation between gene expression and methylation (Figure 6L). We used ssGSEA to assess the association between the immune infiltration level and genes (Figure 7G), resulting in similar results to the above. In general, increased expression of oncogenes and decreased expression of anti-oncogenes were positively related to a decrease in immune infiltration, especially for NK cells, neutrophils, mast cells, macrophages, stromal score, immune score and ESTIMATE score. In terms of the gene set, mRNA expression of gene set was positively related to macrophage, and was negatively associated with central memory T cell, CD4 naive T cell, T helper type 2, T helper type 17 and CD8+ T cell (Figure 7H). Higher infiltration of natural regulatory T cells and gamma delta T cells was observed in the SNV gene set than in the WT gene set (Figure 7I). In addition, higher infiltration of dendritic cells, natural regulatory T cells, macrophages, monocytes, T helper type 1 cells, and B cells, induced regulatory T cells and the overall infiltration score of 24 immune cells were detected in the CNV gene set than in the WT gene set (Figure 7J). Notably, patients with higher risk scores were associated with significantly lower levels of neutrophils, NK cells, T helper type 1, and mast cells (Figures 7K–N).

Of particular concern was the relationship among genes, TMB and MSI. We found that CDKN3 mRNA expression was positively related to TMB score (Supplementary Figure 3). Negative correlations were observed between the mRNA expression of GUCY2C and PLP1 and the TMB score (Supplementary Figure 3). There was a positive correlation between the risk score and TMB score (Figure 8A), Patients and the diagnostic accuracy of risk score was higher than TMB score for PCA progression (Figure 8B). with higher TMB scores were more prone to progress than those with lower TMB scores (Figure 8C). Likewise, we observed similar results regarding the correlation between the MSI score and risk score and the impact of the MSI score on PFI (Figures 8D–F). Furthermore, the expression distribution of immune checkpoint genes in tumor tissues and normal tissues in PCA patients was determined, and we found that CD274 (also named PD-L1), LAG3, PDCD1LG2 (also named PD-L2), TIGIT and SIGLEC15 were significantly downregulated, while CTLA4 was significantly upregulated (Figure 8G). We subsequently analyzed the correlations between genes and these immune checkpoint genes. We observed that anti-oncogenes presented a significantly positive correlation with PD-L1, PD-L2, TIGIT and SIGLEC15 (Figure 8H). The correlations between gene expression and drug sensitivity in the pancancer analysis of GDSCs and CTRP are shown in Figure 8I and Figure 8J, respectively. 5-Fluorouracil, GSK1070916, KIN001-102, XMD13-2, MK-1775, BRD-K70511574, ISOX, JQ-1, Merck60, PHA-793887, Repligen 136, apicidin, crizotinib, vorinostat, Compound 23 citrate, GW-405833, and PAC1-1 showed relatively good effects.

CLOCK, PER (1, 2, 3), CRY2, NPAS2, RORA, and ARNTL showed a higher correlation with anti-oncogenes, while CSNK1D and CSNK1E presented a greater relationship with oncogenes (Figure 8K). Overall, patients with higher risk scores showed lower mRNA expression of PER1, PER2, and CRY2 and higher expression of CSNK1E (Figures 8L–O). There was a positive correlation between MSI score and CSNK1E (Supplementary Figure 3). In addition, a negative correlation was observed between TMB score and PER (1,2) and CRY2 (Supplementary Figure 3). Specifically, OPCML mRNA expression showed a significantly positive relationship with PER2, CRY2, and RORA. GUCY2C tended to rise with increasing CLOCK, PER2, PER3, CRY2, NPAS2, and RORA (Supplementary Figure 3). An obvious positive correlation was observed between PLP1 mRNA expression and PER1 and PER2 (Supplementary Figure 3). We detected a significantly positive relationship between CLOCK, PER (1,2,3), CRY2, NPAS2, RORA, ARNTL and PTGS2 mRNA expression (Supplementary Figure 3). In contrast, TUBB3 showed an increased tendency with the upregulation of CSNK1D and CSNK1E. SLC25A27 mRNA expression rose with increasing CSNK1E (Supplementary Figure 3).