We systematically searched the transcriptome data of five autoimmune nephropathy from 36 datasets in the Gene Expression Omnibus (GEO) database (Edgar et al., 2002) based on the review articles (Kant et al., 2022; Boesen and Kakalij, 2021). These datasets included bulk RNA datasets of 5 Lupus Nephritis (LN) datasets (GSE112943, GSE200306, GSE32591, GSE81622, GSE99967), 7 IgA Nephropathy (IgAN) datasets (GSE175759, GSE37460, GSE104954, GSE141295, GSE116626, GSE35489, GSE115857), 10 Membranous Nephropathy (MN) datasets (GSE216841, GSE200828, GSE200818, GSE182380, GSE197307, GSE133288, GSE115857, GSE108113, GSE104954, GSE104948), 8 Focal Segmental Glomerulosclerosis (FSGS) datasets (GSE200828, GSE200818, GSE182380, GSE197307, GSE133288, GSE108113, GSE104954, GSE104948) and 3 Anti-neutrophil Cytoplasmic Antibody-associated Vasculitis (ANCA-AAV) datasets (GSE108113, GSE104954, GSE104948), as well as three scRNA datasets of IgAN (GSE171314), MN (GSE171458), and LN [GSE121893 (Arazi et al., 2019)] (Supplementary Data S1). It should be noted that we excluded cases and controls that were less than five in bulk datasets and those with obviously non-active LN in GSE99967. We also excluded one IgAN patient with microproteinuria in GSE17131.

For the immune datasets, we collected transcriptome and genomic data with clinical annotations from three ccRCC cohorts: IMmotion151 (Buttner et al., 2022) (atezolizumab + bevacizumab versus Sunitinib; first-line) from the European Genome-phenome Archive (EGA) dataset, CheckMate (nivolumab versus everolimus; latter-line) from the Supplementary Materials provided by Braun et al. (2020), and JAVELIN (avelumab + axitinib versus sunitinib; first-line) from the Supplementary Materials provided by Motzer et al. (2020a). We also gathered four clinically annotated transcriptome immune cohorts of urothelial carcinoma (UC) [IMvigor210 (Balar et al., 2017), treated by atezolizumab, n = 208], NSCLC [POPLAR (Fehrenbacher et al., 2016), treated by atezolizumab, n = 81], RCC [IMmotion150 (McDermott et al., 2018), treated by atezolizumab alone or in combination with bevacizumab, n = 162], multi-cancer [PCD4989g (Herbst et al., 2014), treated by atezolizumab, n = 206]. Patients with clearly four treatment outcomes of “CR, PR, PD, SD” were enrolled except JEVELIN which with known survival time records were enrolled.

