Development and validation of a nomogram with an autophagy-related gene risk model for predicting survival in patients with melanoma

Background: Melanoma is the third most common skin malignant tumor in the clinic, with high morbidity and mortality. Autophagy plays an important role in the development and progression of melanoma. We aimed to establish an autophagy-related genes(ARGs) expression based risk model for individualized prognosis prediction in patients with melanoma. Methods: Differentially expressed autophagy-related genes (DEARGs) in melanoma and normal skin samples were screened using TCGA and GTEx database. These DEARGs were used to perform KEGG functional enrichment analysis and GO analysis. Univariate and multivariate Cox regression analyses were performed on DEARGs to identify the optimal prognosis-related genes. These prognosis-related DEARGs were used to construct a risk score model, and the predictive effect of this risk model on the prognosis of melanoma patients was tested by the Kaplan-Meier curve, log-rank test, and ROC curve. Method of univariate and multivariate analysis were used to confirmed that the risk model of independent predictive value relative to other clinical variables, and build a nomogram based on the independent prognostic factors in the univariate analysis to predict overall survival(OS) in patients with melanoma, we used internal validation and calculation of concordance index (C-index) to test prediction effect of the nomogram. We also used the t-test to analyze the relationship between risk factors (risk genes and risk score) and clinical variables in the risk model. Results: We screened and finally obtained 6 optimal DEARGs (risk gene) through univariate and multivariate Cox analysis to construct the risk model: EIF2AK2(HR=0.403, P=0.007), IFNG(HR=0.659, P=0.003), DAPK2(HR=0.441, P=0.022), PTK6(HR=1.609, P=6.04E-05), BIRC5(HR=2.479, P=0.001), and EGFR(HR=1.474, P=0.004) were selected to establish the prognostic risk score model and validated The results of GO analysis showed that the gene function of the was concentrated in the functions of gland morphogenesis, protein insertion into membrane, and autophagy. The results of KEGG enrichment analysis showed that the function of the DEARGs was concentrated in the autophagy–animal, p53 signaling pathway, and platinum drug resistance. Kaplan-Meier survival analysis demonstrated that patients with high risk scores had significantly poorer overall survival (OS, log-rank P=6.402E−11). The model was identified as an independent prognostic factor. Finally, a prognostic nomogram including the risk model, T-stage, N-stage, and radiotherapy was constructed, and the calibration plots indicated its excellent predictive performance. Conclusion: The autophagy-related six-gene risk score model could be a prognostic biomarker and suggest therapeutic targets for melanoma. The prognostic nomogram could help individualized survival prediction and improve treatment strategies.

The results of GO enrichment analysis showed that the gene function of the DEARGs was concentrated in the functions of gland morphogenesis, protein insertion into membrane, and autophagy. The results of KEGG enrichment analysis showed that the function of the DEARGs was concentrated in the autophagy-animal, p53 signaling pathway, and platinum drug resistance. Kaplan-Meier survival analysis demonstrated that patients with high risk scores had significantly poorer overall survival (OS, log-rank P=6.402E−11). The model was identified as an independent prognostic factor. Finally, a prognostic nomogram including the risk model, T-stage, N-stage, and radiotherapy was constructed, and the calibration plots indicated its excellent predictive performance.

Conclusion:
The autophagy-related six-gene risk score model could be a prognostic biomarker and suggest therapeutic targets for melanoma. The prognostic nomogram could help individualized survival prediction and improve treatment strategies.

