Five public AKI datasets from the Gene Expression Omnibus (GEO) were downloaded and processed: GSE126805 (n = 83; normal tissues = 41, AKI tissues = 42), GSE139061 (n = 48; normal tissues = 9, AKI tissues = 39), GSE30718 (n = 47; normal tissues = 19, AKI tissues = 28), GSE90861 (n = 46; normal tissues = 23, AKI tissues = 23), and GSE43974 (n = 391; normal tissues = 188, AKI tissues = 203). The former four datasets were integrated and the batch was removed to construct a training cohort, whereas GSE43974 was treated as an independent test cohort. GSE126805, GSE139061, GSE30718, and GSE90861 were used to test whether these biomarkers could distinguish AKI tissues from normal tissues. Approval and informed consent from the institutional review board were not required for the AKI cohorts from the public databases. To reveal the connection between AKI and ccRCC, ccRCC multi-omics information was retrieved from the GDC TCGA portal.

To remove batch effects derived from the study design, sequence platform, and technological replication, we filtered only normal and AKI tissue expression matrices and clinical characteristics from four cohorts: GSE126805, GSE139061, GSE30718, and GSE90861. The batch effect was removed using the default function in the sva package. A principal component analysis (PCA) was used to visualize the efficacy of the batch removal.

We used the limma package to identify differentially expressed genes (DEGs) between normal and AKI tissues in the merged expression matrix. DEGs were identified using a significance threshold of p < 0.05 and an absolute log-fold change >1.5. An enrichment analysis of the DEGs was performed using the clusterProfiler and ggpplot2 packages, and these were visualized using an enrichment plot. Additionally, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene set enrichment analyses (GSEA) were performed to better understand the biological roles of these DEGs. The significantly differentially expressed terms or pathways were identified by a threshold value of p <0.05 and q-value <0.05.

After identifying the DEGs from the combined expression profile of the combined expression matrix, we sought to identify the most relevant AKI-related biomarkers. We adopted four machine-learning algorithms, including Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression, Random Forest (RF), eXtreme Gradient Boosting (XGBoost), and support vector machine (SVM), to select reliable biomarkers to distinguish AKI from normal tissues. The detailed parameters adopted for these were as recommended in each of the R packages. The detailed parameters of LASSO logistics were as follow: alpha = 1, maximum number of iterations = 8,000, tol = 1e-4). The detailed parameters of Random Forest were as follow: n_estimators = 80, criterion = gini, min_samples_split = 2, min_samples_leaf = 1. And for XGBoost, the detailed hyperparameters were as follow: n_estimators = 500, max_depth = 6, reg_alpha = 0, colsample_bytree = 1. Finally, the detailed parameters of SVM were as follow: tol = 1e-3, max_iter = −1. Furthermore, the receiver operating characteristic (ROC) curve was adopted to further evaluate the accuracy of the filtered biomarkers in the training and testing cohorts using the R package pROC. The GSE90861 dataset was used as a validation cohort to verify the specificity and sensitivity of the biomarkers retrieved from the test cohort.

Single-sample gene set enrichment analysis (ssGSEA) and bulk sequence-based deconvolution algorithms from the GSVA and ESTIMATE packages were used to detect different immune cells and components between normal and AKI tissues. The gene sets used for ssGSEA consisted of 28 types of immune cells, including activated B cells, activated CD4 T cells, activated CD8 T cells, activated dendritic cells, CD56^bright natural killer cells, CD56^dim natural killer cells, central memory CD4 T cells, central memory CD8 T cells, effector memory CD4 T cells, effector memory CD8 T cells, eosinophils, gamma delta T cells, immature B cells, immature dendritic cells, macrophages, mast cells, MDSCs, memory B cells, monocytes, natural killer cells, natural killer T cells, neutrophils, plasmacytoid dendritic cells, regulatory T cells, T follicular helper cells, type 1 T helper cells, type 17 T helper cells, and type 2 T helper cells. The estimated algorithms contained three scores: ESTIMATE, immune, and stromal scores. The correlations between the seven biomarkers and immune cells or their components were calculated using the R package ggcor.

To further reveal the heterogeneity of AKI and the potential relationship between AKI and ccRCC, we performed an unsupervised cluster analysis based on seven novel biomarkers (CCNL1, NFKBIZ, HBB, TRIB1, SOCS3, HSPA6, and EGR1) using the Consensus Cluster Plus package. The optimal cluster number was identified based on the PCA algorithm and cumulative distribution function curve. Differences in estimates, immune scores, and immune components were also qualified and compared between the ccRCC subtypes. A detailed description of the parameters employed are described in previous studies (15, 16).

Figure 1A shows a flowchart of our study. Figure 1B shows the raw PCA results for the four microarray databases (GSE126805, GSE139061, GSE30718, and GSE90861). The four different colors represent different datasets. Every dataset was discrete without any intersection. Figure 1C shows the combat PCA results following batch removal.