Mutation and clinical data of ccRCC from The Cancer Genome Atlas (TCGA) database were downloaded from cbioportal (https://www.cbioportal.org/study/summary?id=kirc_tcga). This study generates no new data, so no ethics approval is needed. The Supplementary Data S1 of this study contains comprehensive information on all of the data.

All expression array data were log2-transformed using the GEO protocol. RNA-seq count data were normalized to FPKM for each gene based on the gencode v22 and log2-transformed, subsequently. The analysis of DEGs was performed using limma (R package, v3.2.3), with significant differences in gene expression in a single dataset being determined by p. value <0.05 and abs(logFC) > 1. Disease-associated genes were required to satisfy the requirement of being significant in three or more datasets, with the ANCA-AAV threshold set at two for only four datasets, exceptionally.

In the single-cell dataset, data were processed and analyzed using Seurat (R package, v3.2.3), and batch effect correction was performed using harmony (R package, v0.1.1). Significance thresholds in the single-cell analysis were set at p < 0.05 and avg_logFC > 0.25. The differential genes between immune-mediated kidney disorders and health from a single-cell perspective were determined to be significant either in the holistic condition or in more than two clusters.

Only four samples in the Bi dataset had information on the known response of ICB treatment, and the samples treated with atezolizumab + bevacizumab in IMmotion151 were mapped to the Bi dataset using Scissor (R package, v2.0.0) corresponding to the parameter’s alpha = 0.05 and cutoff = 0.25. For the results of the mapping only the cells with consistent response were retained, i.e., those whose outcome was the same as either response or non-response in the Bi. and in IMmotion151 cells were included in the downstream analysis. The mapped cells were still analyzed for DEGs between the response vs. non-response groups using Seruat, and significant DEGs were identified using p. value <0.05 and avg_logFC > 0.25.

All significant genes in each dataset were enriched pathways using the enrichPathway function from ReactomePA (R package, v1.30) and the profiling of 22 immune cells computed by CIBERSORT (Bulk RNA-seq) and SCINA (scRNA-seq), respectively. Autoimmune nephropathy-associated pathways need to meet the above criteria consistent with autoimmune nephropathy genes. Differential immune cells were screened in bulk using the Mann-Whitney U test and in scRNA-seq using the Fisher’ test. Immune cells associated with kidney disease also needed to meet the following criteria: (i) significant in three or more bulk datasets; (ii) significant in single-cell conditions; and (iii) significant in either bulk or scRNA for the same disease.

The log2FPKM matrices of IMmotion151, IMmotion150, JAVELIN, and PCD4989g were individually 0.75 quantile normalized, and combined with the CheckMate data, followed by batch correction using the ComBat function in the SVA (R package, v3.34). The expression matrix of genes for pre-feature selection was subset to autoimmune nephropathy-related genes. Pathway expression matrix analysis was performed using the ssGESA (GSVA, R package, v1.34) and subset to immune-mediated kidney disease pathways.

For the mutation matrix of all ICB datasets, only non-synonymous mutations were taken into account in this study. Candidate single-mutation variables that were significantly related to response (Fisher’ exact test, p. value <0.05) and survival (Cox-proportional hazards regression model, Wald test, p. value <0.05) in IMmotion151. Co-mutation genes are defined as pairs of genes that have either mutated or not. We filtered candidate co-mutation variables by response or survival in IMmotion151 and retained only those that were significantly associated with overall survival with KIRC in TCGA.

Five model variables require feature selection for subsequent analysis: genes expression, pathways expression, immune cell proportion, genes mutation, and co-mutation. Recursive feature reduction using caret (R package, v6.0-93) with 10-fold cross-validation was performed on all candidate variables. ML models were created using a variety of algorithms, including support vector machine (SVM), Naïve Bayes (NB), random forest (RF), k-nearest neighbors (KKNN), AdaBoost Classification Trees (AdaBoost), eXtreme Gradient Boosting (XGBoost), Neural Network (NeuralNet). The caret package (R package, v1.7.5.1) provides the first five methods, nerualnet package (R package, v1.44.2) provides the NeuralNet algorithm, and XGBoost obtain from the xgboost package (R package, v6.0-93).

To evaluate the performance of the models, we compared it with the ICB-derived gene model and biological functional characterization models sequentially. For ICB-derived genes selected by: (i) top 25 response significant genes (TopRes); (ii) top 25 survival significant genes (TopPFS); (iii) combined top 25 significant genes in response and survival (survival and response multiplied by p. value and sort in ascending order then took the top 25 genes, TopSig). The ICB-derived model built by XGBoost were identical to the TIs-ML model in procedure. The set of ICB-related functional gene lists were collected from multiple studies and is detailed. Expression scores for above functional signatures were calculated by the ssGSEA function from GSVA (R package, v1.34), and the pROC (R package, v1.18) was used to generate the ROC curves by calculate the sensitivity and specificity of the scores in predicting the ICBs’ binary responses. Other ICB response prediction model included that: the GEP score were computed by 18 immune genes expression which formula provided in Keynote-028 (Ott et al., 2019); the TME score were computed by TMEscore (R package, v0.1.4); TIDE score was computed by the TIDEpy python package (https://github.com/liulab-dfci/TIDEpy).

716 genes were identified in the immune-mediated kidney disorders samples from the above analytic processes and excluded the housekeeping genes. We clustered the expression matrix of autoimmune nephropathy genes of IMmotion151 using ConsensusClusterPlus (R package, v1.50), resulting in two distinct consensus that cluster by k-means and distance compute by Euclidean. The cosine similarity was determined by calculating the top 50 genes that were specifically highly expressed in each cluster. For the assignment of the CheckMate and JAVELIN sample clusters, greater cosine similarity with the known IMmotion151 clusters was labeled corresponding clusters (a similarity lower than 0.4 would assign in NA).

R v3.6.3 (https://www.r-project.org) was used to perform all statistical analyses. The two-sided Mann-Whitney U test was used to compare the difference in scores between the response and non-response (as well as other clinical) groups. Fisher’s exact test was utilized to determine whether there is a significant correlation between discrete variables and clinical groups. The logarithmic ranking test was used to assess prognostic differences between subgroups of clinical or model scores for Kaplan-Meier (KM) survival curves. Adjusted p. value of DEGs and GSEA analysis computed by the Benjamini-Hochberg method. By default, continuous variables greater than the median or quantile of 0.75 would be classified as high. The correlation between model score and immunological characteristics or cancer hallmarks was evaluated using Spearman correlation.

All immune-mediated disorders clinical and transcriptome data were collected from Gene Expression Omnibus (GEO) database (Edgar et al., 2002), except a scRNA data of LN, it’ processed data were from https://immunogenomics.io/ampsle, https://immunogenomics.io/cellbrowser/ and https://portals.broadinstitute.org/single_cell/study/amp-phase-1. For immuno-cohorts, clinical and sequencing data of IMmotion151, as well as IMvigor210, POPLAR, IMmotion150, and PCD4989g, were obtained from the EGA dataset with study IDs of EGAC00001001813 and EGAS00001004343, while CheckMate (Braun et al., 2020), JAVELIN (Motzer et al., 2020a) were from Supplementary Materials. Mutation and clinical data of ccRCC from The Cancer Genome Atlas (TCGA) database were downloaded from cbioportal (https://www.cbioportal.org/study/summary?id=kirc_tcga). Details see Supplementary Data S1.

Previously published works suggested a connection between the effectiveness of immunotherapy and autoimmune diseases (Grandi and Tramontano, 2018; Masetti et al., 2021). Figure 1 provides an overview of the conceptual workflow of our study. We adopted five autoimmune nephropathies (Kant et al., 2022; Boesen and Kakalij, 2021) with scRNA or RNA-seq data. Ultimately, we investigated a total of 42 cohorts, including eight lupus nephritis (LN) cohorts, nine IgA nephropathy (IgAN) cohorts, twelve membranous nephropathy (MN) cohorts, nine focal segmental glomerulosclerosis (FSGS) cohorts, and four antineutrophil cytoplasmic antibody vasculitis (ANCA-AAV) cohorts, which were over 1,900 samples. We obtained a total of five cohorts of ccRCC with ICB treatment, which included over 1,000 patients.

To investigate the association between inflammatory signatures of autoimmune nephropathy and immunotherapy, we analyzed the expression array data of patients with LN (GSE32591_glo) (Be et al., 2012). In total, we discovered 649 differential pathways and 30 clusters when comparing LN patients to healthy controls for differential expression analysis (Supplementary Figure S1). Macrophage, and interferon-related pathways were been identified, and those biomarkers have been shown to have a strong correlation with irAE (Jing et al., 2020). The results revealed that “regulation of leukocyte mediated immunity” was the most upregulated pathway class in LN, along with “regulation of viral genome replication” and “regulation of defense response to virus by host”, which are related to human endogenous retroviruses (HERVS). The results suggested that ccRCC was an inflammatory tumor, and inflammatory signals were connected to immunotherapy, which was consistent with previous studies (Hegde and Chen, 2020).

Inflammation is a complication of autoimmune nephropathy, and inflammatory signatures produced by the immune system are related to the immune response. To investigate the connection between inflammatory signals and immune responses, we disaggregated different types of bulk RNA data from over 1,900 patients with immune-mediated kidney disease. After conducting a differential analysis between the cases and healthy controls within each dataset of five diseases, we obtained 1,121 differentially expressed genes (DEGs) in LN 55, IgAN 90, MN 586, FSGS 559, and ANCA-AAV 581 separately (Supplementary Figure S2a). We found that the expression patterns of the top 30 DEGs were generally consistent across distinct datasets, indicating that certain shared molecular processes may contribute to the etiology of diverse autoimmune disorders (Figure 2a and Supplementary Data S2). Genes were enriched in the frequent pathways, such as “interferon gamma signaling” pathway, “interferon signaling” pathway, etc., that are relevant to irAE and immune response (Hu et al., 2023; Head et al., 2019; Ayers et al., 2017). Among them, the immunoregulation-related gene EGR1 (Early Growth Response 1), involved in Interferon alpha/beta signaling, which is a typical pro-inflammatory cytokine (Pozhitkov et al., 2017) and enhances anti-tumor capacity by regulating and activating T cells (Dan et al., 2020), was differentially expressed in all of the five diseases (Figures 2a,b). We further evaluated the aforementioned signatures in the IMmotion151 dataset (Buttner et al., 2022), where 400+ previously untreated ccRCC patients were given atezolizumab (an anti-PD-L1 inhibitor) plus bevacizumab (an anti-VEGF inhibitor), and found 329 genes relevant to response and survival (Figure 2c). Inflammation is a major class of functions enriched for differentially expressed genes in the autoimmune nephropathy patients compared to healthy controls, and such genes show the potential efficacy to predict response or prognosis of ICB therapy.

Furthermore, pathway signatures can replenish the predictions by gene signatures. Within bulk RNA data, we totally enriched a total of 73 differentially expressed Reactome pathways (DEPs, Supplementary Figures S2b,c and Supplementary Data S3). Among them, 33 pathways relevant to response and survival in IMmotion 151 (Figure 2d). In particular, we focused on the canonical pathways most associated with the inflammatory response: interferon pathway; Interleukin pathway; TLR pathway; NF-kB pathway; TCR pathway; BCR pathway, as well as the PD-1 pathway and Apoptosis pathway which are directly related to ICB therapy. The results showed that most of the inflammation-related pathways were significant in at least two or more diseases, especially the TCR pathway and apoptosis pathway were significant in all diseases, and the interleukin pathway and PD-1 pathway were significant in four diseases other than IgAN (Figures 2e-h). Intriguingly, most of these signatures upregulated in responders as compared to non-responders, and prediction ability of pathway signatures were slightly superior to those of gene signatures (Figures 2e-h). Pathways further explain the changes in inflammatory expression patterns of autoimmune nephropathy from a functional perspective and can serve as a credible complementary variable in ICB prediction.

Immune cells may be the direct hub of the immunoreactions in both immune diseases and tumors. We employed CIBERSORT to assess the relative proportions of 22 immune cell types in each specimen (Supplementary Figure S2d). One of the most significantly differentiated cell types, monocytes, was highly concentrated in immune diseases compared to controls. The infiltration and accumulation of them closely related to the pathogenesis of immune diseases (Zheng et al., 2020; Tang et al., 2021), and they can also differentiate into dendritic cells to present antigens, or secrete a variety of cytokines to active T cells during tumor immunotherapy (Ugel et al., 2021) (Supplementary Figure S2d). In summary, our findings suggest that different levels of inflammatory biomarkers derived from bulk RNA data of autoimmune nephropathy can provide valuable insights for predicting ICB outcomes.

Single-cell RNA sequencing can effectively complement bulk RNA sequencing through high-resolution expression data. We collected three single-cell RNA datasets (Supplementary Data S1) from LN, IgAN, and MN to perform a differential expression analysis. Specifically, we identified a total of 1,785 DEGs, which included 440, 1,264, and 616 genes in LN, IgAN, and MN, respectively (Supplementary Figure S3a). Compared with bulk RNA data, the scRNA-seq analysis reported 1,477 new DEGs, and only 308 (17.3%) DEGs were shared (Supplementary Data S4). This result confirmed that scRNA-seq was capable of complementing signatures at single cell levels. The DEGs identified in scRNA alone were mainly enriched in the “Immune System (HSA-168256)” and “Innate Immune System (HSA-168249)” (Supplementary Data S5). Moreover, “interleukin signaling” pathways associated with genes was upregulated, which turned out to be closely correlated with irAE (Jing et al., 2020) and was retrieved by scRNA-seq (Supplementary Figure S3d). In addition, a total of 518 DEGs were detected that were relevant to response and survival in the IMmotion151 dataset (Figure 3a). As a result, based on scRNA-seq data, more inflammatory genes involved in immunotherapy may be distinguished.

Furthermore, the results of DEPs through scRNA-seq data analysis were largely consistent. We additionally discovered 182 DEPs that were associated with autoimmune nephropathy, and 145 of these DEPs had inconsistencies with those discovered through bulk RNA datasets (Supplementary Figure S3b). DEPs, which were only remembered in scRNA data, were also crucial for immunoregulation and immunoreaction in malignancies. For example, “DDX58/IFIH1-mediated induction of interferon-alpha/beta (R-HSA-168928)” is an inflammation-related signal but is associated with immune response in tumor immunotherapy based on recognition of intracellular viral RNA and activation of interferon-stimulated gene expression, et al. (Yoneyama and Fujita, 2007; Loo et al., 2008; Honda et al., 2005) It was effective in distinguishing immune responders from non-responders (Supplementary Figure S3e). Above all, 66 pathway signatures were found using scRNA analysis in the IMmotion151 dataset, which was two times the number of those produced from bulk RNA data (Supplementary Figure S3c). Hence, scRNA-seq data enhanced the accuracy of results by revealing the differences resulting from the immune microenvironment.

Immune cells were annotated by SCINA in the scRNA-seq datasets (Zhang et al., 2019). We counted the immune cells within the 22 cell types that were associated with kidney diseases using Fisher’s exact test with p. value <0.05 (Figure 3b; Supplementary Figures S3f–h). T. cells.CD4. memory.resting was shown to have significantly lower levels in three disease groups comparing to healthy controls (Supplementary Figures S3f–h), which was consistent with previous reports (Okoye and Picker, 2013).

Eventually, we combined bulk and scRNA signatures by the principle: i) it must become significant in the relevant disease; ii) it in scRNA is needed to be significant generally or in at least two clusters; iii) it will be selected if significant in bulk or single-cell for each disease. As a result, we integrated 1,016 genes, 93 pathways, and 16 immune cells of inflammatory signatures from autoimmune nephropathy (Supplementary Figures S3i–k and Supplementary Datas S4, S5). In summary, our test emphasized the importance and necessity of adding single-cell RNA data to bulk RNA for inflammatory biomarkers exploring.

