To compare immune-autophagy-related differentially expressed genes (IAR-DEGs) in KIRC, the Gene GEO database was used. The GSE186645 dataset was selected for subsequent analyses. A total of 1,793 human immune-related genes (IRGs) were downloaded from ImmPort database (https://www.immport.org./home), and a total of 223 human autophagy-related genes were downloaded from the Human Autophagy Database (HADb) (http://autophagy.lu/clustering/index.html). The cutoff conditions were set to an adjusted p-value <0.05, and the absolute value of log-fold change | log2FC| ≥ 1 was statistically significant for the DEGs. ImageGP was used to create volcano maps and venn maps online.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed by ClusterProfiler software package to explore functional annotation and enrichment pathways, with p < 0.05 representing statistically significant differences.

In order to further evaluate the expression and prognostic value of IAR-DEGs in KIRC, differential analysis and prognostic analysis through “survival” package were conducted. Based on the Cox proportional hazards model and Kaplan–Meier model, the hazard ratio (HR) was calculated, with p < 0.05 representing statistically significant differences.

According to the preliminary screening of IAR-DEGs with differentially expressed and prognostic significance, univariate Cox analysis of overall survival (OS) was performed to identify the survival-related IAR-DEGs with a significant prognosis value (p < 0.05). Then, multivariate Cox regression analysis was performed to construct a prediction model based on IAR-DEGs, and the IAR-DEGs were independent prognostic factors. Signatures were established based on the coefficients corresponding to independent prognostic genes. Patients from TCGA-KIRC dataset were divided into low- and high-risk groups weighted by the risk score obtained from the multivariate Cox regression. t-Distributed stochastic neighbor embedding (t-SNE) and principal component analysis (PCA) were used to explore the distribution characteristics of different groups by R packages. Finally, the effectiveness of prognostic indicators was evaluated by the area under the curve (AUC) of “time receiver operating characteristic curve (ROC).”

Based on “survival” package in R software, combined with the clinicopathological characteristics, the correlation between IAR-DEGs and clinicopathological characteristics was analyzed. Through R package “rms,” the nomogram and calibration curve were obtained. Risk scores associated with prognostic models were used as prognostic factors to evaluate 1-, 3-, and 5-year OS.

The relationship between IAR-DEGs expression levels and immune cells was analyzed using the xCell algorithm in the “immunedeconv” R package. The immune score and the effects of gene expression levels on eight immune checkpoint-related genes were also analyzed using the “ggplot2” R package. Finally, TIDE algorithm was used to evaluate two different mechanisms of tumor immune escape using IAR-DEG markers.

The third-order RNA sequencing data of genes were obtained based on TCGA dataset, and the association with ferroptosis-related genes and m⁶A-related genes in “ggplot2” R package was analyzed.

Tumor-associated RNA-seq data were obtained from TCGA-KIRC, mRNAsi was calculated by one-class logistic regression (OCLR) algorithm, and the dryness index was obtained.

Human kidney cell line HK-2 and human KIRC cell lines, 786-O and caki-1, were originally purchased from the cell repository of Shanghai Institute of Life Sciences. The cells were cultured in 1640 Medium (Gibco, Grand Island, NY, USA), containing 10% fetal bovine serum (FBS) (Gibco), penicillin (25 U/ml), and streptomycin (25 mg/ml), with 5% CO₂ environment.

In this study, 19 fresh samples, including tumor tissues and adjacent normal kidney tissues, were collected from patients who underwent laparoscopic radical nephrectomy for KIRC from 2019 to 2020 in the Department of Urology, Zhongda Hospital, and stored at 80°C. All patients were diagnosed with KIRC and did not receive any antitumor therapy preoperatively, and none of them had a history of long-term drug use. The clinical characteristics of 19 clear cell RCC (ccRCC) patients are listed in Table 1. The methodology of this study followed the criteria outlined in the Declaration of Helsinki (revised in 2013), and ethical approval was obtained from the Ethics Committee and Institutional Review Board for Clinical Research of Zhongda Hospital (ZDKYSB077). All patients or their relatives who participated were informed and signed an informed consent form.

Total RNA was isolated with Total RNA Kit (OMEGAbiotec, Guangzhou, China) according to the manufacturer’s instructions. Complementary DNA was synthesized using the HiScript II Q RT SuperMix (R223-01) reagent kit (Vazyme Biotech Co., ltd., Nanjing, China). The qRT-PCR was performed using the SYBR green PCR mix (vazyme). The specific primers set for mIR-DEGs and GAPDH are listed in Supplementary Table S1. Data were normalized to GAPDH expression levels using the 2^−ΔΔCt method.

All specimens were fixed in 10% neutral formaldehyde solution and embedded in paraffin. Envision two-step dyeing and DAB color development were used. Primary antibodies BID (ab32060, Abcam, Cambridge, UK), NAMPT (ab236874, Abcam), and BIRC5 (ab76424, Abcam) were used in this study.

Total proteins from HK-2 and human KIRC cells lysed in radioimmunoprecipitation assay (RIPA) (KeyGen, Nanjing, China) buffer were extracted and quantified by bicinchoninic acid (BCA) assay (KeyGen, China). Proteins were analyzed by 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE), and the gels were transferred onto polyvinylidene fluoride (PVDF) membranes. Then, bovine serum albumin (BSA)-blocked PVDF membranes were incubated with specific primary antibodies BID (1:1,000; ab32060), NAMPT (1:1000; ab236874), BIRC5 (1:5,000; ab76424), and CANX (1:2000; ab133615) overnight at 4°C, followed by incubation of secondary antibodies for 1 h. Finally, bands were visualized using an enhanced chemiluminescence (ECL) kit (vazyme, China).

