Identifying novel circadian rhythm biomarkers for diagnosis and prognosis of melanoma by an integrated bioinformatics and machine learning approach

Data download and preprocessing

Firstly, 4 melanoma-related transcriptomic datasets (GSE15605, GSE114445, GSE46517 and GSE65904 datasets) were downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/), and 1 melanoma-related gene expression profile data was obtained from TCGA database (https://tcga-data.nci.nih.gov/tcga/). A total of 1471 circadian rhythm related genes (CRGs) were obtained from MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/) and Genecards database (https://www.genecards.org/) (Supplementary Table 1). Then, preprocessing was performed on each expression profile data. The R software package “limma” was used for background correction, median normalization and gene symbol conversion of each dataset.

GSE15605 and GSE114445 datasets were based on GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array). GSE15605 dataset contained 16 normal skin samples and 46 primary melanoma samples, and GSE114445 dataset contained 6 normal skin samples and 15 primary melanoma samples. GSE46517 dataset was based on GPL96 platform (Affymetrix Human Genome U133A Array), containing 8 normal skin samples and 31 primary melanoma samples. GSE65904 data was based on GPL10558 platform (Illumina HumanHT-12 V4.0 expression beadchip), containing tumor tissue samples from 214 melanoma patients. TCGA-melanoma dataset included tumor tissue samples from 458 melanoma patients.

In this study, GSE15605 and GSE114445 datasets were combined as a new expression profile as the training set, which was used for the selection of key CRGs in melanoma. Separate GSE15605, GSE114445 and GSE46517 datasets were used as validation sets to validate the expression of key CRGs and evaluate the diagnostic ability of key CRGs. The TCGA-melanoma cohort was used as the training set for prognostic evaluation and prognostic model construction of key CRGs. GSE65904 dataset was used as the validation set for verification of the prognostic model.

Screening of melanoma-related circadian rhythm genes (CRGs)

Dataset merging

Firstly, GSE15605 and GSE114445 datasets were merged as the training set for screening melanoma-related key CRGs. To eliminate batch effects between different datasets, batch correction was performed on the merged dataset using the “sva” package in R software. Then the box plots, density plots of gene expression and Uniform Manifold Approximation and Projection (UMAP) analysis were performed to validate the effect of batch correction.

Differential expression analysis

Differential expression analysis was performed using R software package “limma”. Differential expression genes (DEGs) between melanoma and normal skin samples in the training set were screened with the threshold of “|logFC|>1.2, P<0.05”.

Weighted gene co-expression network analysis (WGCNA)

The R package “WGCNA” was used to perform weighted gene co-expression network analysis (WGCNA) on the expression matrix of the above DEGs. Specifically, we first constructed an adjacency matrix to describe the association strength between nodes. Then, the optimal soft thresholding power β was chosen to transform the adjacency matrix into topological overlap matrix (TOM), making the constructed network follow scale-free topology and be closer to real biological networks. Next, module partitioning analysis was performed using hierarchical clustering and dynamic tree cut algorithm to determine gene co-expression modules. The module eigengenes (MEs) of each gene module were then calculated and the connectivity between different modules was analyzed. The correlation between each module and clinical trait (melanoma and normal skin) was then calculated. Finally, gene significance (GS) and module membership (MM) were calculated, and hub genes closely related to melanoma pathogenesis were screened in the modules with GS>0.6, MM>0.6 and P<0.05.

Screening of melanoma-related CRGs

The hub genes obtained from the WGCNA network were compared with the CRGs downloaded from the MSigDB database to identify the common genes, which were the melanoma-related CRGs that affect the occurrence and development of melanoma.

Functional enrichment analysis and protein-protein interaction (PPI) analysis

To explore the potential biological functions and signaling pathways of circadian rhythm genes (CRGs) in melanoma pathogenesis, we performed functional enrichment analysis on melanoma-related CRGs using DAVID database. Specifically, Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis were performed on the above melanoma-related CRGs. For GO functional annotation, we selected “Gene Ontology” as the analysis type in the DAVID database and chose “Biological Process”, “Molecular Function”, and “Cellular Component” as the annotation types. We selected the species as human and submitted the list of CRGs. For KEGG functional enrichment analysis, we used the DAVID database and selected “KEGG Pathway” as the analysis type. Finally, we visualized the analysis results.