Obtaining immunopredictive markers in ccRCC patients undergoing immunotherapy may enhance the specificity of immune response predictions. Additionally, the confidence of the indicators can be improved by combining bulk and single-cell RNA data from immunotherapy. Therefore, we mapped the bulk RNA expression data of IMmotion151, divided according to responders (complete response and partial response; CR + PR) and non-responders (stable disease and progressive disease; SD + PD), to the scRNA-Seq data obtained from Bi et al. (2021), which included four patients treated with anti-PD1-based therapies (Figures 4a,b; Supplementary Figure S4a). This combination approach of bulk and single cells revealed a strong correlation within the responders between the two datasets (Figure 4c) as well as the two PR patients from Bi et al. (2021) (Supplementary Figure S4b; Pearson correlation = 0.84, p. value = 9.93e-05). Cells in Bi. that had consistent response status (responders or non-responders) with IMmotion151 were kept for further study. Additionally, we got 1711 DEGs and 240 DEPs involving the immune system, hemostasis, extracellular matrix organization, and signal transduction, between mapped effective and ineffective cells. (Figures 4d,e, and Supplementary Datas S4, S5). MT2A, a member of the metallothionein gene family, was obviously upregulated in responders but barely expressed in non-responders, which is consistent with previous findings (Deng et al., 2022), hinting at a potential immuno-predictive biomarker (Figure 4d; Supplementary Figure S4c). We further explored DICs using a similar approach to that used in kidney diseases (Figure 4f; Supplementary Figure S4d) and tested the prediction ability of these gene and pathway signatures within the IMmotion151 cohort (Supplementary Figures S4e–h). The aforementioned findings imply that the immunopredictive markers discovered using the mapping approach exhibit a strong correlation with the prognosis of ccRCC patients undergoing immunotherapy.

