Immune-related lncRNA pairs as novel signature to predict prognosis and immune landscape in melanoma patients.

ABSTRACT: To investigate immune-related long non-coding RNA (irlncRNA) signatures for predicting survival and the immune landscape in melanoma patients.We retrieved gene expression files from The Cancer Genome Atlas and the Genotype-Tissue Expression database and extracted all the long non-coding RNAs from the original data. Then, we selected immune-related long non-coding RNAs (irlncRNAs) using co-expression networks and screened differentially expressed irlncRNAs (DEirlncRNAs) to form pairs. We also performed univariate analysis and Least absolute shrinkage and selection operator (LASSO) penalized regression analysis to identify prognostic DEirlncRNA pairs, constructed receiver operating characteristic curves, compared the areas under the curves, and calculated the optimal cut-off point to divide patients into high-risk and low-risk groups. Finally, we performed multivariate Cox regression analysis, Kaplan-Meier (K-M) survival analysis, clinical correlation analysis, and investigated correlations with tumor-infiltrating immune cells, chemotherapeutic effectiveness, and immunogene biomarkers.A total of 297 DEirlncRNAs were identified, of which 16 DEirlncRNA pairs were associated with prognosis in melanoma. After grouping patients by the optimal cut-off value, we could better distinguish melanoma patients with different survival outcomes, clinical characteristics, tumor immune status changes, chemotherapeutic drug sensitivity, and specific immunogene biomarkers.The DEirlncRNA pairs showed potential as novel biomarkers to predict the prognosis of melanoma patients. Furthermore, these DEirlncRNA pairs could be used to evaluate treatment efficacy in the future.

Introduction