A volcano map (Figure 2A) was created to visually display the expression changes of the DEGs between normal and AKI tissues. A total of 108 genes were identified as DEGs according to adjusted p-values <0.05 and logFC cutoffs >1.5; among these, 101 genes were upregulated and 7 genes were downregulated. Next, we performed GO and KEGG enrichment analyses. As depicted in Figures 2B–D, in terms of biological processes (BP), the DEGs were mostly associated with lipopolysaccharides and molecules of bacterial origin. Regarding cellular components (CC), these were mainly involved in the RNA polymerase II transcription regulator complex and transcription regulator complex. In terms of molecular functions (MF), the DEGs were largely associated with tyrosine proteins, threonine phosphatase activity, and MAP kinases. The 20 most significant KEGG pathway terms are shown in Figure 2E. These genes were primarily associated with the TNF signaling pathway, the IL-17 signaling pathway, lipid and atherosclerosis pathways, the AGE − RAGE signaling pathway in those with diabetic complications, and the MAPK signaling pathway.

GSEA was performed to determine which KEGG pathways were differentially enriched between the AKI and normal groups according to the NES and value of p (p < 0.01) criteria. As presented in Figures 2F–K, the top three significantly upregulated KEGG pathways enriched in the AKI group were pathways related to cytokine-cytokine receptor interactions, IL-17 signaling, and TNF signaling. In addition, oxidative phosphorylation, fatty acid degradation, and peroxisomes were significantly downregulated in the AKI tissues. These data reveal that numerous pathways may be directly or indirectly involved in the regulation of AKI.

Four machine learning algorithms were applied to identify feature genes: LASSO regression (Figures 3A,B), SVM (Figure 3C), and Random Forest (Figures 3D,E). These were combined with a feature selection analysis to determine the top 30 genes with relative importance, and an XGBoost analysis was then performed using these 30 relatively important genes (Figure 3F). A Venn diagram was used to display the overlapping signatures of the four methods (Figure 4G). After applying the four machine learning algorithms, eight signatures were selected, as shown in Figures 3A–G: tribbles homolog 1 (TRIB1), ST6GALNAC3, suppressors of cytokine signaling 3 (SOCS3), NF-κB inhibitor ζ (NFKBIZ), heat shock protein family A 6 (HSPA6), hemoglobin subunit beta (HBB), early growth response 1 (EGR1), and cyclin L1 (CCNL1).

As shown in Figures 3H,I, to further test the diagnostic efficacy of the above genes, ROC analysis was performed to evaluate the accuracy of each gene. The ROC curve showed that the seven biomarkers identified based on the logistic regression were reliable, with the following AUCs in the training set: CCNL1, 0.972; NFKBIZ, 0.937; HBB, 0.726; TRIB1, 0.908; SOCS3, 0.928; HSPA6, 0.905; and EGR1, 0.981. In the validation set, the AUCs were as follows: CCNL1, 0.825; NFKBIZ, 0.868; HBB, 0.821; TRIB1, 0.817; SOCS3, 0.826; HSPA6, 0.781; and EGR1, 0.814. ST6GALNAC3 (AUC = 0.517 in the test set) was not analyzed further in the next workflow.

We analyzed the heterogeneous compositions of the tumor microenvironments (TMEs) in the AKI and normal groups. As depicted in Figures 4A–C, all estimated scores, including the immune, stromal, and estimated scores, were lower in the AKI group compared to the normal group. We then examined the specific correlations between the estimated scores and each of the signatures using Spearman’s correlation analyses, which revealed that TRIB1, EGR1, and CCNL1 were related to the estimate and stromal scores (Figure 4D). As shown in Figure 4E, almost all of these genes were positively or negatively correlated with immune cell infiltration. Of these, EGR1 and eosinophil, HSPA6 and mast cells were significantly and positively relevant. In contrast, TRIB1 and activated B cells, CCNL1 and central memory CD8 T cells were negatively correlated (p < 0.01).

As shown in Figure 4F, the hub gene correlation analysis revealed that the expression levels of CCNL1 and NFKBIZ were mostly positively correlated, and CCNL1 was linked to every other gene. Figure 4G shows that 28 immune cells were differentially concentrated between the AKI and control groups. Among these, activated CD4 T cells, CD56^dim natural killer cells, eosinophils, mast cells, memory B cells, natural killer T cells, neutrophils, T follicular helper cells, and type 1 T helper cells were more enriched in the AKI cluster.

In order to accurately predict AKI risk using these signatures, a nomogram was constructed. As shown in Figure 5A, in the nomogram, each biomarker gene had a parallel score within the range of 0–100 points in terms of its association with the risk of AKI. As shown in Figures 5B–E, the AUC value was 0.919 in the training set and 0.945 in the testing set. In both the training and validation sets, the calibration plot demonstrated a good fit of the constructed nomogram, indicating satisfactory model accuracy.