The mapped DEGs in Bi et al. were intersected with the autoimmune nephropathy DEGs of bulk and scRNA respectively. The genes common to these two groups were merged, and housekeeping genes were excluded. Protocol of DEPs process were same to DEGs. Finally, 716 genes and 92 pathways specific to immune-nephrosis were obtained. Incorporating these markers as intersections with those discovered in the context of renal immune diseases can effectively reduce the number of variables while enhancing the level of confidence. Eventually, 25 genes, 47 pathways, and four immune cells were reserved after feature selection by Recursive Feature Elimination (RFE) in the IMmotion151 training datasets (Supplementary Data S6).

Next, we selected seven ML models to compare AUC for multi-genes and multi-pathways individually. XGBoost obtained consistent predictions for response and prognosis by comparing several methods, which was retained for more in-depth study (Figure 4g; Supplementary Figures S4i–l). We created an ensemble model using XGBoost based on the models of genes, pathways, and immune cells. And the combined inflammatory signatures ML model could make a robust prediction for ICB response and progress-free survival (PFS) in training data from IMmotion151 (Motzer et al., 2022) (IM151), as well as from four additional independent cohorts: CheckMate (Braun et al., 2020) (CM), JAVELIN (Motzer et al., 2020a) (JV101), PCD4989g (Herbst et al., 2014) (PCD), and IMmotion150 (McDermott et al., 2018) (IM150). The inflammatory ensemble model was able to achieve an AUC of 0.98 in the predictive validation set IMmotion151 (Figure 4h), and the model’s scores were significantly able to differentiate PFS in IMmotion151 at p. value <2.18e-20 (Figure 4i). The performance of the inflammatory model continued to be noteworthy in the independent validation set, with AUCs greater than 0.6 and survival log-rank test significantly in both JAVELIN and IMmotion150 (Figure 4h; Supplementary Figure S4). In all datasets, groups with higher scores tended to have a better prognosis than those with lower scores (Supplementary Figures S4m–o).

