Correlation of Autophagy-Related Genes in High Risk and Low Risk Signatures Based on Prognostic Index Model SKCM Patient

Autophagy plays a dual role in tumor development and autophagy-related genes (ARGs) involved in the development of cancer. However, the correlation between these high and low risk ARGs is still unclear. To systematically study the relationship between ARGs and melanoma patients, the expression proles of ARGs were integratedly analyzed based on the TCGA and GTEx dataset. The results suggested 8-ARGs marker can predict the prognosis of skin cutaneous melanoma (SKCM), 5-ARGs (CFLAR, DAPK2, ITGA6, DNAJB9 and RGS19) were low risk index, 3-ARGs (EIF2AK2, EGFR and PTK6) were high risk index. Further validation showed CFLAR, DNAJB9 and PTK6 were signicant related to SKCM. In tumor immunization, CFLAR and DNAJB9 have stronger correlation with tumor immune cells, and there is a signicant positive correlation among the low risk genes. However, there is no signicant correlation between the high and low risk factors. In addition, our results also showed that CFLAR and PTK6 were signicantly negatively correlated with tumor purity; PTK6 was signicantly positively correlated with patient's age, CFLAR was signicantly negatively correlated with patient's age; CFLAR and DNAJB9 were signicantly positively correlated in purity and age. It provides a new prognostic indicator for patients and a new idea for the mechanism of autophagy in melanoma. Expression prole and prognostic value of ARGs. genes, analysis revealed the biological processed and molecular functions involved 36 prognostic-related signaling involved 36 is the most


Introduction
As one of most aggressive type of cancers worldwide, Skin cutaneous melanoma (SKCM) is originated from melanocytes [1][2][3][4][5]. SKCM shows clear differences in incidence, mortality, genomic pro le, and anatomic presentation, depending on the country of residence, ethnicity, and socioeconomic status [6]. Cutaneous melanoma accounts for approximately only 5% of all skin malignancies, but is thought to be the most invasive and lethal form [6]. Once malignant melanoma metastasis occurs, it develops faster than all other solid tumors [7]. Despite recent advances in curative surgery and adjuvant therapy, the mortality is still high due to its late detection and high recurrence rates [8]. Therefore, there is an urgent need to nd new molecular markers to predict the prognosis of SKCM patients.
Autophagy is a multi-step lysosomal degradation process, which is initiated from the formation of autophagosome to deliver the protein and cell organelle to lysosome for degradation. This catabolic process is involved in a variety of autophagy-related genes (ARGs), which has been extensively studied and been proven involved in the development of cancer [9][10][11][12][13]. However, autophagy is suggested to have both tumor-suppressing and tumor-promoting functions during cancer progression, depending on the type and stage of tumor [14]. The role of autophagy in tumorigenesis is still controversial, and its mechanism is still imperfect and uncertain. Autophagy has a complex function in the development of tumor. It is an urgent problem to further explore the potential biological process of autophagy in tumorigenesis and development [15].
A large number of studies have shown that autophagy is also essential for maintaining cellular homeostasis in the skin and reported a correlation between autophagy and SKCM [16][17][18][19]. For example, CX-F9 inhibits the malignant phenotype of melanoma cells and down-regulation of ATG5 contributes to tumorigenesis in early-stage cutaneous melanoma [18,20]. Modulation of autophagy contribute to the regulation of melanocyte biology and skin pigmentation [17]. Recently research showed 7-DEARG signature could be regarded as an independent prognostic signature in clinical setting, and nomogram may supply a more simple and accurate prediction for the prognosis of melanoma [21]. However, there is rarely study revealed correlation between these prognostic signature ARGs. The purpose of this study was to gain insight into the potential clinical utility of ARGs for prognostic strati cation and to investigate the correlation and characteristics of ARGs.
According to the tumor lymph node metastasis (TNM) classi cation system, pathological staging is commonly used to predict the survival rate of patients [22,23]. Other factors, such as age, performance status, and tumor location, also affect patient survival [24,25]. In recent years, nomograms have gradually gained popularity, which use different clinical variables to determine a statistical prognostic model that generates a probability of clinical outcomes for an individual patient [15,26]. Nomograms have been applied in various types of cancers to establish a prognostic nomogram, which can help in prediction of early bone metastases for patients with breast cancer [27]. Therefore, our research obtained prognostic ARGs index by combining TNM and nomograms system model, and explored the correlation signature through GEPIA and TIMER database, which facilitate the development of personalized prognostic information for SKCM patients.