We investigated the expression patterns of AKI-related signatures in ccRCC to comprehensively explore the association between AKI and cancer. TCGA-ccRCC samples were classified into different molecular subtypes according to AKI signature expression levels. For this, we divided the optimal cluster numbers using an unsupervised clustering method. As depicted in Figures 6A–C, considering the relative change in area under cumulative distribution function curve and consensus index, we chose k = 2 as the optimal cluster number. The TCGA ccRCC dataset was classified into two optimal clusters: CS1 and CS2. Compared to CS2 patients, CS1 patients had better overall survival (OS) (p = 0.0026, log-rank test) and progression-free survival (PFS) (p = 0.00033, log-rank test) (Figures 6D,E). In addition, we analyzed the expression levels of AKI genes in normal tissues and in the two ccRCC subtypes using a heatmap. CCNL1, NFKBIZ, SOCS3, and HSPA6 were significantly upregulated in the ccRCC subtype compared to normal tissues, whereas TRIB1 and EGR1 were significantly upregulated in the ccRCC subtypes (p < 0.05; Figure 6F).

Considering the varying characteristics of each cluster, we identified the gene expression profiles of CS1 and CS2. Supplementary Figure S1A presents an enhanced volcano map that displays the differential gene expression between CS1 and CS2, with 19,168 total differentially expressed genes. In addition, the GO analysis demonstrated that DEGs between ccRCC subtypes were annotated in receptor antagonist, receptor inhibitor, peptidase inhibitor, peptidase regulator, and enzyme inhibitor activity in biological process (Supplementary Figure S1B).

Supplementary Figure S1C shows the results of the GSVA algorithm analysis, which demonstrates that G2M-checkpoint, TGFβ signaling, and E2F oncogenes were more enriched in CS2 than in CS1. In addition, regulon analysis was used to detect transcriptome differences, which demonstrated that HNF1A, EPAS1, ZEB2, TFE3, and TP53 were downregulated in CS2, whereas FOXE1 and TBX18 were upregulated in CS1 (Supplementary Figure S1D).

To further elucidate the molecular distinctions between the two ccRCC subtypes, we conducted a comprehensive analysis of the immune-related signatures, signaling pathways, anticancer responses, and other relevant factors in both subtypes. We found that the dysfunction scores and Tumor Immune Dysfunction and Exclusion (TIDE) scores in CS1 were significantly lower than those in CS2 (Figures 7A,B), which indicated lower tumor immune evasion and resistance to checkpoint blockage in this cluster. In addition, the expression levels of nine immune checkpoint inhibitor genes were compared between the CS1 and CS2 cells; among them, CD274 was significantly upregulated in CS1 cells, whereas LAG3 and PDCD1 were significantly downregulated in CS2 cells (Figure 7C).

Among the different immune-related signatures expressed between the two ccRCC subtypes in the heatmap (Supplementary Figure S2A), chemokines, chemokine receptors, major histocompatibility complexes, and immunoinhibitory and immunostimulatory signatures were significantly differentially expressed. The CS1 subtype showed lower expression of CXCL11, CCR4, HLA-DOA, CD244, and IL2RA. Moreover, CS1 showed higher immune cell infiltration than CS2 in the TME infiltration cell-type heatmap (Supplementary Figure S2B). Mismatch repair, nucleotide excision repair, and base excision repair signatures were significantly higher in CS1, whereas EMT1, EMT2, signatures were significantly lower this cluster (Figure 7D).

As the anticancer immune response can be interpreted as a series of immune cell combat processes, we evaluated the tumor immunophenotypes in the two clusters. As shown in Figure 7E, the CS2 cluster demonstrated higher activity associated with priming and activation (step3), T cell recruiting (step4), dendritic cell recruiting (step4), macrophage recruiting (step4), B cell recruiting (step4), Th2 cell recruiting (step4) and recognition of cancer cells by T cell (step6). In contrast, CS1 demonstrated higher activity associated with the release of cancer antigen (step1) and cancer antigen presentation activity (step2) (p < 0.05).

To evaluate the drug responses of the ccRCC subtypes, the estimated IC50 data for each drug were collected from the GDSC database to explore the potential chemotherapy sensitivity of ccRCC. We found that the CS1 subtype was more sensitive to sunitinib, pazopanib, crizotinib, erlotinib, temsirolimus, and axitinib, whereas the CS2 subtype was more sensitive to lisitinib and gefitinib (Figure 8A). A significant difference was also observed in other omics. For instance, the tumor mutation burden was higher in CS2 (Figure 8B), which may have led to poor prognoses in CS2. In addition, Figure 8C shows that CS2 had higher amounts of copy number alterations, copy number losses, and copy number-gained genomes compared to CS1 (p < 0.001).

To further confirm the credibility of the ccRCC subtypes, the nearest template prediction algorithms were applied to reconstruct the subtypes in the TCGA–KIRC cohort and divide this cohort into two different subgroups. As shown in Figures 9A–D, the four Kaplan–Meier curves created using the data from different databases revealed that CS1 had a significantly better survival probability than CS2. In the Cancer Cell database, the value of p was 0.012; in the GSE22541 database, the value of p was 0.038; in the CheckMate OS database, the value of p was 0.015; and in the E-MTAB-3267 database, the value of p was 0.046. The stability and reliability of the subtypes were verified based on these results.