To compare the model efficacy of ccRCC with ICB-derived models and immune-mediated kidney disorders-based models, we build (i) ORR (TopRes); (ii) PFS (TopPFS); (iii) doubly significant genes (TopSig) models (see Method for details). The results indicate that disease-derived (Dise_deriv) gene models perform more effectively and have higher robustness. However, the generalization ability of ICB-derived models’ performance was suboptimal (Figure 4j and Supplementary Data S7). In terms of prognosis prediction, the disease-derived model discriminates prognosis in three datasets, whereas the best ICB-based model, the doubly significant gene model (TopSig), can only discriminate prognosis in two datasets (Figure 4k).

To determine whether diseased-derived models demonstrate superior predictive validity over other biologically sourced models, we compared it to 11 published genetically relevant models (Supplementary Data S8). The results demonstrate that disease-derived models were still the most effective at predicting both response and prognosis, while other models were less effective or could solely predict prognosis but not response (Figures 4l,m). For instance, LRRC15. CAF.sig (Dominguez et al., 2020) could predict the prognosis among IMmotion151, JAVELIN, and IMmotion150, but response predictive efficacy AUC was all lower than 0.6. In conclusion, the model based on inflammatory biomarkers is obviously superior to immunotherapy-based signatures models or published features with tumor biology, both in distinguishing immune response and survival for ccRCC with ICB.