Result
Identi cation of prognostic ARGs in melanoma RNA-seq and clinical data from 472 SKCM tissue samples and 1 non-tumor samples were downloaded from TCGA. Only 210 ARGs showed expression values in TCGA database. Finally obtained 49 upregulated and 45 down-regulated ARGs according to FDR<0.05 and [log2(fold change)]>1 ( Figure S1A and 1B) and a scatter plot was visualized to show the expression pattern of 94 differentially expressed ARGs between melanoma and non-tumor tissue ( Figure S1C). In addition, functional enrichment analysis of 94 differentially expressed ARGs provides a biological understanding of these genes. GO enrichment shows that the biological process of differential genes is mainly involved in autophagy and process utilizing autophagic mechanism, the molecular functions is mainly involved in ubiquitin protein ligase binding ( Figure S1D). KEGG enrichment shows that pathways of differential genes mainly involved pathways in autophagy, apoptosis, shigellosis and human papillomavirus infection.

Establishment of Autophagy-Related Signature in Training Set and Validation Set
Patients in TCGA dataset was randomly assigned in a 5:5 ratio to training set and validation set with the same proportion of each SKCM stage. 36 differentially expressed ARGs were initially subjected to univaritate Cox proportional hazards regression analysis in the training set and validation set. We constructed autophagy prognostic index (API) to divide SKCM patients into two groups (high-risk and low-risk) with discrete clinical outcomes for overall survival (OS). The result showed 8 ARGs were signi cantly associated with the OS (P<0.05). These include 2 risky genes and 6 potential protective genes. Figure 2 showed autophagy-related signature in training set and validation set. Distribution of prognostic index in TCGA dataset was showed in Figure 2A, survival status of patients in different groups was showed in Figure 2B and heatmap of the expression pro le of the included ARGs was showed in Figure 2C. To determine the performance of the API in predicting clinical outcomes in skin cutaneous melanoma patients, Kaplan-Meier survival curves were plotted to analyze different survival times between high-ris and low-risk groups. Kaplan-Meier analysis showed that patients with high-risk group was signi cantly lower than that in the low-risk group ( Figure 2D).
Univariate analysis (HR=1.186, 95%, CI=1.121-1.254, P<0.001) showed that ARI was signi cantly associated with patient prognosis ( Figure S2A). In addition, after adjusting for clinicopathological features such as age, tumor subtype, tumor size, and lymph node metastasis, API remined an independent prognostic indicator for SKCM patients in multivariate analysis (HR=1.194, 95%, CI=1.129-  Figure S3C). A weighted total score calculated from each variable was used to estimate the 3-and 5-year OS of patient with SKCM ( Figure S3D). The result showed gender, age, stage and risk score were signi cant risk factors for patients with SKCM (P<0.05).
The calibration plots showed well correlation between observed OS and nomogram predicted OS ( Figure  S3Eand S3F).
The survival curve of single gene showed that high expression of CFLAR, DAPK2, EIF2AK2, ITGA6, DNAJB9 and RGS19 were signi cantly related to the increasing survival rate of patients, and low expression of EGFR and PTK6 were closely related to clinical prognosis ( Figure 3A). Overall survival analysis of GEPIA showed expression level of RGS19, EIF2AK2 and EGFR had not related to patient survival ( Figure 3B). Besides, stage pro le suggested that only CFLAR, DNAJB9 and PTK6 had clinical stage relevance ( Figure 3C).

