The new horizon of biomarker in melanoma patients: A study based on autophagy-related long non-coding RNA.

ABSTRACT: Autophagy-related long non-coding RNAs (arlncRNAs) play a crucial role in the pathogenesis and development of the tumor. However, there is a lack of systematic analysis of arlncRNAs in melanoma patients.Melanoma data for analysis were obtained from The Cancer Genome Atlas (TCGA) database. By establishing a co-expression network of autophagy-related mRNAs-lncRNAs, we identified arlncRNAs in melanoma patients. We evaluated the prognostic value of arlncRNAs by univariate and multivariate Cox analysis and constructed an arlncRNAs risk model. Patients were divided into high- and low-risk groups based on the arlncRNAs risk score. This model was evaluated by Kaplan-Meier (K-M) analysis, univariate-multivariate Cox regression analysis, and receiver operating characteristic (ROC) curve analysis. Characteristics of autophagy genes and co-expressive tendency were analyzed by principal component analysis and Gene Set Enrichment Analysis (GSEA) functional annotation.Nine arlncRNAs (USP30-AS1, LINC00665, PCED1B-AS1, LINC00324, LINC01871, ZEB1-AS1, LINC01527, AC018553.1, and HLA-DQB1-AS1) were identified to be related to the prognosis of melanoma patients. Otherwise, the 9 arlncRNAs constituted an arlncRNAs prognostic risk model. K-M analysis and ROC curve analysis showed that the arlncRNAs risk model has good discrimination. Univariate and multivariate Cox regression analysis showed that arlncRNAs risk model was an independent prognostic factor in melanoma patients. Principal component analysis and GSEA functional annotation showed different autophagy and carcinogenic status in the high- and low-risk groups.This novel arlncRNAs risk model plays an essential role in predicting of the prognosis of melanoma patients. The model reveals new prognosis-related biomarkers for autophagy, promotes precision medicine, and provides a lurking target for melanoma's autophagy-related treatment.

Introduction