To further evaluate the role of circadian rhythm mechanisms in melanoma pathogenesis, we performed Gene Set Enrichment Analysis (GSEA), a gene set-based enrichment analysis method that identifies significant differences in biological functions and signaling pathways between two biological conditions. Specifically, we obtained Gene Ontology (GO) biological processes gene sets, including ENTRAINMENT OF CIRCADIAN CLOCK, CIRCADIAN REGULATION OF GENE EXPRESSION, CIRCADIAN RHYTHM, and REGULATION OF CIRCADIAN RHYTHM, as well as the KEGG gene set CIRCADIAN RHYTHM, from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb) as reference gene sets. We then used GSEA software (version 4.0.1) to determine the expression differences of these circadian rhythm-related biological functions and signaling pathways between melanoma and normal skin tissues.

In addition, the protein-protein interaction (PPI) network between these melanoma-related CRGs was further constructed. First, these melanoma-related CRGs were uploaded to STRING database (https://string-db.org/) for protein-protein interaction (PPI) analysis. Then Cytoscape software (version 3.7.2) was used to construct the PPI network and analyze the topological parameters of nodes in the network.

Identification of key CRGs associated with melanoma

We utilized 3 machine learning algorithms including support vector machine-recursive feature elimination (SVM-RFE) analysis, random forest analysis and LASSO regression analysis to screen key genes from the above melanoma-related CRGs.

SVM-RFE is an embedded feature selection approach that ranks features based on their importance [11]. It repetitively trains SVM classifier models and removes the feature with smallest ranking criterion each time until all features are ranked. Here, we utilized SVM-RFE algorithm in R package “e1071” to rank the melanoma-related CRGs by their importance. The CRGs were recursively eliminated until the optimal SVM model with high classification accuracy was obtained. The remaining CRGs were identified as key features.

Random forest is an ensemble supervised learning technique constructed by multiple decision trees [12]. It can estimate the importance of features by permutation and assess prediction error using out-of-bag samples. Here, we used the random forest algorithm in R package “randomForest” to select and rank CRGs by their importance. The top 20 ranked CRGs by importance were selected as key genes under random forest analysis.

LASSO regression is a regularization technique that performs feature selection and regularization to enhance prediction accuracy and interpretability. It penalizes the absolute size of regression coefficients, effectively reducing model complexity and avoiding overfitting [13]. Here, we constructed a LASSO regression model based on melanoma tissue and normal skin tissue samples using R package “rms”. By tuning the regularization parameter λ, LASSO regression selects important features and shrinks the coefficients of unimportant features to zero. In this study, we identified the melanoma-related key CRGs by selecting the features with non-zero coefficients in the LASSO model when λ was set as 0.04.

The commonly selected feature genes by the 3 machine learning methods were identified as key CRGs associated with melanoma.

Expression and ROC analysis of key CRGs

Based on the training set (merged GSE15605 and GSE114445 datasets) and validation sets (GSE15605, GSE114445 and GSE46517 datasets), Wilcoxon rank-sum test was performed to investigate the expression of key CRGs in melanoma tissues and normal skin tissues, respectively.

Then, to further validate the accuracy of key CRGs selection and evaluate the diagnostic value of key CRGs as melanoma biomarkers, receiver operating characteristic (ROC) analysis was performed on the training and validation sets using the “pROC” package in R software. The diagnostic power was reflected by the area under ROC curve (AUC) value. AUC>0.7 indicated that key CRGs had good ability to distinguish melanoma from normal skin tissues.

Prognostic analysis of key CRGs

Firstly, we performed Kaplan-Meier survival curve analysis for each key CRG using the ‘survival’ package in R software based on the TCGA database which integrated transcriptomic data and clinical information of melanoma patients, in order to assess their impacts on the prognosis of melanoma patients. Then, multivariate Cox regression analysis with stepwise variable selection was utilized to construct a prognostic model with key CRGs significantly associated with prognosis of melanoma patients, and determine the relative coefficient of each feature gene. At each step, the most significant CRG was added into the model. CRGs that became non-significant were removed from the final prognostic model. For each patient, the CRGs risk score was calculated as follows:

$$\text{Risk}\text{Score}={\displaystyle \sum _{i=0}^{n}Coefficien{t}_{(mRNAi)}}\times {\text{Expression}}_{(mRNAi)}$$

According to the median value of risk scores, melanoma patients in the training set (TCGA cohort) and validation set (GSE65904 dataset) were classified into high- and low-risk groups. Kaplan-Meier survival analysis was performed to evaluate the survival difference between high- and low-risk groups. ROC analysis at 3, 5 and 10 years was performed using the R package “survivalROC” to evaluate the sensitivity and specificity of the prognostic model.

Then, to further predict the prognosis of melanoma patients, we constructed a nomogram. The nomogram is a statistical prediction model that generates a visualization tool for estimating individualized probabilities of a clinical outcome, such as patient survival, based on a combination of important variables. Here, we analyzed the relationship between different clinical pathological parameters including stage, Tumor (T), Node (N), Metastasis (M) classification and risk score using Spearman’s rank correlation test. Benjamini-Hochberg procedure was applied for multiple testing correction. Subsequently, integrating clinical features (age, sex, stage, T/N/M classifications) and risk score, univariate and multivariate Cox regression analyses were performed. Using the rms package in R software, all independent prognostic parameters were utilized to construct a nomogram to predict 3-, 5- and 10-year overall survival of melanoma patients. The predictive ability of the nomogram was validated by calibration analysis and ROC analysis.

Screening of melanoma-related CRGs

Dataset merging

We first merged GSE15605 and GSE114445 datasets as the training set for screening melanoma-related key CRGs. After batch correction, we plotted the gene expression boxplots and density plots, and performed Uniform Manifold Approximation and Projection (UMAP) analysis to validate its effects. The gene expression boxplots and density plots showed that the data distributions across different datasets became consistent after removing batch effects (Figures 1A, 2B). The UMAP analysis results showed that each sample in GSE15605 and GSE114445 datasets was randomly and discretely distributed (Figure 1C). The results indicated that batch correction could effectively eliminate batch effects between different datasets, making the merged dataset more consistent and reliable.

Screening of melanoma-related CRGs

Through differential expression analysis, 8725 DEGs were identified between melanoma tissues and normal skin tissues (Figure 2A, 2B). These DEGs were further analyzed by WGCNA to screen genes associated with melanoma. In WGCNA, the optimal soft thresholding power β was chosen as 6 (Figure 2C, 2D). Then modules were cut using hierarchical clustering and DynamicTreeCut algorithm. With the maximum module distance set as 0.25, highly similar modules were merged, and eventually 15 independent gene modules were generated, namely greenyellow, saddlebrown, yellow, midnightblue, blue, darkorange, lightyellow, darkred, paleturquoise, darkgreen, cyan, green, salmon, white and grey (Figure 2E). We then analyzed the independence of these gene co-expression modules. As shown in Figure 2F, the distances between modules were all greater than 0.25, indicating these gene modules had good independence. Next, we calculated the correlation between all modules and clinical trait (melanoma and normal skin). The results showed that except grey module, all other modules were significantly correlated with melanoma (Figure 2G). Therefore, these melanoma-related modules were selected for further screening genes closely related to melanoma pathogenesis. With the threshold of GS>0.6, MM>0.6 and P<0.05, 1643 hub genes were identified from the co-expression network of these gene modules (Supplementary Table 2). These hub genes were considered to be closely related to melanoma pathogenesis.

Finally, by comparing CRGs and these hub genes, 125 common genes were identified and considered as CRGs that affect melanoma pathogenesis and development (Figure 2H).

Functional enrichment analysis and protein-protein interaction (PPI) analysis

Here, the biological functions and signaling pathways related to circadian rhythm mechanism in melanoma were identified. GO functional enrichment analysis of the 125 melanoma-associated circadian rhythm genes (CRGs) revealed their involvement in molecular functions (MFs) including identical protein binding, steroid binding, transcription factor activity, sequence-specific DNA binding, and transcription corepressor binding (Figure 3A). These CRGs were enriched in biological processes (BPs) such as rhythmic process, regulation of circadian rhythm, circadian regulation of gene expression, and response to xenobiotic stimulus (Figure 3B, 3E), and were associated with various cellular components (CCs) including chromatin, extracellular exosome, cytosol, cytoplasm, and nucleoplasm (Figure 3C). KEGG analysis showed that besides directly regulating the circadian rhythm pathway, these melanoma-related CRGs may also indirectly participate in the circadian mechanisms of melanoma by modulating other signaling pathways such as cellular senescence, bile secretion, apelin signaling pathway, insulin signaling pathway, and NF-kappa B signaling pathway (Figure 3D, 3F).

GESA analysis further highlighted the significant roles of circadian-associated biological processes and pathways in melanoma. As shown in Figure 4G, 4H, four circadian-related BPs (GOBP ENTRAINMENT OF CIRCADIAN CLOCK, GOBP REGULATION OF CIRCADIAN RHYTHM, GOBP CIRCADIAN RHYTHM, GOBP CIRCADIAN REGULATION OF GENE EXPRESSION) and one KEGG pathway (KEGG CIRCADIAN RHYTHM) were significantly upregulated in the tumor tissues. In addition, we evaluated the interactions between these melanoma-related CRGs by establishing a PPI network. The results showed these melanoma-related CRGs were highly interconnected. Genes including CDK2, AR, CDK4, CEBPA, and CXCL8 occupied central positions in the network control, implying their potential regulatory interactions with numerous genes in the PPI network (Figure 3I).

Identification of key CRGs associated with melanoma

Three machine learning algorithms including LASSO regression, SVM-RFE, and random forest were utilized to select key CRGs from the 125 melanoma-related CRGs. Firstly, SVM-RFE analysis identified 20 feature CRGs whose classification model achieved 92.9% accuracy in evaluating melanoma samples (Figure 4A). Then, random forest algorithm was used to select and rank CRGs by importance. The top 20 CRGs by importance ranking were selected as key genes under random forest analysis (Figure 4B). Further, a LASSO regression model was constructed based on melanoma tissue and normal skin tissue samples. When λ was set as 0.04, 20 feature CRGs were selected and the LASSO model constructed using these CRGs could accurately distinguish melanoma from normal skin (Figure 4C). Finally, by comparing the feature genes identified by the three machine learning approaches, six common key CRGs were identified as melanoma-related key CRGs, including ABCC2, CA14, EGR3, FBXW7, LDHB, and PSEN2 (Figure 4D).

Expression and ROC analysis of key CRGs

The expression of key CRGs was examined in the training set (combined GSE15605 and GSE114445 datasets) and validation sets (GSE15605, GSE114445 and GSE46517 datasets) to investigate their expression in melanoma versus normal skin tissues, as well as evaluate their diagnostic value as melanoma biomarkers. As shown in Figure 5A–5D, compared to normal skin tissues, 4 key CRGs (ABCC2, CA14, LDHB, PSEN2) were significantly upregulated while 2 key CRGs (EGR3, FBXW7) were significantly downregulated in melanoma tissues. ROC analysis demonstrated that in each dataset, all key CRGs had an AUC>0.7, indicating these key CRGs have good ability to distinguish melanoma from normal skin tissues (Figure 5E–5H).

Prognostic analysis of key CRGs

Further prognostic analysis of key CRGs was performed based on the TCGA dataset. Kaplan-Meier survival curve analysis showed all key CRGs were closely associated with prognosis in melanoma patients. Specifically, high expression of 4 key CRGs (ABCC2, CA14, LDHB, PSEN2) was significantly correlated with poor prognosis in melanoma patients (HR>1, P<0.05), while low expression of 2 key CRGs (EGR3, FBXW7) was significantly associated with poor prognosis (HR<1, P<0.05) (Figure 6A–6F).

Subsequently, multivariate Cox regression analysis was conducted on the 6 prognosis-related key CRGs in the TCGA-melanoma cohort to construct a prognostic risk model based on their weighted regression coefficients. The risk score for each patient was calculated as follows: Risk Score = (0.0320×ABCC2) + (0.0005×CA14) + (-0.0146×EGR3) + (-0.1551×FBXW7) + (0.0010×LDHB) + (0.0073×PSEN2). According to the median value of risk scores, the TCGA cohort was classified into high-risk and low-risk groups. As shown in Figure 7A, the expression levels of the 6 key CRGs were displayed along with the risk score distribution and survival status of all patients. Kaplan-Meier analysis revealed that patients in the high-risk group had significantly poorer prognosis compared to the low-risk group (Figure 7B). The ROC results showed that this prognostic model achieved AUC of 0.72, 0.76 and 0.80 for predicting 3-, 5- and 10-year overall survival of patients, respectively (Figure 7C). Similarly, risk scores were calculated for patients in the validation cohort (GSE65904), and no significant difference in prognosis was found between the high- and low-risk groups (Figure 7D). The risk score model achieved AUC of 0.64, 0.71 and 0.71 for predicting 3-, 5- and 10-year overall survival (Figure 7E). These results indicated this prognostic model had good predictive performance.

We next investigated the correlation between the risk score and clinicopathological features of patients in the TCGA cohort, including age, gender, stage, T classification, N classification and M classification. The risk score exhibited significant positive correlations with advanced T classification (P=0.02) (Figure 8A–8D).

Univariate and multivariate Cox regression analyses integrating clinical variables (age, sex, stage, T, N and M classification) and risk score identified T classification, N classification and risk score as independent prognostic factors (Figure 8E, 8F). Therefore, a nomogram was constructed by integrating T classification, N classification and risk score to predict patients’ survival probability (Figure 8G). The calibration curve verified high prediction accuracy of the nomogram, which achieved AUC of 0.83, 0.84 and 0.88 for predicting 3-, 5- and 10-year overall survival, respectively (Figure 8H). Compared with risk score alone, the nomogram had higher AUC values, indicating that the risk score had better predictive potential when combined with clinical factors.