Identification and Construction of a Predictive Immune-Related lncRNA Signature Model for Melanoma

Introduction

Cutaneous melanoma is a tumor characterized by the abnormal proliferation of epidermal melanocytes that are genetically altered by the interaction of genetic, physical, and environmental factors.^1,2 Melanoma is the most invasive form of skin cancer. Risk factors for melanoma can be divided into internal and external factors and include light skin color, UV exposure, nevi, genetic predisposition, and a family history of melanoma.^3–6 In 2019, there were an estimated 95,830 cases of skin melanoma and 7000 deaths from skin melanoma in the United States.⁷ Depending on the features of the tumor (location, stage, and genetic profile), the therapeutic options may include surgical resection, chemotherapy, radiotherapy, photodynamic therapy (PDT), immunotherapy, or targeted therapy.⁸ The development of targeted therapies (such as BRAF and MEK inhibitors) and immunotherapies (such as anti-PD-1 antibodies alone or in combination with anti-CTLA-4 antibodies) has revolutionized the systemic treatment of advanced melanoma, and these new drugs have become pretreatments recommended by many international guidelines for melanoma treatment.⁹ The crosstalk between tumor and immune cells significantly affects tumor invasion, clinical response and treatment outcome.¹⁰ The immune system can detect and destroy tumor cells through immunosurveillance or play a tumor-promoting role via suppressed immune activity and enhanced anti-inflammatory responses; the immune environment varies greatly from one individual to another, and thus, the prognosis and response to treatment also vary among individuals.¹¹ Currently, no prognostic markers for melanoma have been widely recognized, and it is of great significance to study the immune molecular mechanism underlying melanoma and to discover new immune checkpoints.

Many molecules, especially long noncoding RNAs (lncRNAs), play an important role in the progression of melanoma. For example, studies showed that HOTAIR was overexpressed in metastatic tissues, which indicated that HOTAIR could promote the ability of melanoma cells to migrate and invade and that lncRNAs may be involved in the metastasis of melanoma.^12,13 In in vitro studies, knocking out the UCA1 or MALAT-1 lncRNA could reduce the migration of melanoma cells, which indicates that the increased expression of the UCA1 and MALAT-1 lncRNAs might be related to melanoma metastasis.¹⁴ Li et al found that abnormally upregulated expression of a new lncRNA, BANCR, was involved in the proliferation of melanoma cells in vitro and in vivo.¹⁵ However, there are still few studies on the application of immune-related lncRNAs in melanoma.

This study identified immune-related lncRNAs in melanoma and established a model related to melanoma survival based on these immune-related lncRNAs. This study provides a reference for future molecular diagnosis and treatment.

Materials and Methods

Acquisition and Analysis of Melanoma Expression Data

We downloaded the RNA sequencing (RNA-seq) expression and clinical data of melanoma specimens from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/). Moreover, immune-related genes were retrieved from immunological gene sets on the official gene set enrichment analysis (GSEA) website. Immune-related lncRNAs with a P value less than 0.05 and a coefficient greater than 0.4 were identified through correlation analysis between immune-related genes and lncRNAs.

Cox Regression Analysis of the Immune-Associated lncRNA Signature

The combination of clinical survival data and immune-related lncRNA expression data was used to identify prognostic immune-associated lncRNAs by univariate Cox regression analysis. Multivariate regression analysis was also performed to calculate the risk score. The formula of the risk score was as follows: Risk score = coefficient1 * Expression_gene1 + … + coefficientN * Expression_geneN. Melanoma specimens were categorized into high- and low-risk groups based on the median risk score.

Independent Prognostic Analysis in Melanoma Patients

Clinical characteristics, including age, sex, stage, T stage, M stage and N stage, and the risk score were chosen as independent prognostic risk variables according to univariate and multivariate Cox regression analyses. In addition, principal component analysis (PCA) is a dimensionality reduction method often used in image processing. PCA was utilized to display all genes, immune-related lncRNAs, and high- and low-risk score-associated genes.

Statistical Analysis

