3.1 Acquisition of genes related to inflammatory response
The study’s flow chart is depicted in Figure 1. The study population included 472 patients with SKCM from the UCSC cohort and 323 patients with normal skin tissue expression data. The clinical characteristics of these patients are summarized in Table 1(including age, gender, stage, and TNM stage). The screening of prognostic genes associated with inflammation was conducted with data on 200 prognostic genes associated with inflammation received from the GSEA database (additional file 1).
3.2 Construction and validation of a predictive model related to inflammatory response
We removed samples with identical ID values and missing clinical data, leaving 454 samples with survival days and survival status. These samples were randomly divided into a training group (n=228) (additional file 2) and a validation group (n=226) (additional file 3). We initially screened using one-way Cox analysis and identified 15 inflammatory genes linked with OS (P<0.01) (Fig 2A) (additional file 4) and demonstrated their correlation (Fig 2B). We then performed LASSO regression analysis on these genes and screened inflammatory genes with more than 900 replicates in 1000 replacement samples. We then constructed a model containing five genes to predict the prognosis of SKCM patients based on their λ values (Fig 3). The model consists of five key genes, including C3AR1, CXCL10, EIF2AK2, EMP3, ICAM1. The hazard ratio(HR), 95% confidence interval(95%CI) and P value of five hub gene could be seen at Table 2. The risk score is calculated as follows: (-0.0019570292110695 4)×Expression(C3AR1)+(-0.0717393323519231)×Expression (CXCL10)+(-0.178228 168474953)×Expression (EIF2AK2)+(0.082105091032782 3)×Expression (EMP3)+ (-0.00085177603520393 1)×Expression (ICAM1).
We split all patients in the training and validation groups into high- and low-risk groups based on the training group’s median risk score. In the training group, we discovered that patients’ chance of death increased and their survival time reduced (Fig 4A). Scatter plot analysis revealed that patients classified as high-risk had a greater likelihood of dying sooner than patients classified as low-risk (Fig 4B). Kaplan-Meier curves revealed that OS was substantially worse in high-risk individuals than in low-risk patients (P<0.01) (Fig 4C). The risk heat map depicted the expression of many inflammatory genes associated with prognosis in high- and low-risk groups (Fig 4D). We next re-validated our results in the validation group, where we noticed substantial variations in OS between the high-risk and low-risk groups (Fig 4E-4G) and inflammatory prognostic gene expression in the heat map (Fig 4H) proving the model’s accuracy once again. We then used principal component analysis (PCA) (Fig 5A-B) and t-SNE (Fig 5C-D) analysis to demonstrate the discrete distribution of patients into different risk groups. Because patients in the high-risk group may die earlier and have shorter survival times than those in the low-risk group, we can intuitively assume that the model constructed using inflammatory prognosis-related genes can better differentiate the prognosis of patients with SKCM.
3.3 Independent prognostic analysis of OS
We utilized univariate and multivariate Cox regression analyses to determine if clinical parameters (age, gender) and risk score were significant independent predictors of OS. We showed that risk score, gender, and age were independent predictive predictors of OS in univariate and multivariate Cox regression analyses in the training group (Fig 6A-B), validation group (Fig 6C-D), and all samples (Fig 6E-F). Additionally, the ROC areas at 1, 2, and 3 years were 0.754, 0.660, 0.654 (Fig 7A) and 0.670, 0.709, 0.682 (Fig 7B), respectively, for the training and validation groups. Our findings indicated that the risk score model was considerably better in predicting OS in patients with SKCM than other clinical factors, such as age and gender (Fig 7C-D).
3.4 Immune status and tumor microenvironment analysis
To further synthesize the immune cell infiltration in SKCM patients and to explore the differences in immune status between high- and low-risk groups, we applied TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC in the training group, and the results were presented in a heat map (Fig 8). We then used ssGSEA to quantify the enrichment scores of different immune cell subpopulations, related functions and pathways to explore the differences in immune cell enrichment scores between the high-risk and low-risk groups. In the training set we found that the antigen presentation process contained significantly higher levels of aDCs, pDCs, APC-co inhibition, APC-co stimulation, HLA, and MHC class I in the low-risk group compared to the high-risk group. In addition, we found that CD8+T cells, T helper cells, Tfh cells, Th1 cells, Th2 cells, TIL cells, Treg cells, T cell co-inhibition and T cell co-stimulation associated with T cell regulation were all significantly higher in the low-risk group than in the high-risk group. The T helper cells, Tfh cells, Th1 cells, Th2 cells, TIL cells, Treg cells, T cell co-inhibition and T cell co-stimulation were all significantly higher than those in the high-risk group, suggesting that there were differences in T cell regulation between the low-risk and high-risk groups. In addition, we found that the score of CCR and Checkpoint, cytolytic activity, B cells, NK cells, neutrophils, inflammation promoting, para-inflammation, Type I IFN response and other biological processes and cellular contents were significantly higher in the low-risk group. response and other biological processes and cellular content were significantly higher than those in the high-risk group (Fig 9A,9C). We found similar results in the validation set as in the training set (Fig 9B,9D).
To examine the association between risk scores and immunological components in further detail, we conducted a correlation analysis between risk scores and immune infiltration. Six distinct immune infiltration types were identified in human tumors, including C1 (wound healing), C2 (IFN-γ dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-β dominant), although the C5 and C6 immune subtypes were excluded due to their absence in SKCM patients. We evaluated immune infiltration of SKCM in the training set and connected it with the risk score, finding that a high risk score was strongly associated with C3, whereas a low risk score was significantly associated with C2 (Fig 9E).
Tumor stemness can be determined using RNA stemness based on mRNA expression (RNAss) or DNA stemness based on DNA methylation pattern (DNAss). Additionally, we evaluated tumor immunological microenvironment using the immuneScore and stromalScore. Correlation analysis examined the relationship between risk score and tumor stemness and immune microenvironment. The results indicated that the risk score was not significantly correlated with DNAss or RNAss, but was significantly negatively correlated with immuneScore and stromalScore (P<0.05) (Fig 9F).
3.5 Biological function and pathway analysis
We used GSEA to conduct enrichment analysis on KEGG pathways in the high- and low-risk group (Fig 10A). The results analysis revealed that the five pathways with the highest content in the high-risk group were Alzheimer’s disease, aminoacyl tRNA biosynthesis, chemokine signaling pathway, cytokine receptor interaction, and Huntington’s disease (P<0.05, FDR<25%) (Fig 10B). The JAK-STAT signaling route, oxidative phosphorylation, RNA polymerase, toll-like receptor signaling pathway, and viral myocarditis were the five pathways with the highest levels in the low-risk group (P<0.05, FDR<25%) (Fig 10C).
3.6 Prognostic gene expression and sensitivity of melanoma cells to chemotherapy
To make this study more therapeutically relevant, we examined the expression of prognostic genes in NCI-60 cell lines and the association between prognostic gene expression levels and drug sensitivity. The results indicated that all prognostic genes were significantly associated with drug susceptibility to certain chemotherapy agents (P<0.01). The figure depicts the findings of the drug sensitivity analysis with the highest correlation. CXCL10 expression was associated with increased cancer cell resistance to LDK-378, brigatinib, alectinib, and PF-06463922, among others. C3AR1 expression was associated with higher resistance of cancer cells to denileukin diftitox ontak, isotretinoin, carmustine, estramustine, fluphenazine, nelfinavir, megestrol acetate, alectinib, cyclophosphamide, lomustine, and dromostanolone propiona, among others (Fig 11).
3.7 Validation of key genes in melanoma and paraneoplastic tissues
To confirm the changes in expression of 5 important genes (C3AR1, CXCL10, EIF2AK2, EMP3, ICAM1) between SKCM and para-cancerous normal tissue, we first investigated their mRNA expression using quantitative real-time PCR(qRT-PCR). qRT-PCR analysis revealed that prognostic genes were expressed at a higher level in melanoma than in normal tissues adjacent to malignancy (P<0.01) (Fig 12A-E). We then validated the expression of key genes in SKCM and normal skin using the Human Protein Atlas (HPA) database; however, because the HPA database lacked immunohistochemical images of normal and tumor tissues for CXCL10 and EMP3, we validated only the expression of the remaining three key genes, which were all significantly different (Fig 13A -C).