Immunogenic biomarkers are critical in assessing the immune response. However, the efficacy of using a single biomarker to predict the response to immunotherapy in ccRCC is constrained. The calculation of TMB, which includes all nonsynonymous mutations, demonstrated modest predictive capacity with an AUC of 0.572 in the IMmotion151 cohort (Supplementary Figures S5a,b). Using Fisher’s exact test and Cox proportional hazards regression analysis, we identified genes that were differentially mutated in response or prognosis with ICB. In the end, 14 mutation genes were selected as potential immunogenic signatures (Figures 5a,b).

Notably, the highest population frequency of the concurrent mutation gene pairs, EPHB6-GNAS, only in three patients (Supplementary Figure S5c). We calculated co-mutation genes, which are pairs of genes that either one of them were mutated. Hence, we screened for co-mutated gene pairs in the IMmotion151 cohort by analyzing the differences in immune response and prognosis. In total, we identified 1,876 combinations of gene pairs (Figures 5c,d). To reduce selection bias, we filtered additional co-mutation genes based on the overall survival of The Cancer Genome Atlas (TCGA) in KIRC. Ultimately, we discovered 49 co-mutation genes associated with response and survival in patients.

We compared seven different approaches to building machine learning models based on multigene and co-mutation genes. XGBoost consistently demonstrated satisfactory performance in predicting response and prognosis using the IMmotion151 or CheckMate cohorts (Figures 5e,f; Supplementary Figures S5d–g). Overall, predictive models based on the three individual types of immunogenic signatures have demonstrated limited ability to estimate immune response and survival (Figures 5e,f; Supplementary Figures S5d–g). Finally, we established a machine-learning model that integrates all the immunogenic biomarkers. The model achieved an AUC of 0.632 and 0.609 in the testing set IMmotion151 and the independent validation set CheckMate cohort, respectively (Figure 5g). The comprehensive model’s predictive score was capable of stratifying the prognosis of different cohorts (Figures 5h–j). These results revealed that the single immunogenic model was inadequate, and the multi-indicator combination model improved the predictive efficiency of immunotherapy.

