Identification and Validation of Immune-Related lncRNA Signature as a Prognostic Model for Skin Cutaneous Melanoma

Introduction

Malignant melanoma is one of the most aggressive cancers and the fifth most common cancer in the United States.¹ The incidence of skin cutaneous melanoma (SKCM) continues to increase in both males and females.^2,3 In contrast, melanoma cancer mortality has declined with the emergence of targeted therapy and immunotherapy in recent years.⁴ Surgery is the primary treatment option for localized disease with cure rates of 90%.⁵ Other therapeutic approaches such as chemotherapy, radiotherapy, targeted therapy, and immunotherapy are increasingly being used to improve local control. Although the remarkable improvements in targeted therapy and immunotherapy,⁴ the prognosis of metastatic melanoma remains poor due to its intrinsic resistance to chemotherapy or radiotherapy and aggressive clinical behavior.⁶ Therefore, it is crucial to explore new sensitive biomarkers to predict melanoma patients’ prognosis and further guide the proper individual treatment strategies.

Long non-coding RNA (lncRNA) is a class of RNA transcripts more than 200 nucleotides, which take part in gene expression through histone modification, transcriptional, translational, and post-transcriptional regulations.⁷ Although lncRNAs include most transcripts in the mammalian genomes, their biological functions remain largely unknown.⁸ Recent studies indicated that lncRNAs were involved in human diseases such as various cancer.^9–11 For instance, lncRNA DLX6-AS1 promoted proliferation, migration, and invasion of liver cancer via increasing the expression of WEE1.⁹ Jin et al found that overexpression of lncRNA MORT inhibited cell proliferation in oral squamous cell carcinoma through downregulating ROCK1.¹⁰ Various lncRNAs have also been linked with tumor progression in melanoma. MEG3 was found to promote growth, metastasis, and formation of melanoma via regulating the miR-21/E-cadherin axis.¹¹ MALAT1 was found to be highly associated with lymph node metastasis, and UCA1 was related to advanced melanoma.¹² Moreover, lncRNAs regulate genes related to the activation of immune cells, thus leading to tumor immune-cell infiltration while altering the immune microenvironment.^13,14 These findings suggest that lncRNAs could be used as predictive biomarkers for metastatic melanoma.

Therefore, the construction of a model to accurately predict the prognosis may play an important role in improving survival in SKCM patients. The signatures related to tumor immune infiltration reveal promising predictive and prognostic effectiveness in the evaluation, diagnosis, and treatment of cancer.^15–17 The immune-related long non-coding RNAs (IRlncRs) signatures further improve the accuracy of prediction and diagnosis.^18–21 Moreover, the use of dual biomarkers was found to increase the prediction accuracy of diagnostic models for cancers when compared with single genes.²² However, there are few similar models available for SKCM. In this study, a novel modeling algorithm, paring, and iteration were utilized to establish an IRlncRs signature that could be used to improve the diagnosis of SKCM irrespective of the specific expression levels. Then, we validated its predictive value among patients with SKCM, as well as its diagnostic effectiveness, tumor immune infiltration, and chemotherapeutic efficacy.

Materials and Methods

Data Download and Pretreatment

