An immune-related lncRNA risk coefficient model to predict the outcomes in clear cell renal cell carcinoma

Analysis of differential expression of immune-related lncRNAs in ccRCC

The transcriptome and immune gene-related data of ccRCC were obtained from The Cancer Genome Atlas (TCGA) database and The Immunology Database and Analysis Portal (ImmPort). The Ensembl database was used to annotate and distinguish transcriptome data. Using Pearson correlation analysis, with co-expression correlation coefficient >0.4 and P < 0.001 as the identifying criteria, 433 immune-related lncRNAs were identified. We used differential expression analysis, with |log Fold Change| >1.5 and false discovery rate (FDR) <0.05 as the identifying criteria. We obtained 90 differentially expressed immune lncRNAs, and the gene heatmap (Figure 2A) was generated using R software. Sixteen lncRNAs expressions were downregulated, and 74 lncRNAs were upregulated in ccRCC (Figure 2B and Table 1).

Establishment of differentially expressed immune lncRNA pairs and risk coefficient scoring model

By matching the 90 differentially expressed immune lncRNA pairs for multiple cycles, a total of 2663 differentially expressed immune lncRNA pairs were obtained. Next, 27 immune lncRNA pairs were identified using least absolute shrinkage and selector operation (LASSO) regression analysis and Cox univariate regression analysis (Figure 3A). Then, a Cox multivariate regression analysis was performed based on these 27 immune lncRNA pairs, 16 of which can participate in constructing the risk coefficient scoring model (Figure 3B) and the risk coefficient of each immune lncRNA pair was obtained (Table 2).

Evaluation of the prognostic predictive power of the risk model

Above 27 prognostic-related immune lncRNA pairs were used to construct the 1-year, 3-year, and 5-year receiver operator characteristic (ROC) curves of patients (Figure 4A), and the 1-year area under the curve (AUC) was calculated to be the largest AUC of 0.867 (Figure 4B). In addition, the 3-year and 5-year AUC obtained were 0.832 and 0.838, respectively, which also had predictive power. Through the best fit, the cut-off value for distinguishing between high-and low-risk groups of ccRCC patients was 2.822. We included 190 patients in the low-risk group and 360 patients in the high-risk group.

Analysis of the correlation of clinical indicators by the risk model

The relationship between the risk factor score and the risk subgroup patients ccRCC was analyzed via R language software (Figure 5A). According to the time process, the relationship between the patient’s survival status and the risk coefficient score were obtained (Figure 5B), and the Kaplan-Meier curve was constructed according to the survival status of the high- and low-risk groups (Figure 5C). The result was that the survival rate of patients in the low-risk group was significantly higher than that in the high-risk group (P < 0.001).

The heatmap in Figure 6A described the relationship between the level of risk scores and clinically relevant indicators. We found that the survival status of patients with ccRCC (P < 0.001), tumor grade (P < 0.001), tumor stage (P < 0.001), T stage (P < 0.001), M stage (P < 0.001), and N stage (P < 0.01) were related to risk coefficient score significantly. It can be seen from the box plot we constructed that ccRCC patients with a higher risk factor had a higher chance of death (Figure 6B). Furthermore, tumor grade (Figure 6C), tumor clinical stage (Figure 6D), T stage (Figure 6E), N stage (Figure 6F), and M stage (Figure 6G) were also higher.

A Cox univariate and a multivariate regression analysis were performed on the risk score and clinical correlation indicators. Then the R language software’s survival package was used to visualize the data, and forest maps were done (Figure 7A and Figure 7B). It was found that the tumor grade, tumor stage, TNM stage, and risk coefficient score were related to the outcome of the Cox univariate analysis, but in the Cox multivariate analysis, age, gender, and risk coefficient score were independent predictors of outcome. The ROC curve of clinical-related indicators and the 1-year risk coefficient score were compared in the same chart (Figure 7C). The result was that the patient’s risk coefficient score (AUC = 0.867) and tumor stage (AUC = 0.868) had the highest predictive efficacy.