The inflammatory model was more efficient than immunogenic model, but there’s still room for improvement in the inflammatory model. Two models derived from the distinct perspective and integration of the two may further enhance the predictive efficiency. We aimed to develop a machine learning model using multi-omics data to predict response in ccRCC patients receiving ICB treatment. The multi-omics signatures landscape were shown in Figure 6a, Supplementary Figures S6a,b, Supplementary Datas S9, S10. The ML-based multi-omics model was constructed using 354 patients with ICB from IMmotion151 who had RNA-seq and DNA data. The patients were randomly assigned into a training and validation set with a ratio of 2:1. Our TIs-ML model achieved an excellent AUC of 0.997 for distinguishing responders from non-responders in the validation set, significantly surpassing the predictive power of GEP score (AUC: 0.627), TME score (AUC: 0.574), TIDE score (AUC: 0.560), TMB (AUC: 0.572), and PD-L1 expression (AUC: 0.555) (Figure 6b; Supplementary Figure S5b). Positive findings were consistently observed when the model was applied to the CheckMate cohort; our model had an AUC of 0.705, whereas the GEP score had 0.498, the TME score had 0.509, and the TIDE score had 0.573 (Figure 6c). The present study suggests that the TIs-based predictive signatures are highly generalizable and that the TIs-ML model can robustly predict response to ICB.

In addition, individuals in the high-score group (with the threshold set at the median) in the IMmotion151 cohort had a longer PFS than those in the low-score group (Figure 6f; p. value < 0.0001). The CheckMate and JAVELIN101 cohorts (Figures 6g,h; Supplementary Figure S6c) showed a better prognosis in the high-score group (p. value = 0.027 for PFS and 0.54 for overall survival in the CheckMate cohort; p. value = 0.010 in the JAVELIN101 cohort). Similarly, The TIs-ML score outperformed typical biomarkers for predicting prognosis (Supplementary Figures S6d,e). These findings demonstrated that the TIs score could act as a prognostic biomarker for better patient stratification in ccRCC with ICB treatment.

More importantly, we revealed that inflammatory signatures demonstrated a higher importance coefficient than immunogenic signatures (Figure 6e) in the multi-omics model. Autoimmune nephropathy-related pathways and genes played a dominant role in the model (Figure 6d), suggesting that there was a strong correlation between immune-mediated kidney disease and immunotherapy.

