Construction and validation of a model based on immunogenic cell death-associated lncRNAs to predict prognosis and direct therapy for kidney renal clear cell carcinoma

Data sources and access

The University of California, Santa Cruz (UCSC) Xena database (https://xena.ucsc.edu/, until November 1st, 2022) was the primary source of data for this study [9]. The Cancer Genome Atlas (TCGA) database in Xena provided patient information, including transcriptomic data for 607 KIRC cases and clinicopathological information for 979 KIRC clinical patients. Tumor somatic cell mutation data were obtained from the TCGA database (https://portal.gdc.cancer.gov/repository, until November 1st, 2022). Transcriptomic data were obtained in two data formats, HTSeq-Counts and HTSeq-FPKM, which are used for different types of data analysis. We selected patients with complete transcriptional and clinical information and ensured that they were identified in both datasets. Ultimately, a total of 597 samples were analyzed in this study.

Select ICD-associated genes and ICD-related lncRNAs

ICD-associated genes were identified from the GeneCards website (https://www.genecards.org/, until November 1st, 2022). We conducted a search using immunogenic cell death as a keyword, and target genes were selected for subsequent analysis [10]. In addition, we conducted Pearson correlation analysis to investigate the associations among ICD-associated genes and all lncRNAs and acquired ICD-related lncRNAs (filter conditions indicated as |correlation coefficient| > 0.3 and p < 0.05). We ran variance analysis on the results obtained from Pearson correlation analysis using the “DESeq2” package [11], with the retention of lncRNAs with a p value < 0.05 and |log2-fold change > 3|.

All KIRC patients were randomly assigned 1:1 to the training cohort, which was used for model construction, or the test cohort, which was used to verify the model. Then, ICD-associated lncRNAs relevant to patient survival were short-listed in the training cohort using univariate Cox regression analysis. We used least absolute shrinkage and selection operator (LASSO) regression analysis to avoid overfitting. Finally, ICD-related lncRNAs were identified through multivariate Cox regression analysis.

KIRC ICD-associated lncRNA prognostic model

We used the results from predictive model building. The format of the risk score was as follows: risk score = coefficient (lncRNA₁) × expression (lncRNA₁) + coefficient (lncRNA₂) × expression (lncRNA₂) + coefficient (lncRNA₃) × expression (lncRNA₃) +......+ coefficient (lncRNA_n) × expression (lncRNA_n). Using the midpoint of the risk score as its threshold, the training group, the test group and all patients were sorted into high- and low-risk cohorts for survival analysis. Kaplan-Meier (K-M) survival analysis utilizing “survival” and “survminer” was conducted to examine the variation in overall survival (OS) between risk categories in the training cohort, test cohort, and total cohort. Receiver operating characteristic (ROC) values at one, three, and five years were extrapolated to evaluate the predictive performance of the signature.

Principal component analysis (PCA) using the expression profiles of all genes, all lncRNAs, lncRNAs associated with ICD, and lncRNAs in the selected prognostic models was conducted to validate the subgrouping effect. In addition, univariate and multifactor Cox independent prognostic analyses were conducted to assess risk scores and clinical data to validate the predictive value of the risk model. After that, K-M survival analysis was conducted to probe the model feasibility under dissimilar clinical characteristics.

Creation of a nomogram for predicting patient survival

We used the R programming language “RMS” to create a prognostic nomogram using age, risk group, and tumor stage as factors to predict patient outcomes at 1, 3, and 5 years [12]. We also determined the value of this nomogram in prognosis prediction.

Enrichment analysis

KEGG enrichment analysis was applied to optimize the functions that might be interrelated with the genetic set. We divided the cases into high- and low-risk groups under the simulations constructed from 8 ICD-related lncRNAs, and enrichment analysis was used to filter the pathways with a p value<0.05 and false discovery rate (FDR)<0.25 using GSEA_4.3.2 software (downloaded at https://www.gsea-msigdb.org) [13].

Tumor mutational burden contributes to the prognostic evaluation of tumors

We obtained tumor somatic mutation data from the TCGA database. We manipulated the data using the “TCGAbiolinks” package, sketched waterfall plots using the R programming language “maftools”, and calculated TMB [14, 15].

Evaluation of TME and immuno-infiltration

The R software package “ESTIMATE” was used as an analysis vehicle to calculate any difference in stromal score, ESTIMATE score, immune score, and tumor purity between the high- and low-risk subgroups. We obtained further information through the Tumor Immune Estimation Resource (TIMER) 2.0 website (http://timer.cistrome.org/, until November 1st, 2022) using “TIMER”, “CIBERSORT”, “CIBERSORT-ABS”, “QUANTISEQ”, “XCELL”, and “EPIC”, together with “MCPCOUNTER”, which are six methods used to optimize the correlation between individual immune cell types and risk scores [16]. To detail the constitutive shifts in immune cells, we undertook an immune infiltration analysis of 22 immune cell lineages.

ICD-related signature in immunotherapy versus chemotherapy

To fully address the distinctions between KIRC patients in different risk groups presenting with tumor immune dysfunction and rejection, the Tumor Immune Dysfunction and Exclusion (TIDE) website (http://tide.dfci.harvard.edu/, until November 1st, 2022) was used to retrieve the TIDE scores and related data of KIRC patients [17]. Therefore, we conducted single-sample gene set enrichment analysis (ssGSEA) using the “GSVA” R package. Half-maximal inhibitory concentrations (IC50) of typical antineoplastic drugs were determined using the R software “oncoPredict” package [18]. The IC50 values were compared between different classes utilizing a Wilcoxon signed-rank test.

Validation by RT-qPCR

Twenty-one KIRC samples were obtained from the Second Affiliated Hospital of Nanchang University, and patients were divided into high- and low-risk groups. The study was approved by the Ethics Committee of the Second Affiliated Hospital of Nanchang University, and all participants gave their informed consent. We extracted RNA from each KIRC tissue sample using TRIzol reagent (Life Technologies CA, USA) and randomly selected samples for RT-qPCR analysis. The experiments were performed using BlazeTaq SYBR Green qPCR master mix (GeneCopoeia, Guangzhou, China) and the Applied Biosystems 7500 Fast Real-Time PCR System (Applied Biosystems). All RNAs of every sample were analyzed in three independent experiments. The primers for ICD-associated lncRNAs are shown in Supplementary Table 1. The relative expression of lncRNAs was calculated using the 2^(-ΔΔCt) method.

We used the Human Protein Atlas (https://www.proteinatlas.org/) database to compare the protein expression levels of selected ICD-related genes in KIRC tissues with those in normal tissues.

Availability of data and material

The data sets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Screening for ICD-associated lncRNAs in KIRC patients

The overall study design of this paper is displayed in a flow chart (Figure 1). First, we identified ICD-related genes and constructed a protein–protein interaction (PPI) network associated with 16 of these genes by using the STRING website (Figure 2A). Further Pearson correlation analysis was executed, from which we derived 9096 ICD-associated lncRNAs. Next, we narrowed the pool down to 547 lncRNAs that were shown by differential analysis to be differentially expressed in tumor tissues (Figure 2B).

Based on difference analysis, we rerandomized 526 patients into training and test groups (Table 1). In the training group, a single-factor Cox model was used to screen 180 lncRNAs associated with differential prognosis using p<0.05 as the cutoff value (Supplementary Table 2). Follow-up LASSO analysis revealed 14 differentially expressed ICD-related lncRNAs (Figure 2C, 2D) (Supplementary Table 3), and multifactorial Cox regression analysis was used to further screen the 8 most significant lncRNAs (Supplementary Figure 1A). We subsequently explored the correlations among these 8 lncRNAs (Supplementary Figure 1B) and their correlation with ICD-related genes (Supplementary Figure 1C) that were found to exhibit a marked difference in up- and downregulation in tumor tissues (Supplementary Figure 1D).

Building and validating the prognostic model

Based on the 8 ICD-related lncRNAs, we calculated risk scores for different patients: risk score= coefficient(AP000439.3)×expression (AP000439.3)+ coefficient (RP11.1151B14.5)×expression(RP11.1151B14.5)+coefficient(RP11.479J7.2)×expression(RP11.479J7.2)+coefficient(AC099552.4)×expression(AC099552.4)+coefficient(RP11.19E11.1)×expression(RP11.19E11.1)+coefficient(CTB.33O18.1)×expression(CTB.33O18.1)+coefficient(RP11.339D23.1)×expression(RP11.339D23.1)+coefficient(LINC01192)×expression(LINC01192). Clinical information for all patients in the high- and low-risk clusters is presented in Table 2. In the training cohort, the test cohort, and the total cohort, survival was much worse in the high-risk group than in the low-risk group (all p values<0.001) (Figure 3A–3C), and the time-dependent ROC curve demonstrated that the model’s predictive value was excellent (AUCs were 0.831, 0.775, and 0.796 for the training cohort at one, three, and five years, respectively; 0.710, 0.688, and 0.768 for the test cohort at one, three, and five years, respectively; and 0.775, 0.743, and 0.786 for the total cohort at one, three, and five years, respectively) (Figure 3D–3F). In different cohorts, the manifestations of 8 lncRNAs in patients with different risks (Supplementary Figure 2A–2C), the risk curves (Supplementary Figure 2D–2F), the risk distribution plot (Supplementary Figure 2G–2I), and the scatter plot (Supplementary Figure 2J–2L) showed the forecasting value of the model.

The outcomes of principal component analysis (PCA) (Figure 4A–4D) indicated that the high- and low-risk subgroups showed obvious categorical clustering. Univariate Cox independent prognostic analysis indicated that age (HR=1.68, 95% CL=1.24-2.28, p<0.001), tumor stage (HR=1.92, 95% CL=1.68-2.20, p<0.001), and risk score (HR=1.44, 95% CL=1.34-1.54, p<0.001) were standalone risk factors (Figure 4E). Likewise, multifactorial independent prognostic analysis showed that age (HR=1.57, 95% CL=1.15-2.15, p=0.004), tumor stage (HR=1.73, 95% CL=1.50-1.99, p<0.001), and risk score (HR=1.34, 95% CL=1.25-1.44, p<0.001) were standalone risk factors (Figure 4F).

The expression heatmap containing different clinical data, risk groupings, and prognostic model-related lncRNAs (Supplementary Figure 3A) and the survival curves of patients in a variety of clinical states (Figure 5) illustrate the significance of the model predictions (stage I-II: p<0.001, stage III-IV: p<0.001; female: p<0.001, male: p<0.001; age<=65: p<0.001, age>65: p=0.002).

Clinical OS prediction nomogram

The scatter plot showed a positive correlation with respect to tumor stage and risk scores (p<0.001) (Supplementary Figure 3C). Age >65 years also tended to be a factor for an increased risk score (p=0.096) (Supplementary Figure 3B). Decision curve analysis (DCA) (Figure 6A) and ROC curves (Figure 6B) (risk AUC=0.765, age AUC=0.646, sex AUC=0.500, and stage AUC=0.815) showed that risk had a superior prediction value compared to most clinical information. A prospective estimator of KIRC patients based on age, risk, sex, and disease staging was conceived (Figure 6C), patients were used to verify its effectiveness, and the results indicated good performance (Supplementary Figure 4). The forecasting value of the nomogram combined with the risk model (Supplementary Figure 4D) was higher than that of the nomogram without the risk model (Supplementary Figure 4E).

Enrichment analysis

KEGG enrichment analysis revealed the related functions of variably expressed genes, including ubiquitin-mediated proteolysis, glycolipid biosynthesis, galactose metabolism, ethoxylate, and dicarboxylate metabolism (Supplementary Figure 5). GSEA demonstrated that different pathways were enriched in different gene sets. Pathway functions enriched in the low-risk category in the GOBP library included mitochondrial gene expression, assembly of the mitochondrial respiratory chain complex, mitochondrial translation, neurotransmitter reuptake, and mitochondrial electron transport from NADH to ubiquinone (Supplementary Figure 6). All the pathway information is shown in Supplementary Table 4.

Tumor mutational burden

The TMB in patients was determined and found to be higher in the high-risk subgroup (Figure 7A). Based on this, we generated a waterfall plot of the top 20 mutant genes according to different subclusters (Figure 7B, 7C). The five most commonly mutated genes in the high-risk patients were VHL (45%), PBRM1 (38%), TTN (17%), BAP1 (16%), and SETD2 (16%). VHL (48%), PBRM1 (42%), TTN (15%), SETD2 (8%) and MUC16 (7%) were prone to mutation in low-risk patients.

Tumor immune infiltration status analysis

The TME analysis suggested that the immune score (Figure 8A) and the ESTIMATE score (Figure 8C) were higher in the high-risk segment (p<0.01), but the stromal score (Figure 8B) did not show dramatic differences. Patients at high risk had significantly lower tumor purity (p<0.01) (Figure 8D). We hypothesized that an immunosuppressive microenvironment was present in the high-risk subgroup that weakened antitumor immunity. Additionally, we showed an association between immune cells and the risk score under different algorithms (Figure 8E). Naive B cells (p<0.05), plasma cells (p<0.05), follicular helper T cells (p<0.01), regulatory T cells (Tregs) (p<0.0001), and M0 macrophages (p<0.01) constituted a greater proportion within the high-risk group, and resting memory CD4 T cells (p<0.05), monocytes (p<0.01), M1 macrophages (p<0.0001), and resting mast cells (p< 0.0001) showed comparatively higher expression in the low-risk group (Figure 8F). The relationships among immune cells and risk scores are shown in Supplementary Figure 7. Accordingly, APC coinhibition (p<0.01) and type II IFN response (p<0.001) functions were inhibited in high-risk KIRC patients, while parainflammation (p<0.01) and T-cell costimulation functions were improved (Figure 8G). It can be speculated that tumor development in high-risk KIRC patients is facilitated by both T-cell parainflammation and costimulation.

Benefits of promotional models in the treatment of KIRC

T-cell exclusion and T-cell dysfunction, together with TIDE scores, were consistently higher in the high-risk patients with ICD-related signatures (Supplementary Figure 8). This may indicate that there is more potential for tumors in the high-risk subgroup to exhibit immune escape, complicating treatment. The antitumor drugs with different mechanisms of action to which tumors in the low-risk subgroup are sensitive are shown in Supplementary Table 5, those to which tumors in the higher-risk subgroup are sensitive are shown in Supplementary Table 6, and those without significant intergroup susceptibility differences are shown in Supplementary Table 7. The high-risk group appears to be more sensitive to drugs that act on the PI3K/mTOR signaling pathway and metabolism pathway, which can guide drug selection.

In vitro experimental validation of risk models

The protein expression levels of selected ICD-related genes in KIRC tissues and normal tissues were visualized with immunohistochemical staining images in the HPA database (Supplementary Figure 9A).

RT-qPCR results showed that among the eight ICD-related lncRNAs, AP000439.3 was highly expressed in the low-risk group, and the remaining seven lncRNAs were highly expressed in the high-risk group (Supplementary Figure 9B).