Identification of a signature of histone modifiers in kidney renal clear cell carcinoma

Data acquisition and acquisition of HMs

MRNA expression profiles, clinicopathological data and prognostic information derived from 539 KIRC samples and 72 normal control samples were obtained from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database [18]. A total of 540 histone modifiers (HMs) were derived from the top published journals [19].

Differential expression analysis and functional enrichment analysis of HMs

Differentially expressed HMs in KIRC samples and normal control samples were screened out based on the screening criteria of |logFC| >1 and false discovery rate (FDR) <0.01 by applying the limma R package. To mine the underlying molecular mechanisms of HMs, GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses were implemented by running R.

Cluster analysis

To explore the relationship between the expression of differentially expressed HMs and the KIRC, consensus clustering analysis was performed in KIRC by applying the limma R package and ConsensusClusterPlus R package. In addition, survival differences between different clusters and clinical characteristics were further analyzed by using the survival, survminer and pheatmap R packages.

Construction and validation of a risk model of HMs-based signature

KIRC samples applied in the construction of the risk model must meet the following criteria: (1) samples have full expression information and (2) the overall survival time of the samples were more than 30 days from the time of diagnosis of KIRC. Then, by applying the caret R package the samples meeting the inclusion criteria were divided into a training set and a test set in a ratio of 7:3 according to the caret sampling method.

Univariate analysis was implemented to identify HMs associated with the prognosis of KIRC in the training set by running R. Then, least absolute shrinkage and selection operator (LASSO), a regularization method in regression analysis, was implemented to eliminate HMs that were overfitted with the model by applying the glmnet R package. The risk score in the model was estimated by the following formula: Risk score = (Coef1 × expression of gene1) + (Coef2 × expression of gene2) + … + (Coef n × expression of gene3). Based on the median risk score of this model, the KIRC samples in the training set can be categorised into high- and low-risk groups. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were used to downscale the expression of all genes in the model. The Kaplan–Meier survival curve and receiver operating characteristic (ROC) curve of HMs-based signature were plotted to display the model’s ability to predict the prognosis of the sample by running R. The test set and the entire set were implemented with the same procedures to examine reliability of this risk model constructed from the samples in the training set.

Construction of the nomogram

The relationship between clinical variables and HMs-based signature was explored. Univariate and multivariate analyses were implemented to find independent prognostic factors for KIRC samples. Then, independent prognostic factors were included into the establishment of nomogram which could forecast 3-year and 5-year overall survival (OS) in KIRC, and the calibration curves visualising the comparison between the 3-year and 5-year survival probabilities forecasted by the nomogram for KIRC and the observed actual survival probabilities were then plotted.

Immune characteristics of the HMs-based signature

KIRC is an immunogenic tumor and its development and progression are associated with immune characteristics in the tumor microenvironment. Based on this basis, the relationship between the expression level of 8 HMs in the risk model and the level of immune cell infiltration was analysed in the TIMER database (https://cistrome.shinyapps.io/timer/) [20]. The relationship of HMs with microsatellite instability (MSI) and tumor mutation burden (TMB) was further analyzed given that MSI and TMB can respond to immunotherapy of tumors. Then, the differences in tumor microenvironment (TME) between high- and low-risk groups in the risk model were explored. Since the level of immune checkpoints can reflect the immune status in tumor, the characteristics of expression of immune checkpoints between high- and low-risk groups were further explored. A major cause of immunotherapy failure and tumor progression is immune escape. The tumor immune dysfunction and exclusion (TIDE, http://tide.dfci.harvard.edu/) [21] score is a good indicator to evaluate tumor immune escape [21]. Therefore, differences in TIDE scores between high- and low-risk groups were analyzed to assess the probability of immune escape in different KIRC samples.

Correlation between drug sensitivity and risk signature

To investigate the sensitivity of KIRC samples to common chemotherapy drugs between different risk groups, the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) database [22] and the pRRophetic R package were applied to analyze differences in IC₅₀ of KIRC samples in chemotherapy drugs between different risk groups.