Background
Melanoma is the third most common skin malignant tumor in the clinic, with high morbidity and mortality. There were about 87,000 new cases and 10,000 deaths in the United States in 2017 [1,2]. Melanoma in the early stage can be cured by surgical resection and get longer overall survival(OS), the effect of simple surgical removal in melanoma patients in the middle-late stage is not ideal. Patients with tumor thickness greater than 4 mm(T4N0M0, AJCC stage Ⅱ B~C), regional lymph node metastasis or intermediate metastasis (T1-4N1-3M0, AJCC stage Ⅲ), even after complete surgical resection, had a recurrence rate of 60% ~ 75% and a 5-year OS rate of only 30% ~ 70% [3]. Although molecular targeted therapy and immune checkpoint therapy have shown a good therapeutic effect in metastatic melanoma [4][5][6], it is particularly urgent to further study the molecular mechanism of the occurrence and progression of melanoma and find effective therapeutic targets. Over the past decade, numerous studies have identified new sensitive and effective biomarkers for melanoma.
Autophagy not only plays an important role in the physiological processes such as component renewal, growth, and differentiation of cells, but also acts as a defense mechanism for cells to undergo microenvironment changes due to injury factors, and is involved in the pathological process of various diseases including tumors. In recent years, the research on autophagy and tumor has become a hot topic [7][8][9]. It is currently believed that, on the one hand, autophagy acts as a protective mechanism, so that tumor cells can avoid the damage caused by low nutrition, ionizing radiation and chemotherapy and survive continuously. On the other hand, excessive autophagy effect will cause autophagy sexual death (type Ⅱ programmed death), thus inhibiting the occurrence and development of tumor [10].
In our study, by analyzing the dataset from TCGA database, we aim to study and verify expression characteristics of autophagy-related genes(ARGs), in order to predict the prognosis of melanoma and therapeutic targets in addition, we combined with the risk model(made up of the optimal ARGs) and clinical parameters, establish a reliable nomogram to predict the prognosis of patients with melanoma, this model is more accurate than clinical risk factors predicting ability.

