Co-expression analysis of autophagy genes
We downloaded a total of 232 ARGs from the HADb database. In addition, by screening the mRNA matrix of TCGA, we obtained the mRNA expression matrix related to melanoma autophagy. Furthermore, we set |Correlation coefficient| ≥5 and p-value <0.001 as the cut-off standards, and used the limma package of the R software to analyze the co-expression of the mRNA matrix and lncRNA matrix of melanoma. Finally, we screened 89 mRNAs and 100 lncRNAs in patients with melanoma (Supplementary Figure S1).
Survival analysis of lncRNAs associated with autophagy
Prior to the multivariate COX regression analysis, we performed a univariate Cox survival analysis on 100 lncRNAs. Finally, the Kaplan–Meier survival analysis showed that 39 lncRNAs were closely related to the prognosis of patients with melanoma. The high expression of only two lncRNAs was significantly correlated with shorter overall survival (OS) in patients with melanoma, whereas the remaining 37 lncRNAs did not show such a relationship.
Multivariate Cox regression model analysis
The combination of multiple genes can better predict the prognosis of patients. Thus, we conducted a stepwise COX regression analysis of 37 lncRNA with survival significance, and finally obtained a gene combination (Figure 1, Table 1) composed of seven lncRNAs (HLA-DQB1 antisense RNA 1 [HLA-DQB1-AS1], USP30 antisense RNA 1 [USP30-AS1], AL645929, AL365361, long intergenic non-protein coding RNA 520 [LINC00520], LINC00324, and AC055822). In the multivariate COX regression analysis, patients were divided into high- and low-risk groups according to the median value of their riskScore. According to the riskScore, we used the survivalROC package to draw the survival curve (Figure 2A) and ROC curve (Figure 2B) of this gene combination. Survival results showed that the high expression of this gene combination was significantly associated with poor OS in patients at high risk. The AUC of 0.742 suggests that this gene combination may be used to predict the prognosis of patients with melanoma. The riskScore results showed that the high-risk group had a higher risk coefficient than the low-risk group, rising from left to right as illustrated in Figure 2C. Compared with the low-risk group, the survival diagram results showed that the high-risk group was characterized by more deaths and shorter survival time (Figure 2D). Finally, the expression of seven genes in the high- and low-risk groups was visualized using a heat map (Figure 2E).
In addition, to verify the higher accuracy of the risk model for predicting the prognosis of patients versus clinical information, we first performed a univariate COX analysis based on clinical information, including sex, stage, T, M, and N. The forest plot results of the univariate analysis showed that stage, T, M, N, and riskScore could be used to predict the prognosis of patients (Figure 3A). However, we used a multivariate COX regression analysis and ROC curve to draw the clinical information and determine the index that may be used as the best independent prognostic factor. The forest plot results of the multivariate COX analysis showed that T, N, and riskScore could be used as independent prognostic factors (Figure 3B). However, the ROC results showed that the AUC of riskScore, T, and N was 0.734, 0.707, and 0.658, respectively (Figure 3C). Therefore, the risk model we constructed is more accurate than T and N in predicting patient survival, and is better than other traits (Table 2).
Construction of the co-expression network
We constructed an autophagy-related co-expression network of seven lncRNAs and 16 mRNAs using the Cytoscape software (Figure 4A, Table 3). HR >1 indicated a risk factor, whereas HR <1 indicated a protective factor. We found that, in the co-expression network, only LINC00520 was a risk factor, whereas the other six lincRNAs were protective factors. Of note, the Sankey diagram yielded consistent results (Figure 4B).
GSEA in melanoma
We used GSEA to analyze the GO and KEGG, and study the potential pathogenic mechanism of patients at high and low risks of melanoma. The results of the GO enrichment analysis showed that the BP of the high-risk group was mainly related to the metabolism of melanin. The GO of patients at low risk was related to the positive regulation of autophagy and the activation of the immune response (Figure 5A). The results of KEGG were related to the regulation of autophagy, antigen processing and presentation, and T-cell receptor signaling pathway (Figure 5B). Therefore, we hypothesized that the potential pathogenesis of melanoma in patients at high and low risks is related to the regulation of autophagy and the immune system. It was further demonstrated that the immune response plays an important role in the autophagy of melanoma.