We examined the correlation between the TIs-ML scores and clinical factors after observing a significant connection between the risk score and tumor regression (Figure 7a). A significant difference in TIs score was observed between sarcomatoid and non-sarcomatoid ccRCC, with sarcomatoid ccRCC tending to have a higher TIs score (Figures 7e,f). Our study also identified that the Memorial Sloan-Kettering Cancer Center (MSKCC) (Motzer et al., 1999) and the International Metastatic Renal Cell Carcinoma Database Consortium (IMDC) (Heng et al., 2009; Ko et al., 2015), which are commonly used to predict the prognosis in metastatic RCC using serological biomarkers, were less effective compared with TIs score (Supplementary Figures S7a,b,d,e). The results demonstrated that there was no significant difference in TIs score between primary and metastatic lesions (Supplementary Figures S7c,f), which emphasizes the robustness of our model in predicting response to ICB.

Furthermore, we explored the correlation between TIs scores and well-known immunotherapy biomarkers, including TMB and PD-L1. Responders were more closely related to high scores than high TMB (Figure 7c) and positive PD-L1 (Figures 7b,d). Besides, immunotherapy responders and survival-benefit patients with TMB-low or PD-L1-negative could be bailed by higher TIs-based risk scores, emphasizing that our multi-omics model can complement the prediction ability of conventional biomarkers (Supplementary Figures S7i,j).

We examined the association between risk scores and published ICB-related signatures, and identified a positive correlation between the risk scores and a variety of immune-related hallmark genes. The inflammatory response enhances the efficiency of immune checkpoint therapy by increasing immune cell infiltration and improving tumor immunogenicity (Supplementary Figures S7g,h). The study showed that the TIs score had a positive correlation with activated CD8 T cells and antigen-presenting signatures, as well as novel gene expression signatures such as NK cells, Treg, and B cells in the analysis of the TME (Figures 7g,h; Supplementary Figure S7h). These immune cell signatures are known to play a key role in the tumor-killing process.

Since the model score can predict the response and prognosis of immunotherapy in ccRCC, we tested the specificity of our model by evaluating non-immunotherapy treatment for ccRCC, and also by assessing immunotherapy treatment for other cancer patients. The risk score was ineffective in predicting response or prognosis between the single-agent targeted therapy and antiangiogenic therapy in ccRCC from the IMmotion151, JAVELIN, and CheckMate cohorts (Figures 8a–d). The IMvigor210 (mUC) (Balar et al., 2017) and Poplar (NSCLC) (Fehrenbacher et al., 2016) datasets were also analyzed using the ML-based inflammatory model with negative results for predicting prognosis or response (Figures 8e–h; Supplementary Figure S8). In conclusion, our method can specifically predict ICB outcomes in ccRCC patients, but is not suitable for predicting efficacy in ccRCC patients receiving other therapies or for individuals with other types of cancer undergoing immunotherapy.

To further explore the roles and mechanisms of signatures derived from immune-mediated kidney disorders in guidance ICBs for ccRCC, biological function pathway enrichment analysis was performed for the disease-derived genes.

Pathway enrichment analysis indicated that immune-mediated kidney disease genes have a strong association with biological processes, which plays a crucial role in immune stress and evasion in tumors (Figures 2e, 3c; Supplementary Figure S9a). Furthermore, we performed a consensus clustering analysis using these genes in the IMmotion151 dataset and identified two distinct clusters (Figure 9a). C1 and C2 differentially expressed genes were mainly enriched seven pathways, including “regulation of lymphocyte activation” and “nuclear division,” which plays a significant role in immunogenic cell death and exhibits high expression levels in C1 compared to C2 (Figure 9b). Interestingly, patients in Cluster 1 had a poorer prognosis. However, Immunotherapy improved PFS in the C1 patients while in the C2, neither immunotherapy nor any other therapy had an effect on the prognosis (Figure 9d; Supplementary Figures S9b–d). Consensus cluster analysis suggests that cluster1 patients are more suitable for monotherapy or combination therapy of ICB, while cluster2 could select more economical therapy in ccRCC. However, immune kidney disease signatures are likely only applicable to KIRC and may not be generalizable to most other cancer types in pan-cancer analysis (Figures 9c,e,f; Supplementary Figures S9e–h). In summary, the results demonstrate that kidney disease-related genes have a very strong correlation with immunotherapy at the molecular level and can specifically predict immunotherapy in renal cancer.