Malignant melanoma is the most aggressive skin cancer and has a high mortality rate. There were approximately 100,000 newly diagnosed melanoma cases and more than 7000 melanoma deaths in the United States in 2019,[ and the incidence of melanoma is projected to increase globally. According to the GLOBOCAN data, there were 287,723 new melanoma cases and 60,712 deaths worldwide in 2018,[ compared with 232,130 new cases and 55,488 deaths in 2012.[ Melanoma more commonly affects Caucasian white populations than those with pigmented skin such as Hispanics, Africans, and Asians: an estimated 21.2 per 100,000 population were newly diagnosed with melanoma in Western countries between 2006 and 2010, compared with less than 1 per 100,000 in East Asian countries including China, Japan, South Korea, and Singapore.[ More than half of melanoma patients are stage I at the time of diagnosis and have a 5-year survival rate of 60% to 70%; however, in advanced-stage patients, the 5-year survival rate decreases to 19%[ and the median overall survival is 0 to 9 months.[ Owing to the poor prognosis of patients with metastatic malignant melanoma, current efforts tend to focus on improving early diagnosis.

In the past decade, the prognosis of advanced melanoma patients has been greatly improved by the use of immune checkpoint inhibitors (ICIs), including cytotoxic T lymphocyte associate protein-4 (CTLA-4) inhibitor and programmed cell death protein 1 (PD-1)/programmed cell death ligand-1 (PD-L1) inhibitors. However, 40% of patients show limited or no response to immunotherapy.[ Therefore, the development of immune-related biomarkers is of great significance for evaluating the prognosis of patients and their response to treatment with ICIs. Many studies have assessed the relationship between the prognosis of melanoma patients and levels of various non-coding RNAs.[ However, there is a lack of non-coding RNAs that are suitable as biomarkers to evaluate the effects of immunotherapy in melanoma patients. Inspired by the research of Lv et al,[ we aimed to establish a signature combining long non-coding RNA (lncRNA) pairs to assess the prognosis of melanoma patients.

Long non-coding RNAs (LncRNAs) are transcripts that are more than 200 nucleotides in length and have no protein-coding potential. They can regulate chromatin modification, transcription, and post-transcriptional processing at different levels.[ Furthermore, lncRNAs can affect the proliferation and apoptosis of tumor cells, promote their migration and invasion, and mediate the development of drug resistance, as well as regulating the formation of tumor blood vessels.[ Many studies have shown that lncRNAs are related to the activation or suppression of various immune cells and directly or indirectly involved in tumor immune regulation.[ Wu et al identified a signature consisting of 8 immune-related lncRNAs (irlncRNAs) that could predict prognosis and immunotherapy response in bladder cancer.[ Such irlncRNAs are of great value in evaluating tumor characteristics, prognosis, and treatment.[ Compared to multiple gene combinations, two-gene combination strategy can provide relative expression within a sample, leading to more robust and accurate speculation.[ At the same time, the relative expression of gene pairs takes into account all possible combinations.[ Therefore, developing new irlncRNAs to treat melanoma is of great significance and could improve the quality of life of patients. To the best of our knowledge, there have been few attempts to predict prognosis and the immune landscape in cancer using a two-gene combination strategy composed of irlncRNAs. Therefore, we aimed to construct irlncRNA-pair signatures and evaluate their predictive value with respect to diagnostic effect, chemotherapeutic efficacy, and tumor immune invasion.

Materials and methods

Retrieval and collection of transcriptome data and identification of differentially expressed irlncRNAs

A flowchart of the process is shown in Figure 1. RNA sequencing data for 558 normal skin tissue samples (1 from The Cancer Genome Atlas [TCGA] and 557 from the Genotype-Tissue Expression databases) and 471 melanoma tissues (from TCGA) were obtained from the University of California Santa Cruz Xena browser (https://xenabrowser.net/). The fragments per kilobase million values of all samples were normalized to log2 (Fragments Per Kilobase Million + 1). Merged data from different databases in the genetic Ensembl format were annotated using gene transfer format files (http://asia.ensembl.org). Messenger RNAs (mRNAs) and lncRNAs were distinguished according to the corresponding naming rules in the Gene Transfer Format files. Known immune-related genes were downloaded from IMMPORT (https://www.immport.org/), and co-expression analysis was performed between lncRNAs and these genes. The selection criteria for irlncRNAs were: correlation coefficient greater than 0.4 and P value less than .001. IrlncRNAs were confirmed for subsequent analysis. Differentially expressed irlncRNAs (DEirlncRNAs) i melanoma were identified using the “limma” package in the R software. The screening criteria for DEirlncRNAs were: log (fold change) greater than 1.0 and false discovery rate less than 0.05. Heat maps and volcano plots were used to visualize the results.

Figure 1

Flow chart of the study.

Establishment of two-gene combination strategy and risk assessment model

All possible pairs of DEirlncRNAs (denoted lncRNA-A and lncRNA-B, respectively) were identified and assessed based on the following definitions: (I) if the ratio of lncRNA-A expression to lncRNA-B expression was greater than 1, the pair was defined as lncRNA-A|lncRNA-B = 1; (II) if the ratio of lncRNA-A expression to lncRNA-B expression was less than 1, the pair was defined as lncRNA-A|lncRNA-B = 0; (III) for a given pair, if lncRNA-A|lncRNA-B = 1 in more than 80% of all patients, the pair was considered invalid; and (IV) if lncRNA-A|lncRNA-B = 0 in more than 80% of all patients, the pair was also considered invalid.[ Valid DEirlncRNA pairs were identified by applying these definitions and used to construct a matrix with each entry corresponding to the value of the respective lncRNA-A|lncRNA-B. Combining clinical survival status and survival time data from TCGA, prognostic DEirlncRNA pairs were obtained by univariate Cox regression analysis and least absolute shrinkage and selection operator (LASSO) regression analysis.[ Independent prognostic DEirlncRNA pairs were identified by multivariate Cox regression analysis. Then, receiver operating characteristic (ROC) curves were plotted based on the prognostic related genes’ characteristics, and the corresponding area under the curve (AUC) values were determined. Risk scores were calculated for each patient using the following formula:

where coef (DEirlncRNA pairsi) and expr (DEirlncRNA pairsi) represent the survival correlation coefficient of DEirlncRNA pair and the DEirlncRNA pair matrix of patients, respectively.

Using an optimal cut-off point based on the risk score, we reclassified the melanoma cohort into a high-risk group and a low-risk group for further validation.

Clinical evaluation by risk assessment model

Kaplan–Meier (K–M) analysis was used to describe the survival of melanoma patients. We visualized the relationship between melanoma risk and clinicopathological characteristics using a heat map. Univariate and multivariate Cox regression analyses were used to assess patients’ clinicopathological features (including age, sex, clinical stage, T stage, N stage, and M stage) and risk scores to identify prognostic factors. ROC curves were drawn for comparison.

Analyses of tumor immune cell infiltration and expression of ICI-related immunosuppressive molecules and kinase inhibitors

A detailed Spearman correlation analysis was conducted between tumor-infiltrating immune cells and the risk assessment model, and a lollipop chart was constructed.[ To investigate the relationship between risk score and tumor-infiltrating immune cells, 7 methods were used to evaluate immune infiltrating cells: TIMER,[ CIBERSORT,[ xCell,[ quanTIseq,[ MCP-counter,[ EPIC[ and CIBERSORT-ABS.[ We further analyzed the relationships of the risk assessment model with ICI-related biomarkers (CTLA-4 and PD-1/PD-L1) and kinase inhibitor biomarkers (B-Raf Proto-Oncogene, KIT Proto-Oncogene, NRAS Proto-Oncogene, kRAS Proto-Oncogene, and HRAS Proto-Oncogene).

Correlation analysis between risk models and chemotherapy agents

The “Prophetic” R package was used to identify associations between the risk model and chemotherapy (imatinib,[ cisplatin,[ and paclitaxel[) in melanoma using the TCGA melanoma data-set.

Statistical analysis

All analyses were carried out by R version 4.0.3 and corresponding packages. Differences in quantitative data and normally distributed variables were compared using the t-test, and differences in non-normally distributed variables were compared using the Wilcoxon rank-sum test. Differences were compared for more than 2 groups of variables using one-way analysis of variance and the Kruskal–Wallis test. Prognostic analysis was performed using the K–M survival analysis and Cox proportional hazards model.

Results

Identification of DEirlncRNAs

A total of 652 irlncRNAs were identified (more details are shown in Supplemental Digital Content, Table S1). These included 297 DEirlncRNAs, as shown in the heat map in Figure 2A, of which 245 were upregulated and 52 were downregulated (see volcano map in Fig. 2B).

Figure 2

Establishment of DEirlncRNA pairs. (A) Heat map of immune-related lncRNAs with significant differences. (B) Volcano plot of immune-related lncRNAs; green indicates downregulated genes, and red indicates upregulated genes.

A total of 32442 valid DEirlncRNA pairs were obtained, among which 16 independent prognostic DEirlncRNA pairs were identified by Cox regression analysis and LASSO regression analysis. The results of the univariate Cox regression analysis and LASSO regression analysis are shown in Figure 3A and 3B-C, respectively. The results of the multivariate Cox regression analysis are shown in Figure 3D. The 16 independent prognostic DEirlncRNA pairs were as follows: AC087741.2|AC103591.3, AC087741.2|MIR3142HG, AC087741.2|GSEC, ADIRF-AS1|SPINT1-AS1, AC004847.1|LINC02446, AC027130.1|AL034376.1, AC022034.1|FOXP4-AS1, USP30-AS1|HOTAIR, USP30-AS1|AC055854.1, USP30-AS1|SERPINB9P1, USP30-AS1|AL049555.1, USP30-AS1|AP000759.1, AP003392.4|ATP2B1-AS1, AP003392.4|AC009495.2, AP003392.4|RNF144A-AS1, and AP003392.4|LINC01819. The ROC curves for these 16 DEirlncRNA pairs were shown in Figure 4A. The AUC value and the optimal cut-off point were 0.795 and 2.050, respectively. ROC curves for 3-, 5-, and 10-year survival predicted by the risk assessment model are shown in Figure 4B; the 3-, 5-, and 10-year AUC values were 0.744, 0.786, and 0.794, respectively. A comparison of ROC curves with clinical features is shown in Figure 4C.

Figure 3

Establishment of a risk assessment model using DEirlncRNAs. (A) Univariate analysis of DEirlncRNA pairs. (B, C). LASSO was used to screen overall survival related lncRNAs. (D) Multivariate Cox regression analysis of DEirlncRNA pairs.

Figure 4

Establishment of a risk assessment model using DEirlncRNA pairs. (A) ROC curve for 16 DEirlncRNA pairs and optimal cut-off point. (B) 3-, 5-, and 10-year survival ROC curve for risk assessment model. (C) Comparison of ROC curves with other clinical features.

The risk score and survival outcome from TCGA for each case are shown in Figure 5A-B. K–M analysis showed that patients in the high-risk group had a shorter survival time than those in the low-risk group (P < .001; Fig. 5C). The strip chart (Fig. 6A) and scatter plot obtained by Chi-Squared independence test showed that age (Fig. 6B), T stage (Fig. 6C), and N stage (Fig. 6D) were associated with risk score. However, sex, M stage, and clinical stage were not correlated with risk score (P > .05; Supplemental Digital Content, Figure S1).

Figure 5

Clinical evaluation using the risk assessment model. (A) Risk scores and (B) survival outcomes for each case, and (C) K–M analysis.

Figure 6

Clinical evaluation using the risk assessment model. (A) Strip chart and scatter diagram showing that age (B), T stage (C), and N stage (D) were significantly associated with risk score.

Patients in the high-risk group were more positively correlated with the infiltration of certain immune cells, including resting natural killer cells (CIBERSORT-ABS) and natural killer T cells (xCell). By contrast, they were negatively associated with CD4 + memory T cells (xCell) and CD8 + central memory T cells (xCell), etc. Detailed Spearman correlation analysis was conducted; the resulting lollipop graph is shown in Figure 7. More results of the Spearman analysis of tumor-infiltrating immune cells using different software are shown in Supplemental Digital Content, Figures S2–S6. Comparison of the risk model with ICI-related biomarkers showed that the high-risk group was positively correlated with low expression of CTLA-4 (P < .001; Fig. 8A) and PD-1 (P < .001; Fig. 8B). Further comparison of the risk model with kinase inhibitor biomarkers showed that the high-risk group was positively correlated with low expression of B-Raf Proto-Oncogene (P < .001; Fig. 8C), NRAS Proto-Oncogene (P < .05; Fig. 8D), and kRAS Proto-Oncogene (P < .05; Fig. 8E). The high-risk group was positively correlated with high expression of KIT Proto-Oncogene (P < .05; Fig. 8F) and HRAS Proto-Oncogene (P < .01; Fig. 8G).

Figure 7

Estimation of tumor-infiltrating cells. Estimation of tumor-infiltrating cells by the risk assessment model. Patients in the high-risk group were more positively correlated with the infiltration of certain immune cells, including resting natural killer (NK) cells (CIBERSORT-ABS) and NK T cells (xCell). By contrast, they were negatively associated with CD4+ memory T cells (xCell) and CD8+ central memory T cells (xCell), etc.

Figure 8

Estimation of ICI-related immunosuppressive molecules and kinase inhibitor molecules using the risk assessment model. High risk scores were positively correlated with downregulation of CTLA4 (A), PDCD1 (B), BRAF (C), NRAS (D), and KRAS (E), and with upregulation of KIT (F) and HRAS (G). The model is a potential predictor for chemosensitivity as high-risk scores were related to higher IC50 values for cisplatin (H) and lower IC50 values for imatinib (I) and paclitaxel (J).

Estimation of the correlation between the risk model and clinical treatment

Comparing the high-risk group and the low-risk group, we found that in the high-risk group, high half-maximal inhibitory concentration values were observed for certain chemotherapy regimens, including cisplatin (P = 1.5e−07; Fig. 8H); however, there were low half-maximal inhibitory concentration values for other drugs including imatinib (P = 4.4e−06; Fig. 8I) and paclitaxel (P = .037; Fig. 8J).

Discussion

Malignant melanoma is the most aggressive skin cancer and has a high mortality rate.[ The prognosis of advanced melanoma has improved greatly owing to the development of ICI therapy; however, in 60% of patients, ICIs have limited efficacy or elicit no response.[ Therefore, it is important to explore potential therapeutic targets and prognostic indicators in melanoma. For these reasons, research has increasingly focused on lncRNAs, especially irlncRNAs.[ Previous studies have confirmed that lncRNAs have crucial roles in regulating the occurrence and development of various tumor types.[ In the present study, we aimed to construct a novel signature composed of DEirlncRNA pairs to evaluate the prognosis of melanoma patients.

Many studies have assessed the prognosis of melanoma patients by analyzing coding or non-coding RNAs,[ usually based on quantitative analysis of RNA transcription levels. Guo et al analyzed the expression levels of 16 lncRNAs and constructed a signature to predict prognosis of melanoma patients.[ Zhang et al constructed a lncRNA–mRNA regulation module containing 6 lncRNAs and 4 target mRNAs, again for prognostic assessment in melanoma patients.[ However, a problem with such approaches is that the low frequency of a single differentially expressed gene may lead to a high rate of missed diagnoses.[ Inspired the research of Lv et al[ we established a prognostic signature for melanoma based on combined DEirlncRNA pairs. This process does not require measurement of specific transcriptional expression levels of RNAs; instead, a matrix with entries 0 or 1 is constructed based on comparisons of the expression of DEirlncRNAs. In contrast to previous research, our study aimed to identify DEirlncRNAs and establish the most critical DEirlncRNA pairs. Therefore, we only needed to detect high- or low-expression pairs and did not need to examine the specific expression value of each lncRNA. This novel model has excellent clinical practicality; it not only distinguishes high-risk clinical cases but can also reduce the error-correction work required, making it suitable for use by researchers or clinicians from different backgrounds.

With the development of anti-PD-1 and anti-CTLA-4 agents, there has been constant innovation in the treatment of melanoma.[ Our results showed that PD-1 and CTLA-4 expression was downregulated in the high-risk group. It should be noted that PD-1 and PD-L1 blockade is only effective in patients with high PD-1 or PD-L1 expression.[ Topalian et al found that 17 PD-L1-negative patients had no blocking reaction when treated with a PD-1 inhibitor, whereas 9 of 25 PD-L1-positive patients reacted to ICI treatment (P < .05).[ Shang et al found that low expression of PD-1 in uveal melanoma was an indicator of less infiltration of immune cells.[ Therefore, we speculated that the low expression of PD-1 in the high-risk-score group might be accompanied by less infiltration of immune cells combined with PD-1 inhibitor resistance, resulting in a shorter survival time. Higashikawa et al showed that high expression of CTLA-4 was associated with tumor-infiltrating T cells.[ Ciszak et al conjectured that CTLA-4 inhibitors induced survival-promoting signals in tumor cells of patients with high CTLA-4 expression.[ Therefore, the survival time of patients in the low-risk group with high expression of CTLA-4 was longer than that of patients in the high-risk group. Melanoma is often accompanied by an abnormal Mitogen-Activated Protein Kinase pathway; key molecules of this pathway include RAS and RAF.[ Therefore, in addition to ICI-related biomarkers, we compared the expression of key biomarkers of kinase inhibitors in the high- and low-risk groups.

Different types of tumor-infiltrating immune cells can upregulate or downregulate melanoma immunity in the microenvironment.[ By comparing the correlations of immune cells in the high- and low-risk groups, we found that immune cell infiltration was lower in the high-risk group than in the low-risk group. Immune cells are a major component of the tumor microenvironment and are closely related to tumor cell proliferation, treatment response, and prognosis.[ Previous studies have shown that a higher degree of immune cell infiltration in melanoma is associated with a more significant response to ICIs.[ Therefore, based on the combined results of the survival analysis and measurement of immune cell infiltration and expression of ICI markers in the high- and low-risk groups, we speculated that melanoma patients with low levels of immune cell infiltration and ICI marker expression would have low sensitivity to immune-related therapies, resulting in poor prognosis.

When we compared the patients’ clinicopathological characteristics with the results of the risk model, we found that older patients tended to have higher risk scores (P = .049, Fig. 6B). Patients with stage T4 melanoma also had higher risk scores than those in stages T0–T3 (Fig. 6C). Balch et al reached a similar conclusion in a retrospective study; they found that the prognosis of melanoma patients became worse with increasing age, reflecting the more aggressive tumor biology in older patients.[ However, in contrast to Balch et al, we established a prognostic model of DEirlncRNA pairs, starting from irlncRNAs, to more deeply understand the level of immune regulation in patients in different risk groups. We also we used a more convenient way to build DEirlncRNA pairs, enabling us to judge the prognosis of patients without knowing their precise lncRNA expression levels.

However, the study had some limitations. First, the data were obtained from several public databases, and the clinical information downloaded from these databases was limited and incomplete. For example, we could only analyze and the IC50 values of imatinib, paclitaxel, and cisplatin between the risk model; there was a lack of sensitivity analysis for PD-1 and CTLA-4 inhibitors such as nivolumab and ipilimumab.[ Second, with respect to signature preparation, we adopted a more efficient combination strategy. We determined specific transcription levels of DEirlncRNAs in the form of DEirlncRNA pairs; however, the conversion from a pair of expression levels to a 0 or 1 matrix entry is a one-way process. Finally, irlncRNAs are immune-related regulatory factors in the development of melanoma, and their molecular mechanisms need to be further studied to promote their potential clinical applications.

Overall, in this study, we identified a novel signature composed of irlncRNAs that does not require precise measurement of lncRNA expression levels to predict prognosis in melanoma patients. This signature may help to distinguish which patients are more likely to benefit from anti-tumor immunotherapy.

Acknowledgments

We are very grateful for the contributions of TCGA and GTEx database that provide information on cancer research, as well as all colleagues involved in the study.

Author contributions

Conceptualization: Yu Jin.

Data curation: Zhehong Li, Junqiang Wei, Honghong Zheng, Yu Jin.

Formal analysis: Zhehong Li, Junqiang Wei, Yafang Zhang.

Investigation: Honghong Zheng, Mingze Song, Yafang Zhang.

Methodology: Zhehong Li, Honghong Zheng, Xintian Gan.

Resources: Zhehong Li, Mingze Song.

Supervision: Junqiang Wei.

Validation: Zhehong Li, Xintian Gan, Mingze Song, Yafang Zhang, Yu Jin.

Visualization: Zhehong Li, Junqiang Wei.

Writing – original draft: Zhehong Li, Junqiang Wei, Honghong Zheng.

Writing – review & editing: Junqiang Wei, Yu Jin.