Identification 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 up-regulated 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.
To analyze ARGs involved in skin cutaneous melanoma progression, ARGs were significantly associated with prognosis. The forest map of the hazard ratio indicates that most of these genes are protective factors, 14 high risk genes (ULK, SPNS1, TP63, HSP90AB1, BAK1, GAPDH, ATG9B, EEF2K, BNIP3, PTK6, BIRC5, DAPK1, PARP1 and EGFR) and 22 low risk gens (PRKCQ, CASP3, CALCOCO2, DAPK2, ATG16L2, EDEM1, SERPINA1, FAS, KIF5B, VAMP7, APOL1, CCR2, ITGA6, DNAJB9, GAA, CFLAR and CXCR4) (Figure 1A). Both GO (GO:0006914) and KEGG (hsa04140) analysis showed that these genes are closely related to autophagy process (Figure 1B and 1C). 57% high risk genes were down-regulated in melanoma, 72% low risk genes were up-regulated in melanoma. Given the important clinical implication of these ARGs, we examined that mRNA genetic alterations of these genes and found that amplification is the most common type of mutations (Figure 1D). A total of 13 genes have a mutation rate ≥ 5%, of which TP63 is the most frequently mutated gene (17%).
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 significantly 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 profile 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 significantly 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 significantly 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-1.263, P<0.01) (Figure S2B). In the validation set, the distribution of survival status and risk scores of patients have a similar trend to that in the training set. Consistent with the results in the training set, survival analysis showed a significantly lower OS (P<0.001). The overall analysis showed the area under the curve of the corresponding receiver operating characteristic (ROC) curve for 1 year, 3 years, and 5 years of survival is 0.725, 0.696 and 0.740, respectively. The training set showed the ROC for 1 year, 3 years and 5 years of survival is 0.732, 0.810 and 0.854. The test set showed the ROC for 1 year, 3 years and 5 years of survival is 0.720, 0.583 and 0.634. This indicated that the prognostic index based on ARGs has a certain potential in survival prediction (Figure S2C.
Validation of significant autophagy-related genes index for SKCM patients
To Construct validation of nomogram, we subsequently analyzed the relationship between ARGs prognostic index and clinical features. Significant increase in risk score were in metastasis (Figure S3A), late clinical stage (Figure S3B), in larger tumor size (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 significant 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 significantly 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 profile 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 significant ARGs in melanoma and immunity, we downloaded the score data of six kinds of immune infiltrating 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 significant 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 significantly negative correlated to B_ Cell (Figure 4A). In addition, the results showed that only CFLAR was significantly positively correlated with the overall microenvironment score, CFLAR and DNAJB9 were significantly positively correlated with immune score, CFLAR, DNAJB9 and PTK6 were significantly 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 significant correlation (Figure 4C). SCNA module of target genes were analyzed by TIMER, which provides the comparison of tumor infiltration 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 infiltration of B cell, CD4+T cell, macrophage, and dendritic cell, and the copy number of arm level gain was related to the infiltration of CD4 + T cell. The copy number of DNAJB9 arm level deletion was correlated with B cell, CD4+T cell and dendritic cell infiltration. Copy number of PTK6 high amplification was significantly 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 infiltration, copy number of arm level deletion was correlated with B cell, CD8+T cell, CD4+T cell, and dendritic cell infiltration (Figure S4A). Immune checkpoint genes also play an important role in tumor immunotherapy. Our result found that PTK6 only had a significant correlation with VTCN1, CFLAR was significantly correlated with most immune checkpoint genes, but not with TNFRSF14, CD276, VTCN1, HHLA2, CD70, TNFSF9, CD44. DNAJB9 significantly 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 significantly positively correlated with CFLAR, in SKCM, DNAJB9 was significantly positively correlated with CFLAR (Figure S5B). In addition, our results also showed that CFLAR and PTK6 were significantly negatively correlated with tumor purity, PTK6 was significantly positively correlated with patient's age, CFLAR was significantly negatively correlated with patient's age, CFLAR and DNAJB9 were significantly positively correlated in purity and age (Figure 5).
Tissue profile 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.