Relationship between the expression of targeted ARGs and immunity in SKCM
To reveal the relationship between the expression level of signi cant ARGs in melanoma and immunity, we downloaded the score data of six kinds of immune in ltrating cells of SKCM from timer, and respectively analyzed the correlation between the expression of CFLAR, DNAJB9 and PTK6 and the proportion of these cells. Our results showed that CFLAR was positively correlated with the scores of six kinds of immune. Among them, cells_ Neutrophil had the highest correlation with CFLAR and B_Cell correlation was the lowest. DNAJB9 showed signi cant positive correlation with B_cell, CD8_Tcell, Timer-Neutrophil, Timer_Marcrophage and Dendritic. The correlation of DNAJB9 with Timer_neutrophil was the highest. PTK6 only signi cantly negative correlated to B_ Cell ( Figure 4A). In addition, the results showed that only CFLAR was signi cantly positively correlated with the overall microenvironment score, CFLAR and DNAJB9 were signi cantly positively correlated with immune score, CFLAR, DNAJB9 and PTK6 were signi cantly correlated with stromal score, among which CFLAR and DNAJB9 were positively correlated and PTK6 was negatively correlated ( Figure 4B). Finally, considering the importance of neoantigen for personalized treatment of cancer patients, we further analyzed the correlation between CFLAR, DNAJB9, PTK6 and neoantigen, but unfortunately, none of the three genes had signi cant correlation ( Figure 4C). SCNA module of target genes were analyzed by TIMER, which provides the comparison of tumor in ltration levels among tumors with different somatic copy number alterations for a given gene. The results showed that the copy number of arm-level deletion of CFLAR was related to the in ltration of B cell, CD4+T cell, macrophage, and dendritic cell, and the copy number of arm level gain was related to the in ltration of CD4 + T cell. The copy number of DNAJB9 arm level deletion was correlated with B cell, CD4+T cell and dendritic cell in ltration. Copy number of PTK6 high ampli cation was signi cantly correlated with six kind immune cell lines, copy number of arm level gain was correlated with B cell, CD4+T cell, macrophage, and dendritic cell in ltration, copy number of arm level deletion was correlated with B cell, CD8+T cell, CD4+T cell, and dendritic cell in ltration ( Figure S4A). Immune checkpoint genes also play an important role in tumor immunotherapy. Our result found that PTK6 only had a signi cant correlation with VTCN1, CFLAR was signi cantly correlated with most immune checkpoint genes, but not with TNFRSF14, CD276, VTCN1, HHLA2, CD70, TNFSF9, CD44. DNAJB9 signi cantly correlated with CD200, NRP1, TNFSF4, CD28, CD200R1, HAVCR2, CD80, CD160, TNFSF14, VSIR, CD86 and TNFRSF9. There was a similar correlation between DNAJB9 and CFLAR. While the correlation genes of high-risk signature PTK6 had almost opposite correlation characteristics with low-risk signature DNAJB9 and CFLAR ( Figure S4B).

Correlation between CFLAR, DNAJB9 and PTK6 in patients with SKCM
Firstly, we studied the correlation between high risk signature and low risk signature in tumor, skin-sun exposed, skin-not sun exposed groups by GEPIA. The result showed high risk signature only negatively correlated with low risk signature in the skin-sun exposed group ( Figure S5A). Then we analyze the correlation between CFLAR, DNAJB9 and PTK6 by TIMER. The results showed that in the early stage of SKCM, DNAJB9 was negatively correlated with PTK6 and positively correlated with CFLAR, in the late stage of SKCM, DNAJB9 was signi cantly positively correlated with CFLAR, in SKCM, DNAJB9 was signi cantly positively correlated with CFLAR ( Figure S5B). In addition, our results also showed that CFLAR and PTK6 were signi cantly negatively correlated with tumor purity, PTK6 was signi cantly positively correlated with patient's age, CFLAR was signi cantly negatively correlated with patient's age, CFLAR and DNAJB9 were signi cantly positively correlated in purity and age ( Figure 5).

Tissue pro le and cell localization of CFLAR, DNAJB9 and PTK6
HPA database mRNA level tissue expression level showed that DNAJB9 was widely and highly expressed in all kinds of tissues, but especially low expression in skin. CFLAR was expressed in all kinds of tissues, and PTK6 was highly expressed in stomach, small intestine, colon, duodenum, esophagus and skin, with the highest expression in skin ( Figure 6A). In addition, CFLAR and DNAJB9 were localized in the cytoplasm of A-431 cell line, and PTK6 located in the nucleus ( Figure 6B). The protein level of medium was detected in skin and epithelial cells ( Figure 6C). In melanoma and skin cancer, PTK6 was mainly expressed in cytoplasm and cell membrane, CFLAR and DNAJB9 expressed in cytoplasm, cell membrane and nucleus.