Malignant melanoma is a highly aggressive malignant tumor of the skin that originates from melanocytes, the GLOBOCAN 2018 database demonstrated that nearly 287,723 new melanoma cases were diagnosed and 60,712 deaths worldwide,[ the incidence rate of melanoma increased by 2.7% annually.[ Surgery remains the preferred treatment for early stage melanoma patients.[ However, many malignant melanoma patients suffer from distant metastasis after surgery, and Haydu LE research showed that 37.1% of patients with stage III developed distant metastasis.[ Therefore postoperative adjuvant therapy is a fundamental step in improving the melanoma patients prognosis.[ Early-stage localized melanoma is curable, but metastatic melanoma remains a huge challenge to treat. The challenge for malignant melanoma treatment is the low response rate of conventional regimens and the limitations of targeted therapy.[ Therefore, it is necessary to develop novel prognostic biomarker prediction models and predict the prognosis of melanoma patients.

Autophagy is a process whereby cytoplasmic components are encapsulated in double-membrane vesicles and transported to lysosomes for degradation.[ Autophagy disorders and a variety of diseases often go hand in hand, including neurodegenerative diseases,[ cardiomyopathy,[ type II diabetes,[ and cancer.[ More and more studies have shown that autophagy plays a pivotal part in developing melanoma.[ In the complex network of intracellular gene regulation, regulating the expression of autophagy-related proteins is a fundamental way to regulate cellular autophagy, and long non-coding RNAs (lncRNAs) are a crucial part of regulating autophagy-related proteins.[ LncRNA is defined as a transcript of more than 200 nucleotides that is not translated into protein.[ Current studies have demonstrated that lncRNAs have many functions, including transcriptional regulation in cis or trans, organization of nuclear domains, and regulation of protein or RNA.[ Therefore, an in-depth study of the role of autophagy-related long non-coding RNAs (arlncRNAs) in autophagy-related will help realize arlncRNAs as a new target for individualized tumor diagnosis and treatment.[ Many studies believe that precise regulation of cellular autophagy can be a new direction for tumor therapy.[ However, the role of arlncRNAs in melanoma patients remains unclear. Therefore, it is of high value to identify key prognostic arlncRNAs in melanoma patients.

This study screened out prognostic arlncRNAs by analyzing transcriptome data of The Cancer Genome Atlas (TCGA) database. Meanwhile, we establish an arlncRNAs-based prognostic risk model to explore the diagnostic and prognostic value of arlncRNAs in melanoma patients.

Materials and methods

Flowchart of the study process

The design idea and workflow are shown in Figure 1.

Figure 1

The design idea and work flow.

Data acquisition and collation

Melanoma sequencing data sets were downloaded from TCGA database (https://portal.gdc.cancer.gov/). We annotated the data by gene transfer format files from Ensembl (http://asia.ensembl.org). Meanwhile, the corresponding clinical data were also obtained from the TCGA database, and patients with complete follow-up data were selected for further analysis. Autophagy-related genes (ARGs) were obtained from the human autophagy database (http://autophagy.lu/clustering/index.html).

Identification of arlncRNAs

Transcriptome profiling data were divided into message RNAs (mRNAs) and lncRNAs by gene transfer format files. According to Pearson correlation analysis and the standard of correlation coefficient >0.3 and P value < .001,[ we established the association between ARGs and lncRNAs through mRNA-lncRNA co-expression network by the “limma” R package. We constructed and visualized the autophagy-related mRNA-lncRNA with prognostic value using the Cytoscape program and “ggalluvial” R package. Finally, a Sankey plot and co-expression network were performed.

Establishment of an arlncRNAs risk model to evaluate the risk score

According to P value <.001, univariate Cox regression analysis was first carried out. Multivariate Cox regression analysis based on the Akaike information criterion (AIC = 2182.17) search for the relative optimal fitting model. Depending on the optimal fitting model results, the arlncRNAs risk model was built and the risk score was calculated for each patient by the following formula.

coef (arlncRNAi) and expr (arlncRNAi), respectively represent the arlncRNAs-survival correlation coefficient and arlncRNAs expression of patients.

Melanoma patients were divided into high- and low-risk groups by the median risk score. The 2 groups’ performances were evaluated by Kaplan–Meier (K–M) analysis and the time-dependent receiver operating characteristic (ROC) curves.

Independent prognostic analysis

The clinical characteristics data from TCGA, including age, gender, TNM stage, and clinical-stage from the TCGA database, were put in order. Patients with complete clinical characteristics were included in independent prognostic analysis. Univariate and multivariate Cox regression analyses were performed to analyze the survival prognosis relationship with clinical characteristics and risk score. The P value <.05 was regarded as statistically significant. Finally, the violin plot and the ROC curve were performed.

Principal component analysis (PCA)

The core idea of PCA is the mapping of N-dimensional features to K-dimensional (K < N). In our study, the N dimensions represent the number of prognostic lncRNAs, K dimensions represent 2 and 3 dimensions, which can be visualized. Two- and three-dimensional PCA, respectively maps the point of prognostic gene expression level by PC1-PC2 and PC1-PC2-PC3. We assessed the prognostic model's validity by comparing PCA of all autophagy-related mRNAs, all arlncRNAs, and 9 prognostic arlncRNAs.

Gene set enrichment analysis (GSEA)

GSEA software (version 4.1.0 downloaded from https://www.gsea-msigdb.org/gsea/index.jsp) was used to explore Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between lncRNAs expression and dysregulation of tumor signaling pathways in high- and low-risk group. GSEA was performed by the gene sets “c2.cp.kegg.v7.2.symbols.gmt.” Normalized enrichment score (NES) was calculated, and nominal P value <.05 and false discovery rate q-value <0.25 were considered significant.

Results

Construction of the arlncRNAs prognostic model

A total of 457 melanoma patients with complete follow-up data were enrolled in our study. A total of 904 arlncRNAs were obtained by establishing an mRNA-lncRNA co-expression network (The entire set of 904 arlncRNA were shown in Table S1, Supplemental Digital Content), among which 31 arlncRNAs were significantly associated with the patients’ survival by Univariate Cox regression analysis (see Fig. 2A). Then, through multivariate Cox regression analysis, 9 prognostic arlncRNAs were screened by defining the AIC value =2182.17. They are USP30-AS1, LINC00665, PCED1B-AS1, LINC00324, LINC01871, ZEB1-AS1, LINC01527, AC018553. 1, and HLA-DQB1-AS1. A visual co-expression network (see Fig. 2B) and a Sankey plot (see Fig. 2C) of 9 autophagy-related mRNAs-lncRNAs were performed. Finally, an arlncRNAs prognostic model was established, and melanoma patients were divided into high- and low-risk group by the median risk score (median risk score = 1.0463).

Figure 2

Identification of prognostic arlncRNAs with significant value in melanoma. (A) The forest plot showed the Hazard ratio, 95% confidence interval and P value of prognosis arlncRNAs by univariate Cox proportional hazards analysis. (B and C) co-expression network and sankey of 9 autophagy-related mRNAs-lncRNAs. ArlncRNAs = autophagy-related long non-coding RNAs.

The prognostic value of the arlncRNAs risk model

Risk curves (see Fig. 3A) and scatter plots (see Fig. 3B) were used to illustrate the relationship between arlncRNAs prognostic model and patients’ survival status. Heat map (see Fig. 3C) shows that the expression of ZEB1-AS1, USP30-AS1, LINC00324, LINC01871, and HLA-DQB1-AS1 were high in the low-risk group and the expression level of LINC00665, LINC01527, AC018553. 1, and PCED1B-AS1 were low in the low-risk group. K–M analysis showed that patients in the low-risk group had a longer survival time than those in the high-risk group (P value = 9.631e-12, see Fig. 4A). The K–M analysis curves of the arlncRNA, USP30-AS1 (see Fig. 4B), LINC00665 (see Fig. 4C), PCED1B-AS1 (see Fig. 4D), LINC00324 (see Fig. 4E), LINC01871 (see Fig. 4F), ZEB1-AS1 (see Fig. 4G), LINC01527 (see Fig. 4H), AC018553. 1 (see Fig. 4I), and HLA-DQB1-AS1 (see Fig. 4J), were presented, respectively. A total of 351 patients with complete clinical characteristics were enrolled in an independent prognostic analysis. The clinical features are detailed in Table S2, Supplemental Digital Content. Combined with the clinicopathological characteristics of the patients, the time-dependent ROC curve (Fig. 4A) of arlncRNAs risk model showed that the value of area under curve (AUC) was between 0.7 and 0.8. The hazard ratio (HR) of the risk score and 95% confidence interval (CI) were 1.586 and 1.438–1.750 (P value <.001) in univariate Cox regression analysis (see Fig. 5B), and 1.535 and 1.374–1.716 (P value <.001) in multivariate Cox regression analysis (see Fig. 5C). The multi-index ROC curve (see Fig. 5D) showed that the AUC value of arlncRNAs risk model was 0.743, which was higher than other clinical characteristics.

Figure 3

The prognostic value of the arlncRNAs risk model. (A) The risk score of 457 patients. (B) The scatter plot of each patients’ survival status. (C) The heatmap showed the expression levels of 9 arlncRNAs in the high- and low-risk groups. ArlncRNAs = autophagy-related long non-coding RNAs.

Figure 4

Kaplan–Meier survival analysis. A, Kaplan–Meier survival analysis of the high- and low-risk groups. B–J, Kaplan–Meier survival analysis of the 9 arlncRNAs’ expression, USP30-AS1 (B), LINC00665 (C), PCED1B-AS1 (D), LINC00324 (E), LINC01871 (F), ZEB1-AS1 (G), LINC01527 (H), AC018553. 1 (I), and HLA-DQB1-AS1 (J), respectively. ArlncRNAs = autophagy-related long non-coding RNAs.

Figure 5

Assessment of the prognostic arlncRNAs risk model. (A) the time-dependent ROC curve of arlncRNAs risk model. (B) univariate Cox regression analyses of arlncRNAs risk models and melanoma patients’ clinical characteristics. (C) the ROC curves of arlncRNAs risk models and melanoma patients’ clinical characteristics. (D) multivariate Cox regression analysis. ArlncRNAs = autophagy-related long non-coding RNAs; ROC = receiver operating characteristic curve.

Correlation of the expression of the 9 arlncRNAs with clinicopathological characteristics

As showed in Figure 6, 5 arlncRNAs (USP30-AS1, LINC00665, PCED1B-AS1, ZEB1-AS1, and LINC01527) were significantly associated with T stage, and 2 arlncRNAs (PCED1B-AS1 and LINC01527) were significantly correlated with patients’ age. None arlncRNAs had statistical significance with gender, N stage, M stage, and clinical stage (P value >.05).

Figure 6

The correlation of the expression of the nine arlncRNAs with patients’ clinical characteristics. (A) age (≤65 and >65). (B) gender (female and male). (C) T stage (T0, T1, T2, T3, and T4). (D) N stage (N0, N1, N2, and N3). € M stage (M0 and M1), (F) clinical stage (stage I-II and stage III-IV), ∗P < .05, ∗∗P < .01, ∗∗∗P < .001.

Two-dimensional PCA analysis was performed between low- and high-risk groups based on the expression of all autophagy-related mRNAs (see Fig. 7A), all arlncRNAs (see Fig. 7B), and 9 prognostic arlncRNAs (see Fig. 7C). Three-dimensional PCA analysis was performed between low- and high-risk groups based on the expression of all autophagy-related mRNAs (see Fig. 7D), arlncRNAs (see Fig. 7E), and 9 prognostic arlncRNAs (see Fig. 7F). The results showed that the distribution of low- and high-risk groups in the other 2 patterns was disorderly, and there was no clear rule to follow, while the nine arlncRNAs risk model could divide melanoma patients into two parts.

Figure 7

Principal component analysis (PCA). (A-C) two-dimensional PCA were performed between low- and high-risk groups based on the expression of autophagy-related mRNAs, all arlncRNAs and 9 prognostic arlncRNAs. (D-F) three-dimensional PCA were performed between low- and high-risk groups based on the expression of autophagy-related mRNAs, all arlncRNAs and 9 prognostic arlncRNAs. arlncRNAs = autophagy-related long non-coding RNAs, PCA = principal component analysis.

Five signaling pathways were enriched in the high-risk group (Fig. 8A). The top 20 signaling pathways with the highest absolute value of NES were enriched in the low-risk group (Fig. 8B,C). These pathways were closely related to the progression of autophagy and melanoma. The detailed data for the NES, nominal P value, and false discovery rate q-value were shown inTable S3, Supplemental Digital Content.

Figure 8

Gene set enrichment analysis. (A) 5 signaling pathways were enriched in the high-risk group. (B-C) the top 20 signaling pathways with the highest absolute value of NES were enriched in the low-risk group. NES = normalized enrichment score.

Discussion

Autophagy plays a double-edged role for tumor cells.[ In the early stage of tumor genesis, autophagy plays a major role in inhibiting tumor growth.[ Autophagy disorders can also cause genome-instability and inflammation, which promote tumor growth.[ Many shreds of evidence indicate that autophagy is activated, the development and progression of melanoma cells were affected.[ D’Arcangelo D and his colleagues distinguished three ARGs wipi1, BAG1, and pex3 from 5 different databases, and they found these genes can be used as melanoma biomarkers.[ Recently, many studies have proven that the autophagy process of tumor cells can be affected by autophagy-related lncRNAs via to lncRNA-microRNA axis.[ Similarly, autophagy-related lncRNAs have also been found in melanoma cells and mediate the autophagy process. Peng et al found that lncRNA-ZNNT1 was over-expressed in Uveal melanoma cells and promoted autophagy by up-regulating autophagy-related gene 12, meanwhile the knockdown of ZNNT1 weakens autophagy induced by PP242.[ In summary, the arlncRNA mediates the autophagy process of melanoma through different pathways and affects the prognosis of melanoma patients. Generally, the above studies suggest that ARGs and arlncRNAs can affect the prognosis of melanoma patients, and many regulatory pathways of arlncRNAs have been validated. However, the role of arlncRNAs in melanoma patients remains unclear. Therefore, we comprehensively analyzed arlncRNAs in melanoma patients from the TCGA.

In this study, a total of 9 arlncRNAs (USP30-AS1, LINC00665, PCED1B-AS1, LINC00324, LINC01871, ZEB1-AS1, LINC01527, AC018553. 1, and HLA-DQB1-AS1) were identified as prognostic biomarkers in melanoma. We further incorporated the nine arlncRNAs into a novel signature parameter for evaluating prognoses of melanoma patients. Patients were divided into high- and low-risk groups according to the medium risk score of the arlncRNA risk model. K–M analysis showed that patients with high- and low-risk groups had significant differences in survival curves (P value = 9.631e–12). Furthermore, the novel parameter composed of nine arlncRNAs was indicated an independent prognostic factor by multivariate Cox analysis (P value <.001). We further verified the predictive effect of the novel parameter by the time-dependent ROC curve, and the results indicated that the arlncRNAs risk model had a high predictive performance.

So far, only 4 arlncRNAs (LINC00665,[ PCED1B-AS1,[ LINC00324,[ and ZEB1-AS1[) have been biologically validated in other tumors. However, as far as we know, these 9 arlncRNAs have not been fully validated in melanoma. Therefore, the 9 arlncRNAs have great potential and research value in exploring autophagy-related therapeutic targets and studying autophagy-related signaling pathways. We performed GSEA enrichment analysis and KEGG pathway analysis to explore the possible role of the arlncRNAs in different signaling pathways. The results revealed certain differences in high- and low-risk groups’ signaling pathways. Signaling pathways related to cell metabolism were enriched in the high-risk group, while autophagy- and immunity-pathways were mainly enriched in the low-risk group. KEGG pathway analysis reflected that patients with high metabolism of melanoma cells had a poor prognosis, while enhancement of immunity and autophagy improved the prognosis of melanoma patients. These results are consistent with current understanding; that is, stronger immunity and autophagy improve the prognosis of melanoma patients, while high metabolism suggests a poor prognosis.[

More interestingly, the nine arlncRNAs were not associated with the N stage and M stage, which indicated that lncRNAs did not affect the prognosis by impacting lymph node metastasis and distant metastasis. We were surprised to find that 5 prognostic lncRNAs (USP30-AS1, LINC00665, PCED1B-AS1, ZEB1-AS1, and LINC01527) expressed differently in different T stages of melanoma patients. It indicated that autophagy was closely related to the depth of tumor invasion in melanoma patients. Wang et al showed that autophagy could inhibit the epithelial-mesenchymal transition (EMT) of RAS mutant tumor cells.[ As a developmental process, EMT makes the quiescent cells lose their epithelial features and obtain their interstitial phenotypes, which leads to tumor cell activation and infiltration.[ Therefore, combined with our conclusion, we speculate that autophagy affects EMT, thus affecting the invasion depth of melanoma cells. However, this conjecture still needs to move forward a single step by biological verification.

The current study has several limitations. First of all, our research needs to be further validated in other independent databases to determine the arlncRNAs risk model's effectiveness. Second, the clinical information downloaded is limited and incomplete. Therefore, some patients have to be excluded, and the reduction of sample size may cause some bias in statistical analysis. At the same time, clinical information is absent, such as surgery, radiotherapy, and chemotherapy. Finally, the molecular mechanism of arlncRNAs, which regulates the occurrence and development of melanoma, needs to be further studied in order to promote the clinical application of phagocytic therapy.

Conclusion

In conclusion, our study shows that prognostic arlncRNAs can accurately predict melanoma patients’ survival outcomes. We established nine arlncRNAs risk models, dividing the population into high- and low-risk groups. We also proved that the risk model could be an independent prognostic factor for melanoma patients. Our study also shows that these nine arlncRNAs are potential prognostic and diagnostic biomarkers and promising targets for the treatment of melanoma.

Acknowledgments

We are very grateful for the contributions of TCGA 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.

Formal analysis: Haiying Cao.

Funding acquisition: Yu Jin.

Investigation: Zhehong Li, Yu Jin.

Methodology: Zhehong Li.

Project administration: Honghong Zheng.

Resources: Honghong Zheng.

Software: Zhehong Li, Junqiang Wei, Yafang Zhang, Mingze Song, Haiying Cao.

Validation: Junqiang Wei, Honghong Zheng, Yafang Zhang, Mingze Song, Haiying Cao, Yu Jin.

Visualization: Zhehong Li, Junqiang Wei.

Writing – original draft: Zhehong Li, Junqiang Wei.

Writing – review & editing: Yu Jin.