Identification of KIRC Tregs-related genes (KTRGs)
Through the ssGSEA analysis, we figured out the Tregs activity score of KIRC. Based on the optimal cut-off value, the KIRC samples were separated into high- and low-Treg activity groups. We observed a poor prognosis that occurs in high Tregs activity groups (Fig. 1A), which lures us to further explore the potential reason and mechanism of Tregs on KIRC prognosis. The first thing we did was to determine the Tregs-related genes. To do this, we calculated the Pearson correlation coefficients (PCC) between each gene expression and Tregs score based on the whole genome expression profile of KIRC(TCGA). Those with absolute PCC > 0.4 and P-value < 0.05 were identified as Treg-related genes (Fig. 1B). Then we carried out differential expression analysis by comparing the tumour with the normal expression profile of KIRC with a filtering threshold based on absolute log2(Fold Change) > 2 and adjusted P-value < 0.05, differential expression genes (DEGs) were finally screened (Fig. 1B). 76 KIRC Tregs-related genes (KTRGs) were eventually determined by overlapping Treg-related genes and DEGs. The heatmap revealed that these genes' expressions could rightly discriminate between tumour and normal samples (Fig. 1C). When the KTRGs expression profile was conducted dimensionality reduction by t-distributed Stochastic Neighbor Embedding (t-SNE), we also observed tumour and normal samples harboured in different areas in a two-dimension axis, indicating that the expression of these genes could distinguish the tumour and the normal tissues well (Fig. 1D). We then presented the mRNA expression profile of KTRGs to clinical survival data, and a univariate Cox regression analysis was subsequently performed. 22 genes of 76 KTRGs were estimated to be associated with the overall survival of KIRC. Among these, SPI1, TYMP, LILRB1, IL20RB, ITGAX, C1QL1, FCER1G, SLAMF8, TGFBI, CCL5, and CXCL13 were risky genes to the prognosis of KIRC patients, high expression of which represent a worse survival probability (Fig. 1E and Figure S1). While OGDHL, MTURN, HSD11B2, PPARGC1A, LDHD, FREM1, NE3C2, TRIM2, SLAMF8, AGPAT9, and WDR72 were protective genes to OS of the patients and low expression of which represent a worse clinical outcome. In the protein-protein interaction (PPI) network, HSD11B2, NR3C2, and PPARGC1A show a protein interaction, SLAMF8, FCER1G, ITGAX, LILRB1, SPI1, CCL5, and CXCL13 show a protein interaction, while others do not show a relationship (Figure S2A). In mRNA levels, WDR72, MTURN, NR3C2, AGPAT9, TRIM2, PPARGC1A, OGFHL, LDHD, FREM1, and HSD11B2 show positive correlation with each other and negative correlation with other genes (Figure S2B). We further analyzed the epigenetic characterization, and we found that the mutation rate of all KTRGs is low (i.e., mutation rate < 5%), but the copy number variation (CNV) of some is exuberant (Figure S2C and D). Among these, TGFBI, SLAMF8, MTURN, LILRB1, LDHD, ITGAX, IL20RB, and HSD11B2 show high amplification (Amp) rate, TRIM2, OGDHL, NR3C2, IL20RB, FCER1G, and ALDH6A1 show high deletion (Del) rate. The expression of IL20RB, TRIM2, NR3C2, CXCL13, SLAMF8, and FCER1G are affected by copy number variation, while other genes do not show a significant effect.
Clinical relevance and mechanisms of KTRGs
Apart from the relationship between prognosis and KTRGs mentioned, we matched the clinical stage and grade data and mRNA expression profile of KTRGs. Analysis of variance (ANOVA) was performed to estimate the difference of genes in different groups (i.e., Grade I - IV). Gene expression levels were shown by heatmap. We found OGDHL, PPARGC1A, LDHD, AGPAT9, ALDH6A1, WDR72, HSD11B2, FREM1, MTURN, NR3C2, and TRIM2 decrease little by little with the advance of grade and stage, whereas TGFBI, IL20RB, C1QL1, CCL5, CXCL13, TYMP, ITGAX, LILRB1, SLAMF8, SPI1, and FCER1G, gradually increase with development of disease (Fig. 2A and B). Through analyzing the correlation between 50 hallmark pathway activities and KTRGs expression, we found several genes (i.e., TGFBI, SPI1, SLAMF8, LILRB1, etc.) expression positively correlated with cancer-related pathways such as G2M checkpoint, E2F targets, Hedgehog signalling, Wnt pathways, P53 pathways, etc. PPARGC1A, OGDHL, ALDH6A1, and AGPAT9 have a positive correlation with adipogenesis, fatty acid metabolism, heme metabolism, bile acid metabolism, androgen response (Fig. 2C). PROGENy algorithm (Pathway RespOnsive GENes) concludes several cell-specific signalling pathways [15]. Consistently, TGFBI, SPI1, SLAMF8, LILRB1, FCER1G, CXCL13, CCL5, and C1QL1 show a positive correlation with most cancer-related pathways. WDR72, TRIM2, PPARGC1A, FREM1, ALDH6A1, and AGPAT9 show a negative correlation with most cancer pathways but a positive correlation with androgen pathways (Fig. 2D). All in all, these KTRG expressions are related to tumour progression and might function differently by regulating specific metabolisms.
KTRGs risk signature construction
To further identify the potential important KTRGs, we matched the KTRGs expression profile with KIRC clinical data, and then the data was submitted to the Boruta algorithm. Three types of KTRGs, including confirmed (i.e., important) genes, tentative genes (i.e., probably important), and rejected (i.e., unimportant) genes, were finally determined. Among these, CXCL13, WDR72, TRIM2, ALDH6A1, NR3C2, PPARGC1A, and IL20RB were defined as confirmed genes, SLAMF8, FCER1G, LDHD, MTURN, TYMP, and OGDHL were defined as tentative genes, while others were considered as rejected genes (Fig. 3A). Then, we separated the KIRC(TCGA) datasets into two cohorts, training and testing cohorts, according to the ratio of 7:3. We dismissed the rejected genes and abstracted the expression profile of confirmed and tentative genes from the training cohort to conduct lasso-cox regression analysis. Five genes were finally screened, including NR3C2, WDR72, ALDH6A1, CXCL13, and IL20RB. Interestingly, the genes identified by the lasso-cox algorithm are all confirmed genes, indicating the importance of these genes for prognosis under these two algorithms. Based on the lasso coefficient per gene, a risk score was finally calculated according to the equation:
KTRGs risk score = 0.01755748*CXCL13–0.15617289*WDR72–0.04032272*ALDH6A1–0.16222042*NR3C2 + 0.01788896* IL20RB
According to the ideal cut-off value estimated by the surv_cutpoint function from the survminer R-package, the KIRC training cohort was divided into high and low-expression groups. Survival curve analysis shows that high risk group patients are faced with the problem of low survival probability (Fig. 3B). From ROC curve analysis, we perceived the area under curves (AUC) at 3-/5-/7-year were 0.6915, 0.7025, and 0.7050, indicating that the KTRG risk score is highly accurate in predicting the prognosis of KIRC patients. The heat map shows that NR3C2, WDR72, and ALDH6A are highly enriched in the low-risk group, while CXCL13 and IL20RB are highly enriched in the high-risk group. The same analysis was performed in KIRC testing cohort, and we observed a similar result (i.e., For ROC, AUC at 3 years = 0.7212, AUC at five years = 0.7392, AUC at 7 years = 0.7107). We utilized other two RCC datasets to further validation, and the results bear resemblance to KIRC datasets (i.e., For E-MTAB-1980 ROC, AUC at three years = 0.7309, AUC at five years = 0.7455, AUC at 7 years = 0.6919; For KIRP ROC, AUC at three years = 0.5940, AUC at five years = 0.6196, AUC at seven years = 0.6097) (Figure S3A and B). These results revealed that the KTRG risk score is practical and accurate for predicting outcomes of patients with kidney cancer.
Clinical relevance and robustness of KTRGs risk score
To further explore the relationship between KTRG risk score and clinical variables, we matched the clinical variables (i.e., pathological stage, histological grade, Gender, and Age) with a 5-gene KTRG expression profile in two RCC cohorts. We found that KTRG risk groups were not correlated with patients' age in the KIRC cohort but correlated with stage and grade, appearing that more advanced stage and grade harboured in high-risk groups (Fig. 4A). CXCL13 and IL20RB are highly enriched in high-risk groups, indicating that the expression of these two genes might promote disease progression. We then divided the clinical data into different clinical sub-groups and matched them with the corresponding mRNA expression profile of KTRGs. Intriguingly, we found that, whatever the subgroup is, patients' outcomes in the high-risk group are always poorer, indicating that the KTRG risk score is robust and steady, which could predict the prognosis of patients in every sub-classes of kidney cancer (Figure S4). The same analysis was performed in the E-MTAB-1980 cohort, and we observed the same results (Fig. 4B and Figure S5).
Clinical nomogram model construction
Next, we utilized the clinical variables (i.e., age, gender, stage, and grade) to adjust the KTRG risk score. In multivariate Cox regression analysis, we found that the KTRG risk score still acts as a risky factor for the prognosis of patients, which is not affected by other clinical variables, suggesting that our established risk score is an independent predictor (Figure S7). Then, we conducted the stepwise regression analysis based on the minimum Akaike information criterion (AIC) value in the KIRC(TCGA) dataset; four variables, including risk score, age, grade, and stage, were selected to construct a clinical nomogram model to predict the 3-/5-/7-year survival probability (Fig. 5A). The calibration curve and ROC curve revealed that the model possesses a high discriminative accuracy (AUC at three years = 0.8081, AUC at five years = 0.761, AUC at seven years = 0.7394) (Fig. 5B and C). The nomogram risk score was further calculated based on the established model, and the high and low model risk scores showed extremely significant differences in overall survival probability (i.e., after 100 months, the survival probability of patients with a high model risk score is proximately zero, while in low model risk score are fifty per cent or so) (Fig. 5D). The result reveals that the nomogram risk model could stratify the patients with different prognoses more effectively and discriminatively. This is indeed the case when observing decision curve analysis (DCA). The nomogram model shows a maximum benefit in predicting a 3-/5-/7-year prognosis of patients compared to other single variables (Fig. 5E). The same analysis recurred in E-MTAB-1980 datasets, and the results are consistent with KIRC(TCGA) dataset (i.e., for ROC curve, AUC at three years = 0.8974, AUC at five years = 0.8805, AUC at seven years = 0.8241; for survival curve, the survival rate of patients in high-risk group is forty per cent while in low-risk group are eighty per cent about, after fifty mouths) (Figure S6).
Differentially enriched pathways between high and low KTRGs risk groups
Given the differential outcomes of the high and low KTRG risk groups, we try to explore the potential distinct mechanisms between these two groups. We first submitted the mRNA expression profile of KIRC(TCGA) to gene set enrichment analysis (GSEA) software, and then the potential cancer-related hallmark pathways were further explored. We found several pathways highly enriched in the high KTRGs risk group, and the top three are the EMT pathway, Allograft rejection, and IL6 JAK STAT3 signalling (Fig. 6A). Several metabolism-related pathways are highly enriched in the low KTRGs risk group, such as heme metabolism, fatty acid metabolism, and adipogenesis (Fig. 6B). Further analysis combining the pathways activity score and risk score to conduct Pearson correlation coefficient, the results represent that KTRGs risk score highly positively correlated with MYC TARGET V2, IL6 JAK STAT3, ALLOGRAFT REJECTION, DNA REPAIR, and UNFOLDED PROTEIN RESPONSE, and negatively correlated with FATTA_ACID_METABOLISM, PROTEIN_SECRETION, ADIPOGENESIS, ANDROGEN_RESPONSE, HEME METOBOLISM, etc. (Figure S8). When analyzing the biological process, we found several immune processes highly enriched in high KTRG risk groups, including antimicrobial humoral response, antibacterial humoral response, and regulation of acute inflammatory response to antigenic stimulus (Fig. 6C). Furthermore, the metabolism process, in terms of lipid modification, lipid oxidation, and negative regulation of the amide metabolic process, is mainly enriched in low KTRG risk groups (Fig. 6D). Considering the metabolism-related pathways and processes consistently differential enrichment, we collected some metabolism-related signatures such as energy, amino acid, vitamin, nucleotide, carbohydrates, and TCA and submitted them onto GSEA software. The results are consistent, and all these metabolism signatures are highly enriched in low KTRG risk groups (Fig. 6E and F). These results suggested several hallmark pathways, metabolism signatures, and some immune mechanisms occurred differentially enrichment between the two risk groups.
Immune relevance and immune therapeutic response of risk score
Previous research pointed to the effect of Tregs on the immune system and even on immune therapeutic response. Respecting the differential enrichment of the immune process mentioned above we had analyzed, we further explored the relevance of immune and KTRG risk scores. We found the fraction of some anti-tumor immune cells, such as T cells CD8 and T cells CD4, is higher, while B cell naive, NK cell activated, and M1 macrophages are lower in the high KTRGs risk group. Needless to say, the Tregs fraction is higher in the high-risk group. However, the pro-tumor immune cells - M2 macrophage fraction decreased in the high-risk group (Fig. 7A). For immune cell activity, interestingly, B cells, CD8 + T cells, NK cells, Cytotoxic cells, and Tregs consistently increase in high-risk groups. There is no doubt that the KTRG risk score shows the highest positive correlation with the Tregs activity score but shows a less positive correlation with other cells, such as CD8 T cells and NK cells (Fig. 7B). In addition, the KTRG risk score shows a negative correlation with several immune cell fractions, including monocytes, NK cells resting, and M1/M2 macrophages (Figure S9A). Besides Tregs, the risk score also shows a low positive correlation with the fraction of CD8 T cells and T cells CD4 memory activated. In addition, the immune checkpoint, CTLA4 and PDCD1 (PD1) are highly expressed in the high-risk group, but CD274 (PDL1) is highly expressed in the low-risk group (Figure S9B). The risk score is highly positively correlated with mRNA expression of PDCD1 and CTLA4 but negatively correlated with CD274 (Figure S9C). All in all, the relationship between the KTRG risk score and other immune cells is extremely complex, and we can not predict the potential immune therapeutic response among these results. Therefore, we conducted immune therapeutic response prediction by TIDE analysis. Notably, we found that most non-responders were distributed in the high-risk group, whereas most of the respondents were located in the low-risk group (Figure S9D). In addition, the risk score is higher in no-responders than that in responders (Figure S9E). From this humble evidence, we could speculate that the KTRG risk score could act as a potential predictor for cancer patients who receive immune therapy. We obtained two external immune therapy datasets (Alexandra and GSE78820). After merging the risk score based on their own KTRGs mRNA expression profile and clinical variables (such as survival data and immune therapeutic response data), the subsequent analyses were further performed. Similar to the results above, the high-risk group showed a low survival probability in these two datasets (Fig. 8A and E). Risk score also shows high predicted accuracy (for the Alexandra cohort, AUC at one year = 0.6410, AUC at three years = 0.7868; for the GSE78820 cohort, AUC at one year = 0.7072, AUC at three years = 0.7229) (Fig. 8B and F). Furthermore, the patients that prove to be non-responders (SD/PD) show a higher KTRGs risk score (Fig. 8C and G). In the same way, the patients in the high-risk group primarily harboured in the no-responder group, while those in the low-risk group were mainly distributed in the responder (CR/PR) group (Fig. 8D and H). These results provide potential evidence for our hypothesis that our KTRG risk score could be employed to act as a predictor for predicting immune therapeutic outcomes.