Discussion
Autophagy plays a signi cant role in tumor promotion and suppression, depending on the cell/tissue types and the tumor stages, which hinders the clinical application of autophagy activators or inhibitors. Previous research have demonstrated that autophagy in involved in the regulation of melanocyte biology and skin pigmentation [ 17], ARGs have been reported to predict the survival of SKCM [21]. However, correlation between high risk signature and low signature is still unclear. Therefore, we aimed to identify ARGs associated with OS of patients with SKCM and explore tumor immunization and correlation of signi cant ARGs, thereby helping to develop individualized treatment option based on patient risk.
In present study, we identi ed 36 ARGs, which were dysregulated in patients with SKCM. After univariate and multivariate analysis, 8 ARGs were used to construct risk score signature using their expression levels weighted by corresponding correlation coe cient. The autophagy-related signature could divide patients into high-risk and low-risk groups based on the median risk score. Patients with high-risk score have signi cantly worse OS than patients in low-risk group. Furthermore, the predictive accuracy of autophagy-related signature.
In addition, we established an autophagy-based signature s based on multi-data sets, which showed favorable predictive ability. The results of the TCGA training data set (mainly composed of the US and European populations) are well validated in the validation set of the Asian population composition. Validation in multiple countries enhances the reliability of the results, indicating that ATG scores may be appropriate for patients of different ethnicities. Moreover, our nomogram combining risk score with conventional clinical parameters shown signi cantly improved performance. For example, the nomogram showed a higher predictive ability than the AJCC TNM staging system, consistent with previous studies. Nomogram is of great signi cance to clinicians and patients because that it can predict individual patient outcomes. The proposed autophagy-related signature included 8 ARGs. All genes in the signature were associated with OS of patients with SKCM.
In this present study, CFLAR, DNAJB9 and PTK6 were further con rmed to be the signi cant signature for prognostic factors in melanoma, CFLAR and DNAJB9 were low risk index, PTK6 was high index, which is consistent with previous molecular experiment. It reported that overexpression of PTK6 predicts poor prognosis in bladder cancer patients [30]. And PTK6/BRK is expressed in the normal mammary gland and activated at the plasma membrane in breast tumors [31]. Another research illustrated the clinical potential for PTK6 inhibition to improve treatment of patients with high-risk TNBC [32]. These studies support the results of PTK6, a signi cant high-risk marker identi ed in this study based on ARGs prognostic models. The CaMKII-mediated pathway for CFLAR (c-FLIP) upregulation protects melanoma cells from TRAILinduced apoptosis [33]. c-FLIP is a critical anti-cell death protein often overexpressed in tumors and hematological malignancies and its increased expression is often associated with a poor prognosis [34]. CFLAR (c-FLIP) also played important role on melanomas [35,36]. c-FLIP might play an important role in the obtaining of aggressive biological behaviour and be useful in predicting prognosis [35]. DNAJB9 inhibits p53-dependent oncogene-induced senescence (OIS) and induces neoplastic transformation under oncogenic RAS activation in mouse primary broblasts [37]. These studies support the low risk of CFLAR and DNAJB9 for prognostic factors in melanoma.
In tumor immunization, high and low risk factors seem to show opposite correlation characteristics. There is a signi cant positive correlation among the low risk genes, but there is no signi cant correlation between the high and low risk factors. CFLAR and DNAJB9 have stronger correlation with tumor immune cells. In addition, our results also showed that CFLAR and PTK6 were signi cantly negatively correlated with tumor purity; PTK6 was signi cantly positively correlated with patient's age, CFLAR was signi cantly negatively correlated with patient's age; CFLAR and DNAJB9 were signi cantly positively correlated in purity and age.
The results of tissue expression pro le showed that PTK6 was highly expressed in skin, while DNAJB9 and CFLAR were low expression in skin, which may determine that DNAJB9 and CFLAR are low risk index in prognostic model, and PTK6 is high risk index. The results of localization analysis showed that DNAJB9 and CFLAR had similar localization characteristics, both located in the cytoplasm and PTK6 located in the nucleus. Whether this is related to the high and low risk needs further experimental veri cation. In addition, the localization of PTK6 changed signi cantly in tumor patients, with nuclear transformation into cytoplasm and cell membrane [31,38]. In normal or non-tumorigenic tissues, PTK6 promotes cellular differentiation and apoptosis.