R software 3.6.0 (ATandT Labs Research - Software Tools, RRID:SCR_002937, https://www.r-project.org/) was used for all the statistical analyses involving differential expression, Kaplan-Meier curves, Cox regression analysis, forest plots and PCA. A difference of P<0.05 was considered statistically significant.

Results

Differentially Expressed lncRNAs in Melanoma Patients

The clinical materials and lncRNA expression data of both normal specimens and melanoma specimens were collected from TCGA database. Additionally, the relationship between immune gene expression and lncRNA expression in melanoma patients was determined, and 6359 immune-associated lncRNAs with P < 0.05 and coefficient > 0.4 were extracted. There were 6351 positively regulated immune-related lncRNAs and 8 negatively regulated immune-related lncRNAs (see Supplementary Table 1). Interestingly, the former was much higher than the latter. There were a total of 180 females and 290 males included in the study, with a mean age at diagnosis of 58.22 years. A total of 248 patients survived, and 222 patients died.

Construction of a Model for Melanoma Patients

From the univariate Cox regression analysis, 48 prognostic immune-associated lncRNAs with P < 0.01 were identified (see Supplementary Table 2). Furthermore, 14 independent prognostic immune-related lncRNAs (Table 1) were chosen to construct a survival model based on the multivariate Cox regression analysis: risk score = 0.554 * AL662844.4 + 0.257 * AC009065.3 + (−0.295) * LRRC8C-DT + 0.308 * AC099811.3 + (−0.243) * AC021188.1 + (−0.767) * MALINC1 + 0.254 * AC125807.2 + 0.386 * SPINT1-AS1 + (−0.437) * CCR5AS + (−0.568) * EIF2AK3-DT + (−0.342) * AC022306.2 + 0.349 * AC009495.2 + (−0.334) * AC242842.1 + (−0.081) * AL034376.1. As suggested by the results, LRRC8C-DT, AC021188.1, MALINC1, CCR5AS, EIF2AK3-DT, AC022306.2, AC242842.1 and AL034376.1 had a positive impact on the overall survival of melanoma patients, whereas AL662844.4, AC009065.3, AC099811.3, AC125807.2, SPINT1-AS1 and AC009495.2 had a negative impact on the overall survival of melanoma patients. On the basis of the median risk score value, the melanoma samples were divided into high- and low-risk groups. The high-risk group had a lower overall survival than the low-risk group (Figure 1). Figures show the distribution of the risk score, overall survival status and heatmap of the expression of the fourteen immune-related lncRNAs in melanoma patients (Figure 2).

Table 1 Independent Prognostic Risk Immune-Related lncRNAs

Figure 1 Kaplan-Meier overall survival curves of the risk stratification groups. The tick marks on the curve represent the censored subjects. The number of patients at risk is listed below the curve.

Figure 2 Construction of immune‐related lncRNAs correlated with melanoma prognosis. The distribution of the risk score (A), overall survival status (B) and heatmap (C) of the fourteen immune-related lncRNAs in melanoma patients. The black dotted lines represent the median risk score cutoff value dividing patients into low- and high-risk groups. The red dots and lines represent patients in the high-risk group. The green dots and lines represent patients in the low-risk group.

Prognostic Signature of the Risk Score Combined with Clinical Variables

To assess whether the risk score predicts the prognosis of melanoma independently from other clinicopathological features, we performed a Cox regression analysis. The results of the univariate analyses showed a statistically significant relationship between age, stage, T-stage, N-stage, risk score and survival outcome (P < 0.001) (Figure 3A and Table 2). Multifactorial analyses showed that age, T stage, N stage, and risk score were independent prognostic factors for melanoma patients (P < 0.001) (Figure 3B and Table 2). To assess the sensitivity and specificity of the risk score in predicting the prognosis of melanoma patients, we performed receiver operating characteristic analysis. The area under the ROC curve of risk scores was 0.786, which was larger than the areas under the ROC curves of other clinical variables (Figure 4). Finally, the PCA results are presented in Figure 5. There were no significant differences between the high-risk and low-risk groups in terms of the expression of all genes (Figure 5A) or immune-related lncRNAs (Figure 5B). There was a significant difference between the two risk groups in the expression of the fourteen immune-related lncRNAs (Figure 5C) used in the prognostic model. In conclusion, these results indicate that the immune-related lncRNA signature identified above is an independent prognostic factor of melanoma.

Table 2 Univariate and Multivariate Cox Proportional Hazards Regression Analysis of the 14-lncRNA Signature and Clinical Risk Factors in the Entire Dataset

Figure 3 Univariate and multivariate Cox regression analyses of clinical factors associated with overall survival. (A) Univariate Cox regression analyses. (B) Multivariate Cox regression analyses.

Figure 4 Receiver operating characteristic (ROC) curve. The ROC curve of the risk score had the largest AUC of 0.786 compared with the other clinical variables.

Figure 5 Principal component analysis. PCA of all genes (A), immune-related lncRNAs (B), and the two risk groups (C).

Discussion

Cutaneous melanoma is the most common malignant skin carcinoma and has a high incidence. For patients with early-stage (stage I–IIIB) melanoma, surgery is the primary treatment, while immunotherapy has played an active role in treating unresectable or metastatic tumors.⁸ However, similar to the problems encountered in many tumor treatments, immunotherapy has various immunosuppressive mechanisms, the immune environment varies greatly from one individual to another, the prognosis and response to treatment vary among individuals, and only a small percentage of patients benefit from immunotherapy.^11,16 Hence, in-depth exploration of the immunomolecular treatment of melanoma will provide new treatment and diagnosis options. For example, SNHG5 promoted melanoma development by inhibiting miR-26a-5p expression and inducing TRPC3 expression, suggesting the potential of SNHG5 as a novel target therapy for melanoma.¹⁷ In this study, we identified fourteen immune-related lncRNAs and used these lncRNAs to build a prognostic model. The AUC of the risk score was higher than that of other clinical variables.

Recent sequencing discoveries have demonstrated that there are various mechanisms underlying lncRNA expression in cells, ranging from gene expression to protein translation and stability.¹⁸ LncRNAs have been proven to participate in different processes of melanoma, such as progression and metastasis. Fourteen immune-related lncRNAs, LRRC8C-DT, AC021188.1, MALINC1, CCR5AS, EIF2AK3-DT, AC022306.2, AC242842.1, AL034376.1, AL662844.4, AC009065.3, AC099811.3, AC125807.2, SPINT1-AS1 and AC009495.2, were shown to have a great influence on melanoma. To investigate whether the immune-associated lncRNA signature can be an independent prognostic factor for cutaneous melanoma, we performed univariate and multivariate Cox analyses, which showed that the lncRNA signature is an independent prognostic factor for cutaneous melanoma independent of clinicopathological features. ROC curves and PCA proved the accuracy of the prognostic model.

For melanoma, which is a tumor with a poor prognosis, the establishment of an effective prognostic model and the exploration of its mechanism of action can improve the outcome of treatment. The expression levels of LRRC8C−DT were positively correlated with the abundances of M0 and M2 macrophages but negatively correlated with the abundance of CD8⁺ T cells.¹⁹ M1 macrophages exhibit proinflammatory activity and play an antitumor role, while M2 macrophages exhibit anti-inflammatory activity and promote tumor growth. It has been shown that the antitumor effects of immune checkpoint inhibitors can be enhanced by repolarizing M2 macrophages into M1 macrophages, and CD8 T cells located at the melanoma tumor invasive margin can predict the clinical response to PD-1 blockade therapy and are positively correlated with the response to pembrolizumab.^20,21 AC021188.1 has been found to be associated with prognosis and was included in 5-disease prognostic signature lncRNAs in head and neck squamous cell carcinoma (HNSCC).^22,23 Bida et al revealed the effect of the lncRNA MA-linc1 on cell cycle progression and cancer growth and found that its high expression was correlated with lower survival in breast cancer and lung cancer patients, and this might be related to the inhibition of MA-linc1-enhanced anti-mitotic drug paclitaxel-induced apoptotic cell death.²⁴ CCR5 expression sustains myeloid-derived suppressor cell (MDSC) suppression activities, intratumoral Treg infiltration, and melanoma tumor growth and is highly attractive as a means to quench or eliminate unconstrained tumor cell growth.^25,26 A study also showed that the expression of CCR5AS could affect HIV infection and disease progression. For example, inhibiting CCR5AS expression decreased the infection of CD4+ T cells in vitro.^27,28 AC242842.1 is a pyroptosis-related and immune-related lncRNA signature that plays a key role in melanoma and can also be used for predicting prognosis.^5,29 AL662844.4 is an autophagy-related lncRNA that plays a key role in bladder cancer.³⁰ AC125807.2 is a potential prognostic biomarker in lung adenocarcinoma.^31,32 LncRNAs such as AC009495.2 can distinguish acute myeloid leukemia types and significantly change the behavior of acute myeloid leukemia cells.^33,34 Additionally, SPINT1-AS1 plays a role in esophageal squamous cell carcinoma,³⁵ breast cancer,³⁶ renal clear cell carcinoma and other tumors.³⁷

The advantages of this study include that it was a TCGA-based study that included a large sample size and adjusted for the patient’s clinical and demographic characteristics. However, there was a limitation in this study: we did not provide information on experimental data for deep verification, and we will improve this aspect of the study in the future.

Conclusion

This study identified several clinically significant immune-related lncRNAs based on data from cutaneous melanoma patients in TCGA database and established a relevant prognostic model. This signature provides new potential biomarkers for prognosis and immunotherapy selection in cutaneous melanoma.