Correlation analysis between risk coefficient model and immune cell infiltration

XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT were used to estimate the proportion of immune cells in these samples of ccRCC patients based on marker gene and deconvolution algorithm. Pearson correlation test was used to analyze the correlation between the risk coefficient model and tumor immune infiltrating cells with screening criteria P < 0.05, and R language software was used for data visualization (Figure 8). We found that the samples of the high-risk group were positively correlated with the infiltration of NK cells, regulatory T cells, and M1 macrophages in ccRCC and negatively correlated with the infiltration of neutrophils in ccRCC.

Correlation analysis of risk coefficient model with genes

Immune-targeted therapy is one of the most popular drugs for the treatment of renal clear cell carcinoma. We further explored the relationship between the risk coefficient model and genes, and found that among the high-risk patients, the expression levels of CTLA4 (P < 0.001; Figure 9A), LAG3 (P < 0.001; Figure 9B), PDCD1 (P < 0.001; Figure 9C), GAL9 (P < 0.001; Figure 9D), and TIGIT (P < 0.001; Figure 9E) increased. The expression level of PDCD1LG2 increased but not significantly (P > 0.05; Figure 9F). The expression levels of CD274 (P < 0.01; Figure 9G), HAVCR2 (P < 0.01; Figure 9H) decreased in this model. These genes are the potential therapeutic targets for ccRCC.

Correlation analysis between risk coefficient model and targeted therapy drugs

Targeted therapy drugs are the first-line therapy for patients with advanced ccRCC. We also analyzed the relationship between the risk coefficient scoring model and the sensitivity of targeted therapy drugs. The IC₅₀ was used to evaluate the efficacy of drugs. Lower IC₅₀ suggests higher sensitivity. We found that the high-risk group was associated with higher sensitivity of sunitinib (Figure 10A), which was statistically significant (P = 3 e-08); axitinib (Figure 10B), bevacizumab (Figure 10C), pazopanib (Figure 10D), and sorafenib (Figure 10E) were not significantly different in the high- and low-risk groups.

Data acquirement