Conclusion
In conclusion, this study conducted a novel 8 autophagy-related signature to predict OS of patients with SKCM, and studied the correlation between validated ARGs based on a comprehensive analysis of ARGs expression pro les and corresponding, which may help in clinical decision-making for individualized therapy, and identi ed 3 ARGs signi cant to SKCM prognostic. Our ndings indicated that in tumor immunization analysis, high and low risk factors seem to show opposite correlation characteristics. There is a signi cant positive correlation between the low risk signature genes, but there is no signi cant correlation between the high and low risk genes.

Data collection
The gene expression and clinical information of 472 patients were obtained from TCGA (https://cancergenome.nih.gov/) and 81 paracancerous tissues were obtained from GETx https://commonfund.nih.gov/GTEx/ . Use this data in accordance with TCGA and GETx relevant policies to get the gene expression, GENCODE (https://www.gencodegenes.org) Version 28 annotated probes was used to obtain the corresponding gene names. All data were obtained from TCGA and CGGA, approval for our study by the ethics committee was not necessary.

Screening of Differential Genes
The limma package of R was used to analyze glioma and corresponding normal samples, and differentially expressed mRNA were screened out. Screening criteria was based on P<0.05 and Fold Change>1.
Identi cation of prognosis-related genes By using R/Bioconductor's "survival" package, one-to-one univariate survival analysis was performed for the differential gene expression, and the gene expression amount with P value less than 0.05 was selected for further construction of the model. Selection of best genes for modeling.
The prognostic risk score model was established with the following formula: risk score = expression level of Gene1 × β1 + expression level of Gene2 × β2 +…+ expression level of Genen × βn, βis a regression coe cient model calculated by multivariate Cox regression. Subsequently, a prognostic risk score was generated for each patient. All patients were divided into high risk (high risk score) and low risk (low risk score) groups according to the median of their risk score. Then, K-M survival curves were constructed to estimate the prognosis of patients with high or low risk scores, and survival differences between high and low risk groups were assessed by two-sided log-rank test. and the grouping boundary value was calculated by the "survival" package and "survminer" package of R/Bioconductor, and the resulting two groups were separately Kaplan-Meier curves were plotted. At the same time, in order to verify the validity and repeatability of the model, patients were randomly divided into training group and validation group at a ratio of 5:5 using R. In the training group, patients were divided into high-risk group and low-risk group according to the risk score. The grouping boundary value was calculated from the "survival" package and "survminer" package of R/Bioconductor, and Kaplan-Meier curves were drawn for the two groups.

Development of the gene prognostic nomogram
Using 472 training samples from TCGA, we generated a prognostic nomogram using the R package "rms" based on the expression level of genes selected in the previous step. In the software package, the "cph" function is used to construct the COX model. Based on this model, a prognostic nomogram was generated using the "nomogram" function. Line length corresponding to each variable in the prognostic nomogram re ects the contribution of predictive factors to patient prognosis.

Functional analysis
The Bohao Online Enrichment Tool (http://enrich.shbio.com/) was used to perform functional enrichment of differentially expressed ARGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used to assess relevant functional categories. GO and KEGG enrichment pathways with P and Q values less than 0.05 are considered important categories.
Tumor immunization analysis TIMER was applied to download score data of six kinds of immune in ltrating cells. we used R software package estimate to analyze the immune score and stromal score of SKCM samples to explore the relationship between CFLAR, DNAJB9, PTK6 and tumor development in tumor microenvironment [28].

Declarations
Ethics approval and consent to participate  Validation of signi cant autophagy-related genes index (A) Survival curve of 8-ARGs based on TCGA and GETx database, (B) Overall survival of 8-ARGs based on GEPIA tools. Cutoff-High(%): Samples with expression level higher than this threshold are considered as the high-expression cohort. Cutoff-Low(%): Samples with expression level lower than this threshold are considered the low-expression cohort., (C) This feature generates expression violin plots of signi cant ARGs based on patient pathological stage. The method for differential gene expression analysis is one-way ANOVA, using pathological stage as variable for calculating differential expression.

Figure 4
Correlation between high risk ARGs and low risk ARGs in patients with SKCM TIMER online tool was applied to study the correlation of high risk ARGs and low risk ARGs in SKCM adjusted by tumor purity or age.