記住我
ORIGINAL ARTICLE
Year : 2023 | Volume
: 41
| Issue : 2 | Page : 103-110
Transcriptomic profiling and classification of skin melanoma based on ultraviolet response
Dongxing Xiao1, Zhaozhao Guo2, Yuzhen Xiong1, Xinqiang He1, Chong Zhao1, Ni Tang3
1 Department of Plastic Surgery, Longgang Central Hospital of Shenzhen, Shenzhen, Guangdong, China
2 Longyuan Daguan Community Health Service Center, Shenzhen Longgang Central Hospital, Shenzhen, Guangdong, China
3 Department of Dermatology, Longgang Central Hospital, Shenzhen, Guangdong, China
Correspondence Address:
Dr. Dongxing Xiao
Department of Plastic Surgery, Longgang Central Hospital, Shenzhen, No. 6082 Longgang Avenue, Shenzhen, Guangdong Province
China
Source of Support: None, Conflict of Interest: None
DOI: 10.4103/ds.DS-D-22-00178
Background: We aimed to explore the therapeutic biomarker for cutaneous melanoma (CM). Objectives: The objective is to provide a novel direction for improving overall survival (OS) for CM. Methods: We obtained the gene sets related to ultraviolet (UV) reaction from MsigDB database and CM HTSeq-FPKM data from The Cancer Genome Atlas (TCGA). Gene set variation analysis was used to calculate the enrichment scores in each sample. DAVID and Gene Set Enrichment Analysis (GSEA) were used to explore the function of differentially expressed genes (DEGs) between cluster 1 and cluster 2. The ssGSEA was used to analyze the degree of immune infiltration in samples. Weighted gene co-expression network analysis (WGCNA), protein–protein interaction (PPI) network, and mutation analysis were performed to screen the DEGs related to UV response. Results: The samples were divided into the high activity of UV response (cluster 1) and low activity of UV response (cluster 2). We found that cluster 2 was related to poorer OS and had a higher reaction to UV response. Function analysis indicated that the DEGs are involved in angiogenesis, epidermal development, and inflammatory reaction. Furthermore, the cluster 2 had a higher degree of immune infiltration. The results of WGCNA indicated that the genes in the MEyellow module were highly related to UV response, which is involved in the process of angiogenesis, cell migration, and skin development. PPI and mutation analysis indicated that COL5A1 was the risk factor for CM. Conclusion: COL5A1 might be an important biomarker and potential therapeutic target of CM.
Keywords: COL5A1, cutaneous melanoma, immune infiltration, UV Response
#Dongxing Xiao and Zhaozhao Guo contributed equally to this study.
Introduction
Cutaneous melanoma (CM) is a type of the deadliest tumors is still a problem of threatened public health worldwide, especially the populations form European.[1] According to data form GLOBOCAN 2020 database, a total of 325 000 new CM cases (including 151 000 females and 174 000 males) and 57 000 deaths (25 000 females and 32 000 males) were evaluated in 2020.[2] Surgery, target therapy, and immunotherapy are the main therapies for CM.[3],[4],[5] Previous studies have reported that the 5-year overall survival (OS) of CM patients was 97% when at the early stage of CM, while the 5-year OS was only 3% when at the advantage of CM.[6],[7] Thus, the identification novel biomarkers for the early stage of CM is beneficial for precision treatment and improving the OS of patients.
It is well-known that DNA damage induced by ultraviolet (UV) radiation and following mutations are risk factors for CM.[8],[9] In melanocytes, DNA is considered a target, for example, patients with xeroderma pigmentosum have over 1000-fold increased risk of CM induced by the sun.[10] Thus, UV response plays an important role in the procession of CM. As exposure to UV and the incidence of CM continues to increase in the United States, identifying genetic biomarkers related to UV would be beneficial for the diagnosis of at-risk individuals and preventing for patients with CM.[11] Wang et al.[12] found that 4-NQO might be an important biomarker for CM, which has sensitive to UV-induced chromosome breaks.[12] However, only a few biomarkers were explored for CM. It is urgent to identify more biomarkers to improve the prognosis and OS for patients with CM.
Therefore, our study aimed to explore the potential biomarkers for CM using the public database including MsigDB and UCSC Xena database. We first divided the patients into weak UV reactivity (cluster 1) and strong UV reactivity. Then, WGCNA was used to screen the significant modules related to UV response, whose function was further analyzed using DIVIAD. Finally, the construction of the PPI network, survival analysis and mutation analysis were used to identify the risk factors for CM.
Materials and MethodsCollection and preprocessing of the cancer genome atlas-CM
The Htseq-FPKM data of TCGA-CM was downloaded from UCSC Xena database[13] (http://xena.ucsc.edu/). The data was normalized using the expression value of log2 (fpkm + 1). Then, the Ensembl ID was converted to a gene symbol based on the file of gencode.v22.annotation.gene.probeMap. Meanwhile, we obtained the clinical and prognosis data of patients with CM. To exclude the interference of surgery, radiotherapy, and chemotherapy on the survival time of patients, the patients with survival time >90 days were included in this study. Subsequently, a total of 446 patients with CM were screened for this study, the detailed information is listed in [Supplementary Table 1].
Construction of the UV response classifier
We downloaded the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology biological process (GOBP) from MsiDB database[14] (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) and then tailored UV response-related gene sets into the GOBPs, including Hallmark_UV_Response_Up, Response to UV A and Response to UV B. The detailed information was shown in [Supplementary Table 2].
To calculate the enrichment of UV response GOBP in CM samples, we performed the gene set variation analysis (GSVA)[15] with the parameter of method='ssgsea', kcdf='Gaussian', abs. ranking = TRUE. According to the enrichment scores of gene sets, the “ConsensusClusterPlus”[16] in R was used to cluster the samples with different degrees of activity of the UV response. Then, the k-means was used to cluster the samples into 2–6 subtypes, the results indicated that the slope of CDF decline is the smallest when k = 2, suggesting the clustering had the best efficiency [Supplementary Figure 1]. To better understand the differences between subtypes, we used t-SNE dimensionality reduction to analyze the distribution of samples. Furthermore, the Kaplan–Meier (KM) curves were used to analyze the prognosis of cluster 1 and cluster 2. The survival time, survival status, and subtype of each sample are shown in [Supplementary Table 3].
Identification and function annotation of differential expression genes
The limma package[17] was performed to screen the differentially expressed genes (DEGs) between cluster 1 and cluster 2 with the threshold of |logFC|>1 and P < 0.05. After that, the DAVID database was used to analyze the function annotation of the DEGs. Specially, we focused on the biological process with the threshold of P < 0.05. We then performed the GSEA to identify the GOBP and KEGG pathways involved in cluster 1 and cluster 2, respectively. The cutoff set to NOM p-val <0.05. Furthermore, we referred to the 28 immune cells and corresponding 782 marker genes provided by Bindea et al.[18] to analyze the differences in the infiltration degree of each immune cell in cluster 1 and cluster 2 tissues. Thereby, we used ssGSEA to calculate the infiltration degree of 28 types of immune cells in the sample with CM.
Weighted gene co-expression network analysis
It is well know that the gene expression patterns involved in the same pathway or biological process are similar.[19] A co-expression network is constructed using WGCNA that transforming co-expression correlation into topological overlap value or connection weight.[20] According to the background, the WGCNA was performed to establish a co-expression gene network associated with cluster 1 and cluster 2. According to the scale-free network standard, we used the “pickSoftThreshold” function of the WGCNA package to select an appropriate soft threshold β (ranging from 1 to 20). We obtained the best efficiency when β was equal to 9. Then, the blockwiseModules function was performed for one-step network construction and module detection. Finally, we calculated the correlation between module eigengene (ME) and sample clustering. We focused on the MEyellow module highly correlated to cluster 2, the genes of which were selected for the further analysis.
Construction of protein-protein interaction network
The String database[21] (https://cn. string-db. org/) was used to analyze the interaction of genes in MEyellow module at the protein level. The interaction score set to high confidence (0.700). After that, the MCODE plugin of Cytoscape[22] was used for cluster analysis to explore the tightly connected protein nodes. Moreover, the simple nucleotide variation data were downloaded from TCGA (https://portal. gdc. cancer. gov/) to analyze the mutation of the genes of MEyellow in TCGA database. The oncoplot function in maftools package[23] was used to exhibit the top 20 mutation gene types with the highest mutation frequency of MEyellow in TCGA.
Statistical analysis
We used the Wilcoxon rank-sum test to compare the difference in continuous variables between cluster 1 and cluster 2. The continuous-adjusted Chi-square test was performed to compare the difference of categorical variables between the two groups, and the log-rank test to compare the difference in survival time. Univariate cox regression was used to calculate the association of genes with survival time. All statistical analyses were performed using the R language (version 4.1.2), and a P < 0.05 was considered to be statistically significant.
ResultsHigh UV reactivity related to poor prognosis of patients with CM
First, all the patients with CM were divided into two groups including cluster 1 (weak) and cluster 2 (strong) through the ssGSEA enrichment scores of gene sets related to UV response [Figure 1]a. We further found the obvious difference in the biological process of UV response between cluster 1 and cluster 2 through t-SNE dimensionality reduction analysis [Figure 1]b. Besides, to analyze the effect of UV response on the prognosis of patients with CM, KM curves indicated that the CM patients in cluster 1 had a better prognosis than that in cluster 2 (P < 0.009, [Figure 1]c). Interestingly, the ssGAEA scores of gene sets related to UV response in cluster 2 were higher than cluster 1 [Figure 1]d, suggesting that the patients in cluster 2 had a stronger reaction to UV response. Finally, we further explore the expression of gene sets of response to UV A and response to UV B, the results indicated that the expression of genes related to UV response in cluster 2 was higher than cluster 1 [Figure 1]e and [Figure 1]f.
Figure 1: High UV reactivity related to poor prognosis of patients with CM.(a) Consistent clustering heatmap at k = 2; (b) Sample distribution after t-SNE dimensionality reduction; (c) Survival curves of cluster 1 and cluster 2; (d) Boxplot of the ssGSEA score distribution of the UV-response-related gene set for cluster 1 and cluster 2.(e) Boxplot of GOBP_RESPONSE_TO_UV_A gene expression for cluster 1 and cluster 2.(f) GOBP_RESPONSE_TO_UV_B gene for cluster 1 and cluster 2 Expressed boxplot. CM, cutaneous melanoma; GOBP, Gene Ontology biological process.High UV reactivity related to high pathologic T staging
We analyzed the difference of clinical information between cluster 1 and cluster 2, indicating that patients in cluster 2 had a higher degree of pathologic T staging. The other factors including age, gender, stage, and pathologic N and M were not correlated to UV response clustering [Figure 2].
Figure 2: Tumor tissues with high UV reactivity had a higher degree of T staging. (a) Difference in age between cluster 1 and cluster 2. (b) Difference in gender between cluster 1 and cluster 2. (c) Difference between cluster 1 and cluster 2 in staging. (d) Difference between cluster 1 and cluster 2 in T staging. (e) Difference between cluster 1 and cluster 2 in pathologic N staging. (f) Differences between cluster 1 and cluster 2 in pathologic M staging.High UV reactivity related to high degree of infiltration of immune cells in CM tissues
As shown in [Figure 3]a, we screened the DEGs between cluster 1 and cluster 2, and then analyzed the biological function of ones. Interestingly, we found that the up-regulated DEGs in cluster 2 were mainly involved in extracellular matrix organization, angiogenesis, epidermis development, inflammatory response, cell migration, wound healing [Figure 3]b, while the down-regulated DEGs participated in negative regulation apoptotic process, sialylation, and positive regulation of bone resorption [Figure 3]c. The results suggested that high UV reactivity were associated with angiogenesis, epidermal development, and inflammatory response in CM. Meanwhile, we also found that the enhancement of UV response activity has the effect of stimulating immune cell infiltration. We then further performed GSEA analysis to validate the above results. The results indicated that the DEGs in cluster 2 were enriched in keratincati on, leukocyte migration involved in inflammatory response, and angiogenesis involved in wound healing [Figure 3]d, [Figure 3]e, [Figure 3]f.
Figure 3: Immune cell infiltration between cluster 1 and cluster 2. (a) Volcano plot of differentially expressed genes between cluster 2 and cluster 1; (b) Biological process involved in up-regulated genes in cluster 2; (c) Biological process involved in down-regulated genes in cluster 2; (d-f) GSEA analysis between cluster 2 and cluster 1; (g) Boxplot of 28 immune cell ssGSEA enrichment scores between cluster 2 and cluster 1.Furthermore, we calculated the infiltration ratio of 28 types of immune cells in both cluster 1 and cluster 2, respectively using ssGSEA. The results indicated that the ssGSEA scores of immune cells (central memory CD4 T cell, CD56 dim natural killer T cell, central memory CD8 T cell, activated dendritic cell, and regulatory T cell) in cluster 2 were higher than cluster 1 [Figure 3]g. The results further validated that the stronger activity of UV response in CM, the higher degree infiltration of various immune cells was.
Identification of modules and genes related to UV response
We used WGCNA analysis to identify gene modules associated with UV-response classification and found that the network best fit the scale-free network characteristics when the soft threshold β = 9 [Figure 4]a and [Figure 4]b. Then, we use a dynamic hybrid cut method to build a hierarchical clustering tree. Each leaf on the tree represented a gene, and each branch represented a co-expression module in CM [Figure 4]c. We exhibited the correlation coefficient between UV-response classification (cluster 1 and cluster 2) and co-expression gene module. We found that the MEyellow had highly correlation with cluster 2 (R = 0.3, P = 1e-10, [Figure 4]d). Meanwhile, we also found that there is correlation between module membership in the yellow module and gene significance for cluster 2 [Figure 4]e.
Figure 4: The identification of relevant modules and genes for UV-response classification. (a) Screening of soft thresholds; (b) Relationship between soft thresholds and average connectivity; (c) Hierarchical clustering tree was constructed with dynamic mixed cutting method; (d) Correlation between UV-responsive classification and co-expressed gene module; (e) Correlation between Module membership in yellow module and Gene significance for cluster 2.Identification of key genes in MEyellow module
As shown in [Figure 5]a, we found that the genes in the MEyellow module were mainly involved in angiogenesis, extracellular matrix organization, cell migration and skin development. Then, we constructed the PPI network based on the genes in the MEyellow module, indicating that COL5A1, COL5A2, and ADAMTSL2 have the most edges [Figure 5]b. Moreover, we screened the DEGs related to prognosis in MEyellow indicating a total of 13 genes (such as COL5A1, COL1A1, and MEDAG) were closely correlated to the prognosis of patients with CM [Figure 5]c. Finally, we further explored the mutation frequency of genes in the samples with CM. We exhibited the top 20 mutation genes indicating that COL5A1 was not only the mutation gene with a 19% mutation rate but also the prognosis risk factor in CM (HR = 1.11, P = 0.01733) [Figure 5]d. The results suggested that COL5A1 played an important role in the activation of UV response in CM.
Figure 5: The identification of key genes in MEyellow. (a) Biological process of enrichment of all genes in MEyellow; (b) Protein–protein interaction network constructed by all genes in MEyellow; (c) Forest plot of prognosis-related factors in MEyellow module; (d) Waterfall plot of the top 20 genes with the highest mutation frequency in the MEyellow module. Discussion CM is the most serious disease of skin cancer with the feature of low 5-year OS. UV has been reported to be the risk factor for CM. Thus, seeking reliable biomarkers related to UV is a promising therapy for predicting the prognosis of patients with CM and improving their life quality. In our study, we obtained the data from MsigDB and UCSC Xena database, and then screened the gene sets related to UV response. Then, we divided the gene sets into weak UV reactivity (cluster 1) and strong UV reactivity based on ssGSEA scores. The KM curves indicated that strong UV reactivity was correlated to a poor prognosis of CM. The results of WGCNA found that the MEyellow module has the highest correlation with UV response. Further analysis composed of PPI network, survival, and mutation analysis represented that COL5A1 played an important role in the activation of UV response in CM.
In our study, we found that the patients with CM in cluster 2 had poorer OS than cluster 1. Further analysis indicated that patients in cluster 2 had higher T stage and degree of infiltration of immune cells. These results were consistent with the previous studies. Studies have reported that excessive exposure of skin resulted in inflammation, oxidative stress, and DNA damage and mutations, which leaded to increasing incidence of CM and reducing the OS of patients with CM.[24] Besides, Agnieszka et al.[25] indicated that the UV response induced a high degree of infiltration of immune cells and exacerbated the situation of CM. All the results suggested that UV response might increase the degree of immune infiltration and T stage, and then shorten the OS of patients with CM. Further analysis indicated that the up-regulated genes in cluster 2 were involved in extracellular matrix organization, angiogenesis, epidermis development, inflammatory response, cell migration, and wound healing. These functions were related to characteristics of malignant tumors. Thus, we inferred that the genes in cluster 2 might mediate the procession of CM through the above biological processes.
Moreover, we further analyzed the genes in cluster 2 using WGCNA and the PPI network. The results indicated that the genes significantly related to UV response were enriched in the MEyellow module. Furthermore, we found these genes in the MEyellow module were participated in angiogenesis, extracellular matrix organization, cell migration, and skin development. The results suggested that these genes might be important biomarkers in the development of CM. Furthermore, composing of the PPI network and mutation analysis, we found that COL5A1 was the prognosis risk factor in CM with a mutation rate of 19%. The above results suggested that COL5A1 played an important role in the activation of CM response to UV radiation. COL5A1 encodes an alpha chain for one of the low-abundance fibrillar collagens. Previous studies have reported that COL5A1 was considered an important biomarker for predicting the prognosis of various cancers, such as gliomas,[26] ovarian cancer,[27] and papillary thyroid carcinoma.[28] However, few studies have reported the function of COL5A1 in the development of CM. Our study provided a novel sight for exploring the therapies of patients with CM.
Our study still has some limitations. The results of our study need to be validated using an integrated bioinformatics analysis using a large size of samples. Then, functional experiments on COL5A1 were necessary to be performed to improve the meaning of our study. Finally, detailed mechanism of the GO BP involved in the genes related to UV response also needs to be explore in the further.
ConclusionsIn summary, we found that the strong activity of UV response was correlated to poor OS, high degree of immune infiltration, and T stage. COL5A1 might be an important biomarker and potential therapeutic target of CM.
Data availability statement
The datasets used and/or analyzed during the current study are available from the corresponding author via email request.
Financial support and sponsorship
Nil.
Conflicts of interest
There are no conflicts of interest.
References
留言 (0)