Statistics analysis

All statistical analyses were implemented by R (version 4.0.0). Wilcoxon Rank-Sum test was implemented to compare gene expression differences between KIRC samples and control samples. Pearson’s correlation coefficient was applied to estimate the linear relationship between the variables. P-values less than 0.05 were considered significant differences.

Availability of data and materials

Raw data in this study can be obtained from the GDC database (https://portal.gdc.cancer.gov/), TIDE database (http://tide.dfci.harvard.edu/) and GDSC database (http://www.cancerrxgene.org/).

Identification and functional enrichment of differentially expressed HMs

Of the 539 KIRC samples and 72 control normal samples, 21 down-regulated and 27 up-regulated HMs were identified and the results were presented in Figure 1A, 1B. GO analysis (Figure 2A) showed that the GO terms with the highest enrichment in the biological process (BP) classification were histone modification, chromatin organization, peptidyl−lysine modification, chromatin remodeling and histone methylation. In the cellular component (CC) content, condensed chromosome, chromosomal region and spindle were the main enriched GO terms. Histone binding, methyltransferase activity, transferase activity and transferring one−carbon groups were the significantly enriched GO terms in the category of molecular function (MF). The enrichment analysis of KEGG pathways implied that HMs may be closely associated with Lysine degradation (Figure 2B).

Classification of KIRC samples-based HMs signature

When performing consensus cluster analyses on the KIRC samples by clustering variable (k) from 1 to 9, we found that the strongest correlations between samples within groups were found when k = 2 (Figure 3A). Subsequently, the 515 KIRC samples were divided into two clusters. Kaplan-meier (KM) survival analysis between the two clusters showed that OS of KIRC samples in cluster 1 was significantly worse than OS in cluster 2 (p < 0.001) (Figure 3B). We present a heat map depicting the association between the mRNA expression of HMs and the clinical characteristics of KIRC samples, including Age, Gender, Grade and Stage (Figure 3C). We found significant differences between the KIRC samples in Cluster 1 and Cluster 2 in terms of Gender (p < 0.05), stage (p < 0.001) and grade (p < 0.05) of the tumor.

Construction and verification of HMs-based signature

In all, 515 samples from the KIRC met the inclusion criteria for follow-up research. These KIRC samples were further divided into 361 samples as a training set and 154 samples as a test set in a 7:3 ratio by care method sampling. Univariate analysis of differentially expressed HMs in the training set were performed and 21 HMs that significantly correlated with the prognosis of KIRC were identified (Figure 4A). Those HMs that were over-fitted to the model were eliminated, and finally 8 HMs (C17orf49, GLYATL1, HJURP, KAT2A, NCOA7, NEK6, PRDM16, TTK) were identified to build a risk model (Figure 4B, 4C). The risk score can be calculated from the mRNA expression of HMs and the relevant coefficients (Table 1) with the following equation: riskscore = 0.0280971627824297 × C17orf49 expression + (−0.107290283832969) × GLYATL1 expression + 0.520165069344535 × HJURP expression + 0.241951166184651 × KAT2A expression + (−0.00243111668616741) × NCOA7 expression + (−0.13291750024586) × NEK6 expression + (−0.500207346820667) × PRDM16 expression + 0.209483791892391 × TTK expression. As a cut-off value, the median risk score was used to divide the KIRC samples into high-risk and low-risk categories. In the two-dimensional plane, PCA and t-SNE analysis (Figure 5) showed significant dispersion between KIRC samples of high- and low-risk groups. The mRNA expression of the 8 HMs in the high-risk group and the low-risk group was shown in Figure 6A. The correlation between the risk signature and the survival of the KIRC samples were shown in Figure 6B–6D. The KIRC samples in the high-risk group had significantly shorter survival time and significantly higher mortality rates than the samples in the low-risk group. The ROC curve showed an accuracy of 0.706 and 0.748 for the risk score to predict 3- and 5-year OS for the KIRC samples, respectively (Figure 6E). A same analysis was performed on the test set as well as the entire set to verify the stability and reliability of the risk model derived from training set. As shown in Figures 7A, 8A, test and overall sets of KIRC samples of high- and low-risk groups displayed similar expression level of HMs as the training samples. Survival differences between high-risk group and low-risk group in the test (Figure 7B–7D) and overall set (Figure 8B–8D) were consistent with those in the risk model. The accuracy of the test set’s risk score in predicting 3- and 5-year OS for KIRC samples was 0.738 and 0.749, respectively (Figure 7E) and 0.702 and 0.722, respectively in the entire set (Figure 8E). Analyses of the test set and the total set showed good stability and reliability of the risk model.

Relationship between HMs-based signature and clinical variables

The KIRC samples in the risk model can be divided into subgroups based on age, gender, grade, and stage. The heatmap displayed a significant difference in grade (p < 0.01) and stage (p < 0.001) between high- and low-risk groups (Figure 9A). The KM survival curves in each subgroup (age ≤65 or >65, gender = FEMALE or MALE, grade = G1−2 or G3−4, stage = I−II or III−IV) showed that each subgroup of KIRC samples with a high risk score had a worse OS than those with a low risk score (Figure 9B–9E).

Construction of prognostic models

A univariate Cox analysis was conducted using variables including age, gender, grade, stage, and risk score (Figure 10A). Univariate Cox analysis revealed close correlations between age (p < 0.001), grade (p < 0.001), stage and risk score (p < 0.001) for the prognosis of KIRC samples. Further multivariate Cox analyses were implemented using age, grade, stage, and risk score. Cox analysis of the KIRC samples revealed that age (p < 0.001), grade (p = 0.053), stage (p < 0.001), and risk score (p < 0.001) influenced prognosis independently (Figure 10B). In addition, ROC curves of several clinical variables and risk score predicting the OS of KIRC samples shown that risk score (AUC = 0.768) can be as accurate as or better than some traditional clinical variables in predicting KIRC samples’ prognosis (Figure 10C).

As a result of the univariate and multivariate analyses, independent prognostic factors of the KIRC samples were incorporated into the prognostic model. By transforming complex regression equations into visual graphs, nomograms make the results of predictive models more understandable [23]. A nomogram predicting the probability of 3- and 5-year OS was plotted on the basis of this advantage (Figure 11A). The concordance index (C-index) = 0.754 reflected the decent accuracy of the nomogram in predicting the 3- and 5-year OS of KIRC. The calibration curves further showed the consistency between the nomogram-predicted OS probabilities for the KIRC samples and the actual OS for the KIRC samples at 3- and 5-year (Figure 11B, 11C).

Immune characteristics

By searching the TIMER database, we found correlations between expression level of these 8 HMs (C17orf49, GLYATL1, HJURP, KAT2A, NCOA7, NEK6, PRDM16, TTK) and the infiltration level of six common immune cells (B Cell, CD8+ T Cell, CD4+ T Cell, Macrophage, Neutrophil and Dendritic Cell) (Figure 12). The expression of C17orf49 closely correlated with the infiltration of CD8+ T Cell (cor = 0.35, p = 4.32e−14), CD4+ T Cell (cor = 0.295, p = 1.15e−10), Neutrophil (cor = 0.16, p = 5.76e−04), and Dendritic Cell (cor = 0.289, p = 3.34e−10). The expression of GLYATL1 significantly correlated with the infiltration of B Cell (cor = 0.161, p = 5.29e−04) and CD8+ T Cell (cor = 0.152, p = 1.44e−03) only. HJURP showed a close association with immune infiltration, and its expression level significantly correlated with B Cells (cor = 0.261, p = 1.32e−08), CD8+ T Cell (cor = 0.187, p = 8.21e−05), CD4+ T Cell (cor = 0.19, p = 4.27e−05), Macrophage (cor = 0.146, p = 1.88e−03), Neutrophil (cor = 0.281, p = 9.02e−10) and Dendritic Cell (cor = 0.347, p = 2.40e−14). The expression level of KAT2A had significantly negative correlation with the infiltration level of B cells (cor = −0.153, p = 1.01e−03), however, it had significantly positive correlation with the infiltration level of CD4+ T cells (cor = 0.323, p = 1.18e−12) and Neutrophil (cor = 0.129, p = 5.53e−03). The expression of NCOA7 closely correlated with infiltration of CD8+ T Cell (cor = 0.202, p = 1.99e−05), CD4+ T Cell (cor = 0.247, p = 7.79e−08), Macrophage (cor = 0.214, p = 5.02e−06), Neutrophil (cor = 0.245, p = 1.11e−07) and Dendritic Cell (cor = 0.117, p = 1.22e−02). Both NEK6 and TTK show a close correlation with immune cell infiltration. The expression of NEK6 and TTK closely and positively correlated with the infiltration of these six immune cells (B Cells, CD8+ T Cell, CD4+ T Cell, Macrophage, Neutrophil and Dendritic Cell) (p < 0.01). The expression of PRDM16 had significantly negative correlation with the infiltration of B cells (cor = 0.13, p = 5.33e−03), but significantly positive correlation with the infiltration of CD8+ T Cell (cor = 0.1, p = 3.72e−02), CD4+ T Cell (cor = 0.311, p = 8.90e−12) and Macrophage (cor = 0.135, p = 4.21e−03). We further investigated the correlation between the expression of these 8 HMs and the level of MSI and TMB. The results showed that the expression of C17orf49 significantly correlated with the level of TMB (r = 0.11, p = 0.041). The expression of HJURP significantly correlated with the level of MSI (r = 0.15, p = 0.005) and TMB (r = 0.24, p = 1.49e−05). KAT2A likewise showed a close association with MSI (r = 0.17, p = 0.002) and TMB (r = 0.16, p = 0.003). The expression level of PRDM16 negatively correlated with the level of TMB (r = −0.12, p = 0.032). The expression level of TTK showed a meaningful correlation with the level of both MSI (r = 0.17, p = 0.002) and TMB (r = 0.20, p = 2.39e−04) (Figure 13).

Based on the evidence that these HMs were associated with immunity, we analyzed the correlation between HMs-based signature and immune characteristics. Analysis of differences in the TME between high- and low-risk groups showed that the level of immune cell infiltration in the TME was significantly higher in the high-risk group than in the low-risk group (p = 1.1e−05) (Figure 14A–14C). We further analyzed the differences in the expression of immune checkpoints in KIRC samples between high- and low-risk groups, and result showed that the level of expression of IDO2, ICOS, PDCD1, CD70, LAIR1, CD28, CD40, TNFRSF4, CD160, ADORA2A, TNFSF9, LAG3, BTLA, CD48, CD44, CD40LG, TIGIT, TNFSF4, TMIGD2, TNFRSF14, TNFSF14, LGALS9, TNFRSF9, CD86, CD244 and TNFRSF25 were significantly higher in the high-risk group than in the low-risk group (Figure 14D). In addition, a significant finding was that the TIDE scores of KIRC samples in the high-risk group were significantly higher than those of samples in the low-risk group (Figure 14E), revealing that KIRC in the high-risk group was more prone to immune escape.

Drug sensitivity analysis

To improve chemotherapy outcomes in KIRC patients, it is important to understand KIRC’s sensitivity to chemotherapy. Analyses of chemotherapy drug sensitivity in KIRC samples with high or low risk showed that KIRC samples with high risk had lower IC₅₀ values for Sunitinib, Tipifarnib, Nilotinib, Bosutinib, Mitomycin. C, Vinblastine, Camptothecin, Docetaxel and Doxorubicin compared with KIRC samples with low risk (Figure 15), implying that KIRC samples with high risk were more sensitive to these nine chemotherapy drugs. In contrast to the high-risk group, the low-risk group had significantly lower IC₅₀ values for Pazopanib, Bexarotene, and Thapsigargin in KIRC, implying that KIRC samples with low risk were more sensitive to these three chemotherapeutic drugs.