The workflow used throughout this study is illustrated in Figure 1. The RNA-sequencing (RNA-Seq) data of 470 SKCM samples were downloaded from The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/). The corresponding patients’ clinical data, such as age, survival information, and clinical stage, were also derived from TCGA. The clinicopathological characteristics of these patients are listed in Table 1. Gene transfer annotation files were downloaded from Ensembl (http://asia.ensembl.org) to differentiate the messenger RNAs (mRNAs) from lncRNAs for further analysis. The RNA-Seq data files were merged into a matrix file using the merge language script Perl 5.32.1 (http://www.perl.org/). Similarly, we used a symbol script in the Perl language to convert the Ensembl IDs of genes into a matrix of gene symbols.

Table 1 Clinicopathological Characteristics of SKCM Patients from TCGA Database

Figure 1 The process flow of this study.

Furthermore, the immune-related genes (ir-genes) list was downloaded from the ImmPort database (http://www.immport.org). It was used to identify IRlncRs through the co-expression analysis. Ir-genes with correlation coefficients higher than 0.6 and a p-value less than 0.001 were classified as IRlncRs.

Pairing IRlncRs

The IRlncRs were cyclically paired individually, and a 0 or 1 matrix was constructed using the assumption that D is equal to IRlncR B plus IRlncR C; D is defined as 1 if the expression of IRlncR B is higher than IRlncR C, otherwise, D is defined as 0. Then, we further screened the constructed 0 or 1 matrix. If the expression quantity of the IRlncR pairs was 0 or 1, the relationship between pairs and prognosis was disregarded, as the absence of a specific grade could not correctly predict the patient’s survival outcome. The IRlncRs-pairs were considered valid when the expression quantity of the IRlncRs-pairs was either 0 or 1 and accounted for 20% to 80% of the total pairs.

Acquisition of Survival-Related IRlncRs (sIRlncRs)

To get the sIRlncRs, we took the intersection of the survival information and the IRlncRs-pairs using the limma package. The survival time of SKCM patients ranges from 0 to 11,252. Then, univariate Cox regression analysis was performed to identify sIRlncRs from the IRlncRs-pairs. The analysis was carried out through the R packages “survival”, and IRlncRs-pairs with p < 0.01 for survival comparison were considered as sIRlncRs. These sIRlncRs were used for subsequent analysis.

Construction of a Risk Model to Evaluate the Risk Score

The prognosis-related sIRlncRs were identified using the least absolute shrinkage and selection operator (Lasso) Cox analysis with 10,000-rounds cross-validation to prevent overfitting. Then, multivariate Cox regression analysis was performed on these prognosis-related sIRlncRs to establish a prognostic sIRlncRs signature and calculate the coefficients.²³ The forest map was used to visualize the result of multivariate Cox regression analysis. The risk score of the prognostic sIRlncRs signature for each patient was calculated using the following formula:

Risk score = (ExpressionsIRlncRs1×CoefficientsIRlncRs1) + (ExpressionsIRlncRs2×CoefficientsIRlncRs2) + … + (ExpressionsIRlncRs28 × CoefficientsIRlncRs28).

We calculated the prognostic sIRlncRs signature risk score according to a linear combination of the expression level of sIRlncRs weighted by the regression coefficient (β). The β was calculated from the univariate Cox regression hazard analysis using log-transformed hazard ratios (HR).^23–25

Validation of the Constructed Risk Model

We assessed the predictive value of the model using the time-dependent receiver operating characteristic (ROC) curve and calculated the area under the curve (AUC) using the “survivalROC” package.²⁵ The 1-, 2-, and 3-year ROC curves of the model were also plotted. The patients were divided into high-risk and low-risk groups for subsequent analysis according to the optimal cut-off value calculated by the “survminer” package. The Kaplan-Meier survival analysis was used to analyze the survival difference between the two groups. Furthermore, to confirm whether the model was an independent factor for SKCM patients’ survival, we performed univariate and multivariate Cox regression analyses. To confirm the model’s clinical application value, we analyzed the relationship between the model and clinicopathological characteristics with a chi-square test. A band diagram was plotted and labeled according to the p-value as follows: <0.001 = ***<0.01 = * *and <0.05 = *. The Wilcoxon signed-rank test was used to compare the differences in the risk score among groups with different clinicopathological characteristics. The box diagram showed the analysis results. The R packages utilized in these operations were survival, pHeatmap, and ggupbr.

Impact of Tumor-Infiltrating Immune Cells on Risk Scores

To evaluate the relationship between the risk score and immune-cell characteristics, we utilized the widely acknowledged methods to calculate the immune infiltration statuses among the samples from the TCGA project of the SKCM dataset, including TIMER, CIBERSORT, XCELL, QUANTISEQ, MCPcounter, EPIC, and CIBERSORT-ABS. The Wilcoxon signed-rank test was performed to analyze the differences in tumor immune infiltrating cell content between high-risk and low-risk groups of the model. Spearman correlation analysis was used to evaluate the relationship between the risk score and the immune infiltrated cells. The results of both tests were further illustrated through the use of a boxplot diagram. The data were analyzed using the R ggplot2 packages, and a p-value below 0.05 was considered statistically significant.

Clinical Performance of the Model for Identification of Gene Mutations and Treatment

According to the national comprehensive cancer network (NCCN), paclitaxel, vinblastine, and cisplatin are commonly used cytotoxic drugs for the management of SKCM. The clinical performance of the model was estimated by calculating the half-maximal inhibitory concentration (IC₅₀) of these drugs using the SKCM dataset of the TCGA. Wilcoxon signed-rank test analyzed the difference in the IC50 between the high-risk and low-risk groups, using R language loaded with packages pRRophetic and ggplot2. To further explore the clinical performance of the model, we investigated the relationship between this model and immune checkpoint inhibitors including, CTLA4 and PD1. Furthermore, according to the NCCN guidelines, we evaluated the correlation between the risk and specific mutant genes (BRAF, KIT, and NRAS) using R packages limma and ggupbr.

Results

Establishment of IRlncR Pairs and a Risk Assessment Model

A total of 753 IRlncRs were identified following co-expression analysis between known ir-genes and lncRNAs. These are summarized in Table S1. The iteration loop and the 0 or 1 matrix screening identified 9107 valid IRlncR pairs, following the removal of clinical data with non-complete survival information. The IRlncR pairs were added with survival information by intersecting complete tumor survival information and IRlncR pairs. After univariate Cox regression analysis, a total of 2345 sIRlncRs pairs were identified, as illustrated in Table S2. Using modified Lasso regression analysis (Figure 2A and B, and Table S3), 28 sIRlncRs were extracted, 17 of which were included in a Cox proportional hazards model using the stepwise method (Figure 2C). The AUC of the one-year ROC of the 17 pairs resulted in an optimal cut-off point of 1.613 (Figure 3A). The two- and three-year ROCs further validated the optimal cut-off point with an AUC above 0.784 for both curves (Figure 3B). ROC curves were also used to validate the impact of clinical characteristics on the risk score (Figure 3C). A total of 454 cases of patients with SKCM were identified from TCGA. The risk scores were calculated for all cases.

Figure 2 Identification of survival-related IRlncRs using LASSO regression analysis and multivariate Cox regression analysis. (A) LASSO coefficient profiles of the 28 sIRlncRs of SKCM. (B) Plots of the cross-validation error rates. Each dot represents a lambda value along with error bars that represent the confidence interval for the cross-validated error rate. The top of the plot gives the size of each model. The vertical dotted line indicates the value with the minimum error and the largest lambda value where the deviance is within one SE of the minimum. (C) Forest plots of HR of sIRlncRs pairs obtained by multivariate Cox regression analysis. A total of 6 DEGs were found to be prognostic factors. The genes with HR < 1 are protective factors, while the ones with HR > 1 are risk factors in SKCM. Abbreviations: IRlncRs, immune-related long non-coding RNAs; LASSO, least absolute shrinkage and selection operator; SKCM, skin cutaneous melanoma; HR, hazard radio; DEGs, differentially expressed genes.

Figure 3 Primary evaluation of the risk model by sIRlncRs pairs. (A) Plot a curve of AUC value generated by ROCs of 1-year and to identify the optimal cut-off value of the AUC. (B) The 1-, 2-, and 3-year survival receiver operating characteristic curves. (C) A comparison of 1-year ROC curves with other common clinical characteristics showed the superiority of the risk score. Abbreviations: sIRlncRs, survival-related IRlncRs; AUC, area under the curve; ROC, receiver operating characteristic.

Clinical Evaluation by Risk Assessment Model

According to the optimal cut-off point validated previously, we divided the SKCM samples into the high-risk group and the low-risk group. The risk curve and scatterplot were drawn to the risk score and survival status of each SKCM sample (Figure 4A). These results suggested that patients in the low-risk group had better clinical outcomes than those in the high-risk group. Kaplan-Meier curve showed that patients in the low-risk group had better overall survival (OS) than those in the high-risk group (p<0.001) (Figure 4B). The result of univariate analysis exhibited that age (P < 0.001), clinical stage (P < 0.001), T stage (P < 0.001), N stage (P < 0.001) and risk score (P < 0.001) were significantly correlated with OS (Figure 5A). Multivariate analysis further suggested that risk score (P < 0.001) was an independent prognostic predictor for SKCM patients (Figure 5B). Furthermore, chi-square tests found the relationship between the prognosis of SKCM and clinicopathological characteristics, including age, gender, clinical stage, T stage, N stage, M stage. The strip chart (Figure 6A) and scatter diagrams showed that T stage (Figure 6B), N stage (Figure 6C), M stage (Figure 6D), and age (Figure 6E) were significantly related to the risk score.

Figure 4 Risk model for prognosis prediction. (A) From top to bottom are the risk score, survival status distribution of each patient. (B) Overall survival curves of the prognostic signature, in which the blue line represents the low-risk subgroup and the red line represents the high-risk subgroup.

Figure 5 Univariate and multivariate Cox regression analysis of the risk model and clinicopathologic parameters (A-B).

Figure 6 Clinical evaluation by the risk model. A strip chart (A) along with the scatter diagram showed that T stage (B), N stage (C), M stage (D), and age (E) were associated with the risk score. Notes: *P < 0.05 and **P < 0.01.

Association of Tumor-Infiltrating Immune Cells with the Risk Assessment Model

Because the lncRNAs and infiltrating immune cells played essential roles in the tumor microenvironment, we investigated the relationship between the model and the tumor immune microenvironment. We found that the immune score and microenvironment score was significantly higher in the low-risk group than in the high-risk group by the Wilcoxon signed-rank test (Figure 7A and B p<0.001). Furthermore, we showed the differences in infiltrating immune cells between high-risk and low-risk groups in Figure S1–7. A detailed Spearman correlation analysis was conducted, and the bubble chart exhibited that the high-risk group was more negatively associated with tumor-infiltrating immune cells such as B cell, CD4⁺ T cells, monocytes, and CD8⁺ T cells, whereas they were positively associated with M0 macrophages, and resting NK cell (Figure 7C, and Table S4).

Figure 7 Estimation of tumor-infiltrating cells by the risk model. Patients in the high-risk group were more negatively associated with tumor-infiltrating immune cells (A) and microenvironment (B). The bubble chart (C) exhibited the detailed correlation between different tumor-infiltrating immune cells.

Verification of the Correlation Between the Risk Model and the IC₅₀ of Common Chemotherapeutic Drugs

Because cisplatin, vinblastine, and paclitaxel are recommended for cutaneous melanoma treatment by NCCN guidelines, we assessed the correlation between the model and the efficacy of common chemotherapeutics in treating cutaneous melanoma in the TCGA project of the SKCM dataset. The high-risk score was significantly associated with a higher IC₅₀ of chemotherapeutics such as cisplatin (Figure 8A, p<0.001), vinblastine (Figure 8B, p<0.001), and paclitaxel (Figure 8C, p=0.097). Furthermore, we evaluated whether the model could predict patients’ response to immune checkpoint inhibitors including, CTLA4 and PD1. The results showed that patients with low risk-score have better therapeutic effects than those with high risk-score (Figure 9A–D, p<0.01). These data suggested that the model might predict the treatment response to specific chemotherapy agents and immunotherapies.

Figure 8 Assessment of the IC50 of common chemotherapeutic drugs by the risk model. Box plots revealed that the high-risk score was significantly associated with a higher IC50 of chemotherapeutics such as cisplatin (A, p<0.001), vinblastine (B, p<0.001), and paclitaxel (C, p=0.097). Abbreviation: IC50, half-maximal inhibitory concentration.

Figure 9 The association between the risk model and the efficacy of immunotherapies (A–D).

Correlation Between the Risk Model and Specific Mutant Genes in SKCM

Specific gene mutations (BRAF, NRAS, KIT), which are the most common mutant genes in melanoma, had important diagnostic and prognostic implications according to the NCCN guideline. We also assessed the correlation between the risk and specific mutant genes in the TCGA project of the SKCM dataset. The high-risk group had a positive correlation with the presence of particular mutant genes, including BRAF (Figure 10A, p<0.001), KIT (Figure 10B, p<0.001), and NRAS (Figure 10C, ns). The model confirmed that the presence of these mutant genes has important diagnostic and prognostic implications in SKCM.

Figure 10 Verification of the correlation between the risk model and the most frequently mutated genes in SKCM. Violin plots showed that the high-risk group had a positive correlation with the presence of particular mutant genes, including BRAF ((A), p<0.001), KIT ((B), p<0.001), and NRAS ((C), ns). Abbreviation: SKCM, skin cutaneous melanoma.

Discussion

Despite the recent development of immunotherapy and targeted therapy in the management of SKCM, the treatment outcomes for patients with advanced disease are still poor.²⁶ Tumor molecular features impact treatment response and survival in patients with SKCM, even if these patients have similar clinical risk factors.²⁷ Thus, apart from traditional clinical risk factors, it is imperative to identify additional molecular prognostic markers. Previous studies have concluded that the immune system plays both positive and negative roles in regulating tumorigenesis, cancer progression, and metastasis.²⁶ Correspondingly, lncRNAs are involved in regulating the gene expression of the immune system.¹⁴ Therefore, we constructed a model focusing on IRlncRs and tumor-infiltrating immune cells to guide the diagnosis, treatment, and prognosis of SKCM. This study established an excellent model with two-lncRNA combinations, which was inspired by the strategy of ir-gene pairing. The advantage of this technique is that it is not dependent on the expression levels of the gene signature.²⁸

In this study, we downloaded raw datasets of lncRNAs from TCGA and adopted a differential co-expression analysis to pair IRlncRs. Then, IRlncRs-pairs were identified using an improved method of cyclically single pairing along with a 0 or 1 matrix. We further screened for sIRlncRs-pairs using univariate, Lasso, and multivariate Cox regression analysis. The AUC value of the 1-year ROC curve was calculated to estimate the optimal cut-off point to differentiate between the high-risk and low-risk patients.In addition, we assessed the impact of various clinical variables, including clinical-pathological characteristics, tumor-infiltrating immune cells, chemotherapy, and gene mutation, on the prediction model’s accuracy.

A 6-lncRNA signature was set up for predicting prognosis in melanoma patients.²⁹ Generally, lncRNAs possess significant biological functions and characteristics. Our algorithm could initially identify IRlncRs and pair the most significant IRlncRs. Thus, the higher or lower expressed pairs could be detected instead of calculating each lncRNA’s specific expression value. This novel model has a clinical advantage as it facilitates the distinction between high- and low-risk groups. As lncRNAs are correlated with ir-genes, these lncRNAs may be involved in regulating the tumor immune microenvironment or the immune cells. The models of IRlncRs were applied to predict the prognosis of breast cancer,³⁰ hepatocellular carcinoma,³¹ clear cell renal cell carcinoma,²⁰ and pancreatic cancer,³² except for melanoma. In this study, we first established a novel model with lncRNAs signature to predict melanoma prognosis. In this model, several sIRlncRs had been found to play an important role in the malignant phenotypes of various cancer types, such as AC022509.2,³⁰ LINC02544,³³ AC091057.1,³⁴ HLA-DQB1-AS1,³⁴ SEMA6A-AS1,³⁵ AL365361.1,^36,37 and UBA6-AS1.³⁸ For example, UBA6-AS1 promoted GBM’s malignant progression by targeting the miR-760/HOXA2 axis. Therefore, this model could be used to identify novel biomarkers.

To increase the efficacy and accuracy of the risk prediction model, we calculated every AUC value to determine the maximum value for the optimal model and then compared it with other clinical parameters. We also calculated an optimal cut-off point with the highest specificity and sensitivity to distinguish between the high- and low-risk groups instead of just using the median value.Then, we reevaluated the survival outcome by performing the chi-square test and Wilcoxon signed-rank test to identify the impact of clinicopathological characteristics, chemotherapy efficacy treatment, tumor immune infiltration, and mutant gene. These findings suggested that the upgraded algorithm worked well irrespective of the patients’ clinicopathological characteristics.

To investigate the relationship between tumor-infiltrating immune cells and risk scores, we used seven standard acknowledged methods to estimate the immune-infiltrating cell, including XCELL,³⁹ TIMER,⁴⁰ QUANTISEQ,^41,42 MCPcounter,⁴³ EPIC,⁴⁴ CIBERSORTABS,⁴⁵ and CIBERSORT.^46,47 Through the above-integrated bioinformatics analysis, we found that sIRlncRNA pairs were more positively correlated with tumor-infiltrating immune cells such as CD4+T cells, natural killer cells, neutrophils, and myeloid dendritic cells. These immune cells were an important part of the tumor microenvironment and involved in tumorigenesis, invasion, and metastasis.⁴⁸ Our model suggested that the risk scores were associated with sensitivity to chemotherapeutics (vinblastine, cisplatin), immunotherapies and gene mutations (BRAF, KIT), which further validated the model’s prediction accuracy.

However, this study has some shortcomings and limitations. For example, the raw data downloaded from TCGA was relatively insufficient. We did not manage to retrieve datasets simultaneously, including information on lncRNA expression levels, survival outcomes, and clinicopathological characteristics in patients with SKCM. The constructed model required external validation as each sample had a different expression level, which may lead to an unreliable model. To overcome this problem, a 0 or 1 matrix was created to screen all IRlncRs-pairs and hence minimize sample errors due to differences in expression. Moreover, various methods were used to validate and optimize this novel modeling algorithm. Our model was acceptable based on these results despite the lack of external data and experimental validation. However, external validation by other clinical datasets and experiments is recommended.

Conclusion

A novel prognosis model for SKCM was constructed based on 17 pairs of IRlncRs. These IRlncRs might participant in the development and prognosis of SKCM and could be a powerful indicator of the clinical outcome of SKCM patients. However, further research is required to verify the efficacy of this model clinically.