Screening of differentially expressed genes(DEGs)
To download the transcriptome data of 470 patients with melanoma from The cancer genome atlas (TCGA) portal website (https://portal.gdc.cancer.gov/), and clinical information of 470 patients, to obtain the transcriptome data of 44 cases of normal skin tissue from Genotype-Tissue Expression (GTEx) database (https://xenabrowser.net/). All data are processed using R software (https://www.r-project.org/), the clinical data of patients with melanoma is as follows ( Table 1). Transcriptome data from the TCGA database of melanoma samples and transcriptome data from normal skin tissue samples from the GTEx database were merged into the new data matrix. Differential expression analysis was conducted using the "edgeR" package. DEGs were identified in melanoma samples compared to normal samples. Adjusted P value < 0.05 and | log 2 fold change (log2 FC) | >1 were used for the threshold to identify DEGs.

Identification of differentially expressed autophagy-related genes (DEARGs) and enrichment analysis
The 222 ARGs were obtained from the Human Autophagy Database (HADb, http://www.autophagy.lu/), which provides a complete and up-to-date list of human genes and proteins involved directly or indirectly in autophagy. DEARG was screened from DEGs. Then, the "clusterProfiler" package from Bioconductor was used to do the Kyoto

Encyclopedia of Genes and Genomes(KEGG) functional enrichment analysis and Gene
Ontology(GO) analysis of the DEARGs. P-value of < 0.05 was considered statistically significant.

Identification of prognostic DEARGs
In order to screen out the DEARGs associated with the prognosis of melanoma patients, we used univariate Cox regression to screen out the DEARGs. P-value <0.05 was used as the exclusion standard in our study.

Identify prognostic DEARGs for inclusion in the risk model
We thus analyzed the DEARGs associated with the prognosis of melanoma patients via multivariate Cox proportional hazards regression analysis (with forwarding selection and backward selection). Finally, we obtained the optimal DEARGs(risk genes) to be incorporated into the prognostic risk model.

Construction of a risk model in melanoma patients
We used the expression levels of mRNA and estimated regression coefficients of the risk genes to calculate a risk score for melanoma patients. Patients with incomplete clinicalpathological information were excluded. The risk score model was established with the following formula: Risk score = expression level of Gene 1 * β 1 + expression level of Gene 2 * β 2 +…+ expression level of Gene n * β n ; where β is the estimated regression coefficients calculated by the multivariate Cox regression model.
The risk model was used to measure the prognostic risk of each melanoma patient. The median risk score of the melanoma cohort was used as the cut-off value to divide all the melanoma patients into two groups: the high-risk group and the low-risk group. A high risk score indicates a poor prognosis for melanoma patients.
Then, a Kaplan-Meier survival curve was constructed to estimate the prognosis of patients with high-risk scores or low-risk scores, and the survival differences between the high-risk and low-risk groups were assessed by a two-sided log-rank test. The prognostic performance was evaluated by using time-dependent receiver operating characteristic (ROC) curve analysis within 1, 3, and 5 years to evaluate the predictive accuracy of the ARG-based risk model.

Independent prognostic value of the risk model in the melanoma patients
Furthermore, to determine whether the predictive power of the risk model could be independent of other clinical-pathological parameters(including age, gender, pharmacotherapy, radiotherapy, AJCC stage, T-stage, and N-stage) for patients with melanoma, univariate and multivariate Cox proportional hazards regression analyses were performed in the melanoma cohort. We exclude M-stage because patients in M1-stage were too few.
To better predict the 1-year, 3-year and 5-year OS of melanoma patients, we used variables screened by univariate analysis to construct a nomogram. We then internally examine the nomogram and calculate the concordance index(C-index). C-index and AUC values of 0.75 or greater were considered to have excellent predictive value, and values of 0.6 or greater were considered acceptable for survival predictions [11].

The relationship between risk model and clinical characteristics
To examine the relationship between our model and the clinical characteristics of melanoma, we analyzed the relationship between the risk factors from our model (risk score and risk genes) and the clinical characteristics of the entire melanoma cohort (including age, gender, pharmacotherapy, radiotherapy, AJCC stage, T-stage, and Nstage).

Screening for differentially DEGs
We screened 5432 up-regulated genes and 2548 down-regulated genes, a total of 7980 DEGs for melanoma cancer from the TCGA data set, the volcano and heat map is shown in the appendix.

Identification of DEARGs and enrichment analysis
This analysis revealed 31 DEARGs, including 19 genes that were upregulated and 12 genes that were downregulated in melanoma tissues compared with normal skin tissues (Figure1A and Figure1B). The results of GO enrichment analysis showed that the gene function was concentrated in the functions of gland morphogenesis, protein insertion into membrane, and autophagy. The results of KEGG enrichment analysis showed that the function of the DEARGs was concentrated in the autophagy -animal, p53 signaling pathway, and platinum drug resistance (Figure2A and Figure2B).

Identification of prognostic DEARGs
We screened 17 DEARGs associated with the prognosis of melanoma patients through the Kaplan-Meier curves and the log-rank test(P < 0.05) (Figure 3). , and EGFR were identified as high-risk genes (predicting a poor prognosis), while EIF2AK2, IFNG, and DAPK2 were identified as a low-risk gene (serving as a protective factor) in terms of the OS of melanoma patients.

Construction of the prognostic risk model in the melanoma cohort
We used the mRNA levels and estimated regression coefficients of the risk genes to Patients were divided into the high-risk group (n = 223) and the low-risk group (n = 224) according to the median risk score. To determine the difference in prognosis between high-risk and low-risk groups, we created a Kaplan-Meier curve based on the log-rank test.
The prognosis of the high-risk group was worse than that of the low-risk group (p < 0.05) ( Figure 4A). The 3-year and 5-year OS in the high-risk group was 61.5% and 44.8%, respectively, while the corresponding OS in the low-risk group was 81.4% and 72.8%, respectively. We then used a time-dependent ROC curve to test the predictive accuracy of the model. The ROC (AUC) value of the predictive model in 1-year, 3-year and 5-year were 0.752, 0.691, 0.719, respectively (Figure 4B, Figure 4C and Figure 4D). We then ranked melanoma patients' risk scores and analyzed their distribution (Figure 5A). The survival status of each patient in the melanoma cohort is shown in the dot diagram in Figure 5B.
The heat map was used to describe the expression patterns of risk genes in both prognostic groups (Figure 5C).

Independent prognostic value of the risk model in the entire TCGA cohort
Next, univariate and multivariate Cox regression analyses were performed to assess whether the risk score generated by our model could be independent of other clinical characteristics (age, gender, pharmacotherapy, radiotherapy, AJCC stage, T-stage, and Nstage) as prognostic factors for melanoma. Univariate analysis showed that age, AJCCstage, T-stage, N-stage, and risk score were significantly correlated with the prognosis of melanoma patients (Figure 6A). Multivariate analysis showed that T-stage, N-stage, radiotherapy, and risk score were significantly correlated with the prognosis of melanoma patients (p < 0.05)( Figure 6B). These results indicate that the prognosis of risk models can independently predict the level of the prognosis of patients with melanoma.
Then, T-stage(p < 0.001), N-stage(p < 0.001), radiotherapy (p = 0.025), and risk score(p < 0.001) were included in multivariate analysis( Figure 6C). In order to better predict melanoma patients for 1-year, 3-year, and 5-year survival rates, we built a nomogram according to the OS related variables(T-stage, N-stage, radiotherapy and risk score) ( Figure 6D). Internal validation and C-index calculation were used to evaluate the accuracy of the nomogram. The results showed that the nomogram could accurately predict OS at1-year, 3-year and 5-year survival for melanoma patients, with a C-index of 0.722(p = 0.022) (Figure 7A, Figure 7B and Figure 7C).

The relationship between risk model and clinical characteristics
To test our model's ability to predict melanoma progression, we analyzed the relationship between the risk factors from our model (risk score and risk genes) and the clinical characteristics of melanoma patients (age, gender, pharmacotherapy, radiotherapy, AJCC stage, T-stage, and N-stage). We used the X-tile software to determine the optimal cutoff age of 75 (Appendix). The results showed that among these risk factors, the expression values of EIF2AK2 and IFNG decreased with increasing age, the risk scores do the

Discussion
Many studies have shown that autophagy plays an important role in the formation, metabolism and drug resistance of tumors [12]. Correcting the imbalance of the autophagy pathway can inhibit tumor growth and improve the sensitivity to various treatments.
Studies on autophagy and melanoma have been reported, and clinical trials of drugs related to autophagy have also been reported [13][14][15]. The present findings suggest that autophagy plays a driving role in the development and progression of melanoma [16].
There is no doubt that the exploration and discovery of new ARGs will find more new biomarkers and therapeutic targets for melanoma.
In our study, we first identified 31 DEARGs based on the TCGA and GTEx database and then confirmed 6 optimal genes significantly correlated with prognosis in univariate and multivariate Cox regression analyses. The results of gene enrichment analysis showed that the gene function was concentrated in the functions of gland morphogenesis, protein insertion into membrane, and autophagy in GO enrichment analysis. The results of KEGG enrichment analysis showed that the function of the DEARGs was concentrated in the autophagy-animal, p53 signaling pathway, and platinum drug resistance. These results reveal the pathways in which these DEARGs are enriched in melanoma and provide a basis for us to better study the molecular mechanism of melanoma.
Next, we constructed a risk score model based on 6 DEARGs significantly associated with the prognosis of melanoma. We further verify the reliability and stability of the model. The results showed that the model could accurately distinguish patients with different survival outcomes. Univariate and multivariate Cox regression analysis showed that the model could independently predict the prognosis of melanoma patients. The result of the Nomogram showed that combining the model with other clinical characteristics (T-stage, N-stage, radiotherapy and risk score) was also an accurate predictor of the outcome of melanoma patients. Therefore, our model can be used to identify melanoma patients at high risk of death, enabling clinicians to improve patient outcomes through early intervention in their clinical work.
We further analyzed the relationship between factors in our model and clinical characteristics (age, gender, pharmacotherapy, radiotherapy, AJCC stage, T-stage, and Nstage). We found that the various factors in the model were significantly correlated with melanoma's progress. Therefore, our model has high clinical applicability in predicting the development of melanoma. Among these genes, except for EGFR gene related gene detection and drugs have been used in various tumor diseases [17], other genes obviously have great potential to become new melanoma biomarkers and therapeutic targets.
Previous bioinformatics studies in melanoma have been reported in a long time. Hai-zhou Wang et al [18] identified that Plakophilin1 gene is associated with the metastasis in However, there is no bioinformatics analysis yet exploring the relationship between ARGs and melanoma. Compared with the above studies, our research has absolute superiority.
Firstly, we use the TCGA database and the GTEx database to screen DEGs, which is an innovation. Secondly, we used a variety of statistical methods to screen out the genes related to autophagy and significantly related to the prognosis of melanoma patients. We also constructed the nomogram of the clinical characteristics and the risk model to predict the 3-year and 5-year survival rates of melanoma patients. Of course, our study has its flaws, we need further basic tests to verify our conclusions.   Tables   Table 1 Clinical information of the melanoma patients