The GDC Data Transfer Tool was used to download open transcriptome data from TCGA (https://cancergenome.nih.gov/) [31] of ccRCC and normal tissues adjacent to the cancer, which included 539 cases of ccRCC samples and 72 cases of normal tissue samples adjacent to cancer. Then gene transfer format (GTF) files were downloaded using Ensembl (http://asia.ensembl.org) [32] to annotate and distinguish the mRNA and lncRNA of the transcriptome data. Immune-related genes were obtained from ImmPort (http://www.immport.org).

Differential expression analysis of immune-related lncRNAs

Immune-related lncRNAs were identified based on the co-expression strategy and Pearson correlation analysis according to the co-expression correlation coefficient > 0.4 and P < 0.001. DEirlncRNA was selected using the limma package of R language software [33] with |log Fold Change|>1.5 and FDR < 0.05. Obtained lncRNAs were visualized using the heatmap package.

Construction of immune-related lncRNA pairs

DEirlncRNAs were identified using multiple rounds of pairing. Parameter values of 0 or 1 were used for definition, and α was defined as the parameter value. If the expression of lncRNA A was greater than that of lncRNA B in an immune-related lncRNA pair in a sample, then the α value of the lncRNA pair was 1; otherwise, the value of α was 0. If the ratio of the α value (either 0 or 1) of the immune-related lncRNA pair in all samples is less than 80%, it means that the immune-related lncRNA pair was an effective match; otherwise, it needs to be re-paired.

Acquisition of clinical data and establishment of model

First, clinical data related to ccRCC were downloaded from TCGA. Then, the limma package of R software was used to match the immune-related lncRNA pairs in the previous step. We then took the intersection and deleted the repeated clinical data with no follow-up time. A single factor regression analysis was performed on the immune-related lncRNA pairs initially obtained, and the immune-related lncRNA pairs related to the survival status were found. The significance screening criterion was P < 0.01.

To prevent over-fitting, the glmnet package of R language was used to perform LASSO regression analysis [34] on the obtained immune lncRNA pairs, run 1000-repeated random cycles, and immune lncRNA pairs with a matching frequency of more than 100 times identified those with P < 0.05 after the second cross-validation. The best pairing combination was selected to obtain immune lncRNA pairs that can participate in constructing the Cox risk coefficient model. By constructing Cox univariate and multivariate analysis models, the risk coefficient of each immune lncRNA pair related to the outcome was obtained, and the risk score of each patient’s tumor sample was determined. The total risk score of each ccRCC patient sample was equal to the sum of the expression amount of each immune lncRNA pair in the sample multiply risk coefficient. The formula is following:

$$\text{Risk Score}={\displaystyle {\sum }_{i=1}^n{\text{Risk Coefficient}}_i\times \text{I}}{\text{rlncRNA Expression}}_i$$

The Cox analysis results were visualized using the survminer and survival packages of R software.

Construction of the ROC curve with risk coefficient model

ROC curves were constructed using the survivalROC package of the R software, which included ROC for 1-, 3-, and 5- years and the AUC values were calculated to determine the value predicted by the model. We found that the 1-year ROC curve had the largest AUC value. According to the AIC best fit [35], it was possible to distinguish low-risk with high-risk patients by finding a critical value with the largest sum of specificity and sensitivity.

Clinical correlation analysis with risk coefficient model

Survival and survminer packages of the R language software were used to compare the survival differences between the high- and low-risk groups. P < 0.001 indicated a significant difference. A Kaplan-Meier curve was constructed to visualize the data. The relationships between the risk score and the previously obtained clinical indicators (survival status, age, gender, tumor grade, tumor stage, and T, N, and M stages) were analyzed using the chi-square test. The relationships between the risk score and the different subgroups of clinical indicators were analyzed using the Wilcoxon rank-sum test. The limma package and ggpubr package of the R language software were used to visualize the data. To determine whether the risk score can be used as an independent predictor related to the outcome of patients with ccRCC, we performed Cox univariate and multivariate regression analysis on the risk score and clinical correlation indicators, which used the hazard ratio to evaluate. P < 0.05 was the identification criterion, and the survival package of R software was used to visualize the data. To compare the accuracy of the risk score and clinically relevant indicators in predicting survival and outcome, we compared the ROC curves obtained for the 1-year follow-up with the ROC curve of clinically relevant indicators in the same graph.

Correlation analysis of immune cells

To analyze the relationship between risk factor score and immune cell infiltration, the immune cell infiltration data of patients with ccRCC in TCGA was calculated based on CIBERSORT (http://cibersort.stanford.edu/) [36], TIMER (version 2.0; http://timer.cistrome.org/) [37], QUANTISEQ (http://icbi.at/quantiseq) [38], Microenvironment Cell Populations-counter [39], EPIC (http://epic.gfellerlab.org) [40], and XCELL (http://xCell.ucsf.edu/) [41]. The correlation between immune cell infiltration data and risk coefficient score was analyzed using the limma, scales, ggplot2, and ggtext packages of R software, which can be visualized to obtain a bubble chart according to P < 0.05 as the identifying criterion.

Gene correlation analysis

We found that CD274, CTLA4, HAVCR2, LAG3, LGALS9, PDCD1, PDCD1LG2, and TIGIT were abundantly expressed in ccRCC samples. To determine whether these genes differed between the high- and low-risk groups of the risk model, the limma and ggpubr packages of the R language software were used to analyze and visualize the data using violin charts.

Correlation analysis of targeted drugs

To determine whether there was a difference in patients’ response in the high- and low-risk groups of ccRCC patients to targeted drugs, the half-inhibition rate (IC₅₀) of the drug was used as an index to measure drug sensitivity. The data were analyzed and visualized using the limma package, ggpubr, ggplot, and pRRophetic package in R software.