Development and validation of an oxidative stress—associated prognostic risk model for melanoma

Processing of raw data

RNA sequencing samples from 472 individuals were obtained from the TCGA database, which comprised of 471 SKCM samples and one healthy skin tissue sample (https://portal.gdc.cancer.gov/). In addition, to increase the number of healthy samples, we collected 812 RNA sequencing data from healthy skin tissues obtained from the GTEx database (https://gtexportal.org/home/datasets) (Human Genomics, 2015; Gentles et al., 2015). A total of 1399 OS-related genes with relevance score ≥ 7 were collected from the Gene Cards database (https://www.genecards.org). Under —log2 fold change—≥ 2 and a standard false discovery rate (FDR) <0.05, an R package was used to detect genes differentially expressed in SKCM and healthy skin samples (Li et al., 2020). Meanwhile, an average count cutoff of 1 was used to eliminate genes. After univariate and multivariate Cox analyses, 16 OS-related genes associated with SKCM prognosis were obtained for the construction of the risk model. We used the National Center for Biotechnology Information-Gene Expression Omnibus database to download the GSE65904 dataset, a cohort of 214 SKCM patients for external verification. The basic characteristics of these SKCM samples were all displayed in Table 1.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and Gene Ontology (GO)

KEGG enrichment analysis and GO of OS-associated differentially expressed genes (DEGs) were performed to assess their biological functions. The Database for Annotation, Visualization, and Integrated Discovery 6.8 (Huang, Sherman & Lempicki, 2009) was used to perform the analyses.

Establishing protein-protein interaction (PPI) network and screening for important modules

PPI information from OS-associated DEGs was acquired from the STRING platform (http://www.string-db.org/) (Szklarczyk et al., 2019). The PPI network was developed using Cytoscape 3.7.0. The virtual modules and hub genes in the PPI network with a Molecular Complex Detection (MCODE) score and node count >5 (p < 0.05) were selected using the MCODE plug-in (Bader & Hogue, 2003).

Construction of a prognostic risk model

The univariate Cox regression of hub genes found in the PPI network was performed using the ‘survival’ R package. Subsequently, multivariate Cox regression was used to further analyze the above OS-related genes and select genes to build the prognostic risk model. We calculated the risk value of patients with SKCM based on the expression and coefficient values of 16 OS genes and classified them into high- and low-risk queues through the median risk score. The risk score was determined according to the equation below: ${\displaystyle \text{Risk score}=\Sigma \mathrm{expgenei}\times \beta i,}$ where expgenei is the expression level of 16 OS-associated genes and β is the coefficient value for the gene. Receiver operating characteristic (ROC) curves were generated using the ‘timeROC’ and ‘survivalROC’ packages for R and were applied to assess the accuracy of our risk model to predict the overall survival rate of SKCM patients. Separate nomograms were constructed based on clinical characteristics and the 16 OS genes. The calibration chart was used to detect the predictive power of the above nomograms (one based on clinical characteristics and one based on the 16 OS genes) for the overall survival time of SKCM patients. To validate the prognostic performance of the constructed risk model, the analyses described above were conducted using data from the GSE65904 cohort.

Validation of expression levels and prognostic values of hub genes

After clarifying the translational expression level of hub genes in the TCGA cohort, we verified the differential expression of 16 OS genes between SKCM and normal skin tissues using data from the Human Protein Atlas (HPA) (Thul et al., 2017). In the TCGA cohort, a survival analysis of 16 OS genes was performed using the Kaplan–Meier (KM) approach to ascertain whether they were correlated with the prognosis of SKCM patients.

Identification of OS-associated DEGs

Figure 1A shows the workflow of the study. In total, 1399 OS genes were selected to study their differential expression in SKCM and healthy tissues. Of these, 156 were identified as OS-associated DEGs (63 upregulated and 93 downregulated genes) in SKCM (Fig. 1B).

Functional enrichment analysis of OS-associated DEGs

We used the KEGG pathway analysis to analyze the DEGs and found that upregulated genes were mostly correlated with cytokin-cytokin receptor interaction and hunman T-cell leukemia virus 1 infection (Fig. 2A), while downregulated genes were predominantly linked with fluid shear stress and atherosclerosis (Fig. 2B). GO analysis was also performed to further explore the DEGs. As a result, with regard to biological processes, upregulated DEGs were significantly augmented in leukocyte migration and leukocyte chemotaxis (Fig. 3A), whereas DEGs that were downregulated were mainly augmented in response to OS and cellulare response to OS (Fig. 3B). With regard to cell location, upregulated genes were mainly augmented on the external side of plasma membrane and secretory granule membrance (Fig. 3A). Downregulated genes were enriched in the vesicle lumen and cytoplasmic vesicle lumen (Fig. 3B). Finally, it was evident that upregulated OS genes were concentrated in cytokine activity and cytokine receptor binding (Fig. 3A), while downregulated OS genes were mostly implicated in antioxidant activity and heme binding (Fig. 3B).

Creation of a PPI network for OS-associated DEGs and screening of key modules

To further explore the inner relationship of OS-associated DEGs, a PPI network was established with 144 nodes and 834 edges (Fig. 4A). The most meaningful module with 19 nodes and 158 edges was subsequently identified (Fig. 4B). OS-related genes within the key module were primarily involved in chemokine-mediated signaling transduction, leukocyte migration, response to chemokine, and chemokine signaling pathway.

Screening of hub genes and construction of a prognostic risk model

A total of 144 OS-associated DEGs were identified from the PPI network. After univariate Cox regressions, 61 OS genes were identified as genes of prognostic value in SKCM patients (Fig. 5A). The multivariate Cox regression model helped to select 16 hub genes (CDK2, CCR5, NDUFA9, NDUFA13, HLA.DRB1, CXCR3, FOXM1, CCL4, ISG15, FCGR2A, FCGR3A, PIK3R2, SLPI, SELL, PSEN2, and GJA1) that were used to calculate the prognostic risk model (Fig. 5B, Table 2).

Validation of the prognostic value of the risk model

The median risk score of the prognostic risk model was used to separate SKCM patients in the TCGA and GSE65904 cohorts into high- and low-risk queues (Figs. 6A, 6B). The survival time of SKCM patients in the high-risk group was significantly lower when compared to the low-risk group (Figs. 6C, 6D). A ROC curve was constructed to validate the accuracy of our risk model. It indicated that our risk model had a moderate predictive effect using data from the TCGA cohort (area under the ROC curve [AUC] of 1-year survival = 0.742; 3-year survival = 0.723; and 5-year survival = 0.764) (Fig. 6E). Similarly, both the prognostic effect and accuracy were validated in the GSE65904 cohort, in which the AUC of 3-year survival was 0.705 (Fig. 6F). The univariate and multivariate Cox regressions of different clinical characteristics of SKCM patients in the TCGA cohort showed that the gene-based risk score was a robust prognostic parameter for SKCM patients (Figs. 7A, 7B). When compared with other clinical features in the TCGA and GSE65904 cohort, our prognostic risk model showed a better prognostic performance in all AUCs (Figs. 7C–7H), revealing that our gene-based prognostic risk model displayed moderate specificity and sensitivity with regard to predicting SKCM prognosis. The correlation analysis between clinical parameters and risk score indicated that SKCM patients in the T3 and T4 stages or those with primary cancer had a higher risk score. (Figs. 7I, 7J). A heatmap was drawn to show the correlation between the levels of 16 OS genes in the TCGA and GSE65904 cohorts against different clinical parameters, including high- and low-risk groups, TNM stage, age, and gender (Fig. 7K, 7L).

Furthermore, in the TCGA and GSE65904 cohorts, the nomograms of our risk score and different clinical parameters were constructed to predict the overall prognosis of SKCM patients (Fig. 8A, 7B). In parallel, the calibration diagram evidenced that the above nomograms had a good predictive effect on the clinical outcome of SKCM patients (Figs. 8C–8F).

Validating the prognostic value and expression levels of hub genes

In SKCM samples, the expression levels of CDK2, CCR5, HLA-DRB1, CXCR3, FOXM1, CCL4, ISG15, FCGR2A, FCGR3A, SELL, and PSEN2 were considerably increased, while the levels of NDUFA9, NDUFA13, PIK3R2, SLPI, and GJA1 were markedly decreased when compared with that in healthy samples (File S1A). These results were validated by immunohistochemistry analysis from the HPA database (File S1B).

To visualize the interactions between hub DEGs, we constructed a PPI network using the STRING database online tool (Fig. 9A). The genes FCGR2A and CCR5 showed the highest interacting degrees among the hub genes (Fig. 9B). The KM method showed that a high expression level of HLA-DRB1, CXCR3, CCL4, ISG15, FCGR2A, FCGR3A, SELL, and CCR5 were correlated with a significant increase in the overall survival rate in SKCM, whereas a high expression level of PSEN2, CDK2, FOXM1, GJA1, NDUFA9, NDUFA13, PIK3R2, and SLPI were correlated with a significant decrease of the overall survival rate (File S2). Moreover, as shown in File S3, the genes SELL, PSNE2, FOXM1, CDK2, and HDUFA9 were all significantly related to with SKCM ages, while genes HDUFA13 and CCL4 were significantly connected the ganders of SKCM patients. For the TCGA and GSE65904 cohorts, we also constructed the nomograms related to the 16 OS genes to predict the 1-, 3-, and 5-year survival probability of SKCM patients (File S4). The calibration of the nomograms associated with the 16 OS genes presented good consistency between the predicted and observed outcomes.