The statistical analysis was carried out by R software (version 4.0.2). The Perl programming language (version 5.30.2) was used for data processing. Multivariate Cox regression analyses were used to evaluate prognostic significance. When p < 0.05 or log-rank p < 0.05, the difference was statistically significant.

The volcano map shows 1,826 upregulated DEGs and 1,809 downregulated DEGs that we screened in GSE168845 (Figure 1A). Then, 1,793 human IRGs from ImmPort database and 223 human autophagy-related genes from HADb were analyzed by Venn diagram, and five co-expressed genes were obtained: CANX, MAPK1, BIRC5, NAMPT, and BID (Figure 1B). In the GO/KEGG pathway enrichment analyses, we found five co-expressed differential genes enriched in “aging” in biological process (BP); “dendrite cytoplasm,” “neuron projection cytoplasm,” and “plasma membrane projection cytoplasm” in cellular component (CC); and molecular function (MF) enriched in “MAP kinase activity,” “MAP kinase activity,” and “death receptor binding.” Importantly, the five co-expressed genes in KEGG were mainly enriched in “Platinum drug resistance,” “Apoptosis,” and “Apoptosis-multiple species” (Figure 1C).

Through a screening in TCGA-KIRC database, we compared the expression levels of CANX, MAPK1, BIRC5, NAMPT, and BID in normal kidney tissues and renal clear cell tumor tissues, and we found that their expression levels in tumor tissues were upregulated (Figure 2A). And Kaplan–Meier model analysis shows that the expression levels of the above five DEGs are significantly related to the prognosis, the high expression of CANX and MAPK1 is associated with a good prognosis (Figures 2C,F), and the high expression of BID, BIRC5, and MAPK1 is associated with a poor prognosis (Figures 2B,D,E). Univariate Cox regression analysis (Figure 3A) and multivariate Cox regression analysis (Figure 3B) were used to further explore the correlation between the five DEGs and prognosis, showing that CANX, BIRC5, NAMPT, and BID are independent prognostic factors for KIRC.

We used lasso Cox regression to construct a prognostic model of DEG-related risks, Risk Score = (−0.4879) * CANX + (0.3075) * NAMPT + (−0.3041) * BIRC5 + (0.694) * BID (Figure 4A, Figure 4B). According to the median risk score (50%), patients were divided into high-risk and low-risk groups. It can be seen in the t-SNE and PCA heat maps that BID, BIRC5, and NAMPT are highly expressed in the high-risk group, and CANX is low in the high-risk group (Figure 4C). If HR = 2.333, the prognosis model can be considered as a risk factor model. The median survival time of the high-risk group was significantly lower than that of the low-risk group (Figure 4C). We used ROC to evaluate the prognostic prediction efficiency of the model, and the results showed that the AUC was 0.73 (1-year OS), 0.685 (3-year OS), and 0.697 (5-year OS) (Figure 4C).

Regarding the correlation between CANX, BID, NAMPT, BIRC5, and clinicopathological characteristics in the risk prognosis model, our results show that the immune-autophagy-related DEGs associated with T stage, N stage, M stage, and pathological stage are BIRC5 and BID (Figures 5A–D), there is no age-related gene, and the gene related to the patient’s gender is BIRC5 (Figures 5E,F). We used the nomogram to predict 1-, 3-, and 5-year OS in the entire TCGA cohort (Figure 5G). We also found that the 1-, 3-, and 5-year OS on the nomogram is consistent with the calibration curve of predicted probability, and the 1-year OS is the highest (Figure 5H).

We used Gene Set Enrichment Analysis (GSEA) to analyze the KIRC patient data in TCGA-KIRC database. The results showed that both BIRC5 and BID mediate ion channel transport (Figures 6A,B), and both NAMPT and CANX mediate the channel of NABA secretion (Figures 6C,D).

The KIRC population in TCGA-KIRC database was divided into immune-autophagy-related DEG low-expression group (G1) and immune-autophagy-related DEG high-expression group (G2), and the correlation between the expression of immune-infiltrating cells and immune-infiltrating cells was analyzed. The results show that CANX, NAMPT, BIRC5, and BID are highly correlated with the expression levels of a variety of immune infiltrating cells, and the expression levels of monocytes, myeloid dendritic cells, and CD8⁺ effector memory T cells are significantly correlated with CANX, NAMPT, BIRC5, and BID. It suggests that these cells may be related to the progression of KIRC (Figure 7).

Based on the original intention of this study to have a positive effect on the targeted drug therapy of KIRC, we also statistically analyzed the correlation between the expression level of immune checkpoints in KIRC tissues and the expression of immune-autophagy-related DEGs. The results showed that CD274, HAVCR2, LAG3, and PDCDILG2 were significantly correlated with BID (Figure 8A); CD274 and PDCDILG2 were significantly correlated with BIRC5 (Figure 8C); CTLA4, LAG3, PDCD1, PDCDILG2, TIGIT, and SIGLEC15 were significantly correlated with NAMPT (Figure 8E); and CD274, CTLA4, LAG3, PDCD1, and TIGIT were significantly correlated with CANX (Figure 8G). It can be seen that both CD274 and PDCDILG2 have appeared three times. It is speculated that they are sensitive immune checkpoints for KIRC treatment and diagnosis.

In addition, the response of CANX, BID, NAMPT, and BIRC5 with different expression levels to immune checkpoint inhibitors was predicted based on Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (Figures 8B–H). The results indicated that the p-values of all immune-autophagy-related genes except NAMPT were <0.05, which indicated that the immune checkpoint inhibitors were effective against KIRC with high expression of CANX, BID, and BIRC5, and the survival period was prolonged after immune checkpoint inhibitor treatment.

Following the same analysis method, the results in Figure 9 show that the expression of immune-autophagy-related DEGs is correlated with the expression levels of multiple ferroptosis-related genes, and NCOA4, EMC2, NFE2L2, HSPB1, SAT1, and DPP4 are significantly correlated with CANX, NAMPT, BIRC5, and BID.

In addition, we analyzed the correlation between m⁶A methylation-related genes and immune-autophagy-related DEGs by the same method and found that CANX, NAMPT, BIRC5, and BID were significantly correlated with multiple methylated genes (Figure 10). We further verified that m⁶A-related genes were differentially expressed in kidney cancer and normal tissues and were statistically significantly associated with patient prognosis (Supplementary Figures S1, S2). In particular, METTL14, VIRMA, ZC3H13, TYHDC2, YTHDF3, YTFDF2, IGF2BP2, and RBMX were significantly associated with four immune-autophagy-related DEGs.

By OCLR scores, we found that, except for BID, the expression levels of CANX, NAMPT, and BIRC5 were significantly different from the dryness degree of KIRC (Figure 11). These results suggested that CANX, NAMPT, and BIRC5 may influence the degree of similarity between KIRC cells and stem cells and thus affect the BP and degree of dedifferentiation of tumors.

To detect the expression of four genes (CANX, BID, NAMPT, and BIRC5) in KIRC, we performed the qRT-PCR in KIRC cells and clinical tissue samples. We verified the expression levels of four genes in normal kidney cell lines (HK-2 cells) and two KIRC cell lines (786-O and caki-1). The results showed that the expression levels of four genes were significantly increased in KIRC cells compared with normal kidney cells (Figures 12A–D). In addition, Western blotting results showed that protein levels of NAMPT and BIRC5 were expressed at increased levels in RCC cell lines 786 and caki-1, but there was no significant difference in protein levels of BID and CANX (Figure 12J). BID, NAMPT, and BIRC5 were detected with the same results in tumor tissues and with adjacent normal kidney tissues, while CANX was not significantly different (Figures 12E–H). Then we detected the protein expression of BID, NAMPT, and BIRC5 in the tissues by immunohistochemistry (IHC). Results demonstrated that NAMPT and BIRC5 were significantly increased in KIRC tissues compared with adjacent normal kidney tissues. However, BID was negative in most tissues (Figure 12I).